Hatree Fock超详细推导 您所在的位置:网站首页 轨道能量公式推导 Hatree Fock超详细推导

Hatree Fock超详细推导

2024-06-14 10:58| 来源: 网络整理| 查看: 265

第三章 Hartree-Fock theory

目录第三章 Hartree-Fock theory3.1 序言3.2 绝热近似与BO近似3.3 HF理论3.4 Slater行列式的能量3.5 Koopmans' Theorem3.6 Roothaan-Hall 方程3.7 HF程序

3.1 序言

独立粒子近似:忽略电子之间的作用或采用平均场

HF理论:全波函数由电子波函数作为基底展开。由于电子为费米子,全波函数必须满足反对称性,通过引入Slater行列式解决。通过变分原理,在波函数正交归一化条件的限制下,求解泛函的极值得HF方程,迭代求解即可得波函数的解。

3.2 绝热近似与BO近似

体系的哈密顿算符为:

\[H_{total} = T_n+H_e+H_{mp}=T_n+(T_e+V_{ne}+V_{ee}+V_{nn})+H_{mp} \]

其中\(H_e\)即电子的哈密顿算符,\(H_{mp}\)称作质量极化(当体系中有超过两个粒子时,体系的质心运动和内部运动无法严格分离,故引入该项)

\[H_{mp}=-\frac{1}{2M_{tot}}(\sum_i^{N_{elec}}\nabla_i)^2 \]

我们假设电子的薛定谔方程为:(注意对于电子来说,原子核的位置\(R\)可视为常数)

\[H_e(R)\Psi_i(R,r)=E_i(R)\Psi_i(R,r),i=0,1,2\cdots n \]

上述\(n\)个方程对应\(n\)个本征态

物理量算符都是厄米算符,故对于哈密顿算符有:

\[\int\Psi_i^*H\Psi dx=\int \Psi_jH^*\Psi_i^*dx\\ \rightarrow \langle\Psi_i|H|\Psi_j\rangle=\langle \Psi_j|H|\Psi_i\rangle^* \]

同一算符对应不同本征值的本征态是相互正交的:(注意这里的\(i\)和\(j\)对应的是某一个本征态)

\[\int\Psi_i(R,r)^*\Psi_j(R,r)dr=\delta_{ij}=\begin{cases}1,&i=j\\0,&i\neq j \\\end{cases} \]

假设体系的总波函数可以展开为不同本征态下电子全波函数和原子核位置函数的线性组合:

\[\Psi_{tot}(R,r)=\sum_{i=1}^{\infty}\Psi_{ni}(R)\Psi_i(R,r) \]

代入体系的薛定谔方程有:

\[\sum_{i=1}^{\infty}(T_n+H_e+H_{mp})\Psi_{ni}(R)\Psi_i(R,r)=E_{tot}\sum_{i=1}^{\infty}\Psi_{ni}(R)\Psi_i(R,r) \]

因为:

\[T_n=-\sum_A\frac{1}{2M_A}\nabla^2_A=\nabla^2_n \]

代入上述方程:

\[\sum_{i=1}^{\infty}(\nabla^2_n+H_e+H_{mp})\Psi_{ni}\Psi_i=E_{tot}\sum_{i=1}^{\infty}\Psi_{ni}\Psi_i \]

我们认为\(H_e\)和\(H_{mp}\)只作用在电子波函数上,同时\(\Psi_i\)是对应\(i\)本征态下电子薛定谔方程的准确解

即\(H_e\Psi_i=E_i\Psi_i\),则方程可变为:

\[\sum_{i=1}^{\infty}\bigg\{\nabla_n^2\Psi_{ni}\Psi_i+\Psi_{ni}H_e\Psi_i+\Psi_{ni}H_{mp}\Psi_i\bigg\}=E_{tot}\sum_{i=1}^{\infty}\Psi_{ni}\Psi_i\\ \sum_{i=1}^{\infty}\bigg\{\nabla_n(\Psi_{ni}\nabla_n\Psi_i+\Psi_i\nabla_n\Psi_{ni})+\Psi_{ni}E_i\Psi_i+\Psi_{ni}H_{mp}\Psi_i\bigg\}=E_{tot}\sum_{i=1}^{\infty}\Psi_{ni}\Psi_i\\ \sum_{i=1}^{\infty}\bigg\{\Psi_{ni}(\nabla_n^2\Psi_i)+2(\nabla_n\Psi_i)(\nabla_n\Psi_{ni})+\Psi_i(\nabla_n^2\Psi_{ni})+\Psi_{ni}E_i\Psi_i+\Psi_{ni}H_{mp}\Psi_i\bigg\}=E_{tot}\sum_{i=1}^{\infty}\Psi_{ni}\Psi_i \]

根据不同本征态下的本征波函数相互正交的性质,我们在方程两边分别左乘\(\Psi_j^*\),并在全空间积分:

上述方程左边括号中的第三项和第四项以及方程右边仅在\(i=j\)时不为零

\[\nabla_n^2\Psi_{nj}+E_j\Psi_{nj}+\sum_{i=1}^{\infty}\bigg\{2\langle\Psi_j|\nabla_n|\Psi_i\rangle(\nabla_n\Psi_{ni})+\langle\Psi_j|\nabla_n^2|\Psi_i\rangle\Psi_{ni}+\langle\Psi_j|H_{mp}|\Psi_i\rangle\Psi_{ni}\bigg\}=E_{tot}\Psi_{nj} \]

上述括号中的即是不同电子态的耦合项,其中前两项分别为一阶和二阶非绝热耦合项,第三项为质量极化。

非绝热耦合项也就是说体系中会涉及不同电子态之间的作用(不同electronic surface之间存在耦合)

C6E7050C32FA144B099DFF40D1F86EBE

这里我们引入绝热近似,即限制体系仅在一个电子态上,即只剩\(i=j\)的项

这里还说除了空间简并的波函数(对应同一本征值的多个线性无关的本征函数),一阶非绝热项的对角元素均为零,意思就是忽略掉该项,为何???

再忽略掉质量极化项,方程变为:

\[(\nabla_n^2+E_j+\langle\Psi_j|\nabla_n^2|\Psi_j\rangle)\Psi_{nj}=E_{tot}\Psi_{nj}\\ (T_n+E_j(R)+U(R))\Psi_{nj}(R)=E_{tot}\Psi_{nj}(R) \]

其中上述的\(U(R)\)称为对角修正项,远小于\(E_j(R)\),约是电子和原子的质量比。我们可看到此时方程即原子的薛定谔方程,因为原子的势能面基本上由电子能量决定

注意此时要区分一下绝热近似和BO近似:

(1)BO近似认为原子的势能面(PES)与对角修正项无关,即忽略掉\(U(R)\)项,因而方程可以进一步简化为:

\[(T_n+E_j(R))\Psi_{nj}(R)=E_{tot}\Psi_{nj}(R) \]

(2)绝热近似则认为对角修正项不可忽略

(3)当求解电子薛定谔方程,存在多个解较为接近时,BO近似就不适用了

(4)对于大多数体系,BO近似误差很小,我们可用diagonal Born-Oppernherimer correction(DBOC)参数衡量:(小原子影响大,大原子影响小)

\[\Delta E_{DBOC}=\sum_{A=1}^{N_{nuc}}-\frac{1}{2M_A}\langle\Psi_e|\nabla_A^2|\Psi_e\rangle \]

(5)当绝热近似不成立时,我们要重新回到体系总波函数的位置,此时要分别用高斯函数作为基底函数构成的行列式波函数来描述原子核和电子的波函数,再经过一定程度的近似,我们称这些近似方程为Nuclear Orbital plus Molecular Orbital(NOMO),最后求解得到的能谱就能包括原子核和电子的状态了。

3.3 HF理论

前面我们提到要求解电子的薛定谔方程得到电子的能量\(E_j\),因此引入HF近似

我们还要考虑电子的自旋,其中自旋函数为\(\alpha、\beta\),满足正交归一化条件:

\[\langle\alpha|\alpha\rangle=\langle \beta|\beta\rangle=1\\ \langle\alpha|\beta\rangle=\langle \beta|\alpha\rangle=0 \]

要得到电子能量的近似解,要采用变分原理。由于波函数的近似解对应的能量高于实际解的能量,因为我们可以认为当能量最低时,对应的为解趋近于精确解。

近似波函数对应的能量可表示为哈密顿算符的期望值和波函数归一化积分的比值:

\[E_e=\frac{\langle\Psi|H_e|\Psi\rangle}{\langle\Psi|\Psi\rangle} \]

若波函数已经归一化,则归一化系数为1,得:

\[E_e=\langle\Psi|H_e|\Psi\rangle \]

补充知识:

全同粒子:质量、电荷和自旋等性质完全相同的粒子。在量子力学中,粒子的状态是用波函数描述的,当两个粒子的波函数在空间重叠时,无法区分,所以有全同粒子的不可区别性原理。

自旋:微观粒子的一种性质。自旋为0的粒子从各个方向看都一样,自旋为1的粒子旋转360°后一样,自旋为2的粒子旋转180°后一样,自旋为1/2的粒子旋转2圈后一样。自旋为半整数的费米子服从泡利不相容原理,玻色子不遵从。

对全同粒子体系波函数引入交换算符\(\hat{P}_{ij}\),即交换两个粒子的位置,有:

\[\hat{P}_{ij}\Psi(\cdots,q_i,\cdots,q_j,\cdots)=\Psi(\cdots,q_j,\cdots,q_i,\cdots) \quad (j\neq i) \]

由全同粒子的不可区别性,可知交换后的状态与原来的状态是不可区分的(即交换后波函数与原波函数线性相关),因此有:

\[\hat{P}_{ij}\Psi=C\Psi \]

而:

\[\hat{P}_{ij}\hat{P}_{ij}\Psi=\hat{P}_{ij}(C\Psi)=C\hat{P}_{ij}\Psi=C^2\Psi=\Psi \]

所以有:

\[C^2=1\\ \therefore C=±1 \]

即:

\[\hat{P}_{ij}\Psi=+\Psi \]

称为交换对称波函数

或:

\[\hat{P}_{ij}\Psi=-\Psi \]

称为交换反对称波函数

交换对称性与粒子的自旋有确定的关系,自旋为半整数的粒子,波函数是交换反对称的,称为费米子。自旋为整数的粒子,波函数是交换对称的,称为玻色子。而电子是属于费米子。

全同粒子的波函数是求解薛定谔方程得到的,原始解未必具有确定的交换对称性,需要进行对称化或反对称化。

电子的坐标可以用空间坐标和自旋坐标表示:

\[q_i=\{r_i,\omega_i\} \]

考虑无耦合体系,即体系电子总波函数是单粒子波函数的乘积:(独立粒子近似)

\[\Psi_k(q_1,\cdots,q_N)=\phi_1(q_1)\cdots\phi_N(q_N) \]

以二粒子为例,未进行对称化的波函数为:

\[\Psi_k(q_1,q_2)=\phi_1(q_1)\phi_2(q_2) \]

若两个单粒子波函数为不同的函数,则对称化的波函数是:

\[\Psi_k(q_1,q_2)=\frac{1}{\sqrt{2}}[\phi_1(q_1)\phi_2(q_2)+\phi_1(q_2)\phi_2(q_1)] \]

反对称化的波函数是:

\[\Psi_k(q_1,q_2)=\frac{1}{\sqrt{2}}[\phi_1(q_1)\phi_2(q_2)-\phi_1(q_2)\phi_2(q_1)] \]

对于\(N\)个电子的体系,对应的反对称化波函数可用Slater行列式表示:

\[\Psi_k(q_1,\cdots,q_N)=\frac{1}{N!}\begin{vmatrix}\phi_{k1}(q_1)&\phi_{k2}(q_1)&\cdots&\phi_{kN}(q_1)\\ \phi_{k1}(q_2)&\phi_{k2}(q_2)&\cdots&\phi_{kN}(q_2)\\ \vdots&\vdots&\ddots&\vdots\\ \phi_{k1}(q_N)&\phi_{k2}(q_N)&\cdots&\phi_{kN}(q_N)\end{vmatrix} \]

当存在任何两个电子的位置和自旋相同时(\(q_i=q_j\)),行列式的值为0,因此有泡利不相容原理:不可能由两个或更多的电子处于完全相同的量子态中。

同时电子的不同自旋轨道相互正交:

\[\langle\phi_i|\phi_j\rangle=\delta_{ij} \]

3.4 Slater行列式的能量

我们定义一个反对称算符\(A\),取Slater行列式的对角项乘积的总函数为\(\Pi\),将\(A\)作用于\(\Pi\)表示为所有的置换结果的总和,则可以将电子全波函数的Slater行列式表示为:

\[\Psi_k=A[\phi_{k1}(q_1)\phi_{k2}(q_2)\cdots\phi_{kN}(q_N)]=A\Pi \]

其中\(A\)为:

\[A=\frac{1}{\sqrt{N!}}\sum_{p=0}^{N-1}(-1)^pP=\frac{1}{\sqrt{N!}}\bigg\{1-\sum_{ij}P_{ij}+\sum_{ijk}P_{ijk}-\cdots\bigg\} \]

根据置换的轨道的个数共可分为\(N\)种情况,总共得到\(N!\)个波函数

\(A\)满足的性质:

(1)\(A\)和\(H\)两算符对易,\(AH=HA\)

(2)\(AA=\sqrt{N!}A\),证明如下:

\[A^2=\frac{1}{N!}\sum_{p=0}^{N-1}\sum_{q=0}^{N-1}(-1)^p(-1)^qPQ \]

上述两个变换叠加虽然进行了\(N!N!\)次变换,但是仍然仅有\(N!\)个结果。将两步变换合并为一个变换\(R\),仍要保留次数为\(N!N!\),因此可进一步推得:

\[A^2=\frac{1}{N!}\sum_{R}^{N-1}(-1)^RR\sum_{P}^{N-1}1=\frac{1}{N!}\sum_{R}^{N-1}(-1)^RRN!=\sum_R^{N-1}(-1)^RR=\sqrt{N!}A \]

前面定义的电子的哈密顿算符为:

\[H_e=T_e+V_{ne}+V_{ee}+V_{nn} \]

分别为:

\[T_e=-\sum_i^{N_{elec}}\frac{1}{2}\nabla_i^2\\ V_{ne}=-\sum_{A}^{N_{nuc}}\sum_i^{N_{elec}}\frac{Z_A}{|R_A-r_i|}\\ V_{ee}=\sum_i^{N_{elec}}\sum_{j>i}^{N_{elec}}\frac{1}{|r_i-r_j|}\\ V_{nn}=\sum_A^{N_{nuc}}\sum_A^{N_{nuc}}\frac{Z_Az_B}{|R_A-R_B|} \]

其中一、二项为单电子项,第三项为双电子项,故令:

\[h_i=-\frac{1}{2}\nabla_i^2-\sum_A^{N_{elec}}\frac{Z_A}{|R_A-r_i|}\\ g_{ij}=\frac{1}{|r_i-r_j|}\\ \Rightarrow H_e=\sum_i^{N_{elec}}h_i+\sum_i^{N_{elec}}\sum_{j>i}^{N_{elec}}g_{ij}+V_{nn} \]

波函数的能量可进一步推导:(哈密顿算符\(H\)和置换算符\(A\)对易)

\[E=\langle\Psi|H|\Psi\rangle=\langle A\Pi|H|A\Pi\rangle=\langle\Pi|AHA|\Pi\rangle=\langle\Pi|HAA|\Pi\rangle\\ =\sqrt{N!}\langle\Pi|H|A\Pi\rangle=\sum_p(-1)^p\langle\Pi|H|P\Pi\rangle \]

分三项考虑:

(1)原子核-原子核排斥项:

\[\langle\Psi|V_{nn}|\Psi\rangle=V_{nn}\langle\Psi|\Psi\rangle=V_{nn} \]

(2)单电子项:

无轨道交换:

\[\langle\Pi|h_i|\Pi\rangle=\langle\phi_1(q_1)\phi_2(q_2)\cdots\phi_N(q_N)|h_i|\phi_1(q_1)\phi_2(q_2)\cdots\phi_N(q_N)\rangle\\=\langle\phi_1(q_1)|\phi_1(q_1)\rangle\cdots\langle\phi_i(q_i)|h_i|\phi_i(q_i)\rangle\cdots\langle\phi_N(q_N)|\phi_N(q_N)\rangle\\ =\langle\phi_i(q_i)|h_i|\phi_i(q_i)\rangle=\bar{h}_i \]

有轨道交换:

\[\langle\Pi|h_i|P_{ij}\Pi\rangle=\\\langle\phi_1(q_1)\cdots\phi_i(q_i)\cdots\phi_j(q_j)\cdots\phi_N(q_N)|h_i|\phi_1(q_1)\cdots\phi_j(q_i)\cdots\phi_i(q_j)\cdots\phi_N(q_N)\rangle\\=\langle\phi_1(q_1)|\phi_1(q_1)\rangle\cdots\langle\phi_i(q_i)|h_i|\phi_j(q_i)\rangle\cdots\langle\phi_j(q_j)|\phi_i(q_j)\rangle\cdots\langle\phi_N(q_N)|\phi_N(q_N)\rangle\\ \because \langle\phi_j(q_j)|\phi_i(q_j)\rangle=0\\ \therefore \langle\Pi|h_i|P_{ij}\Pi\rangle=0 \]

综上,仅当无轨道交换时,单电子项才存在

(3)双电子项:

同理,对于双电子项,如果发生三个或以上的电子轨道的交换,会导致不同的电子轨道重叠积分为零,因而只有不发生轨道交换以及发生轨道交换的恰好是\(g_{ij}\)算符对应的\(i、j\)轨道时,才不为零

无轨道交换:

\[\langle\Pi|g_{ij}|\Pi\rangle=\langle\phi_1(q_1)\phi_2(q_2)\cdots\phi_N(q_N)|g_{ij}|\phi_1(q_1)\phi_2(q_2)\cdots\phi_N(q_N)\rangle\\ =\langle\phi_1(q_1)|\phi_1(q_1)\rangle\cdots\langle\phi_i(q_i)\phi_j(q_j)|g_{ij}|\phi_i(q_i)\phi_j(q_j)\rangle\cdots\langle\phi_N(q_N)|\phi_N(q_N)\rangle\\ =\langle\phi_i(q_i)\phi_j(q_j)|g_{ij}|\phi_i(q_i)\phi_j(q_j)\rangle\\ =J_{ij} \]

则\(J_{ij}\)项被称为库伦积分

\(i、j\)轨道交换:

\[\langle\Pi|g_{ij}|P_{ij}\Pi\rangle=\\\langle\phi_1(q_1)\cdots\phi_i(q_i)\cdots\phi_j(q_j)\cdots\phi_N(q_N)|h_i|\phi_1(q_1)\cdots\phi_j(q_i)\cdots\phi_i(q_j)\cdots\phi_N(q_N)\rangle\\ =\langle\phi_1(q_1)|\phi_1(q_1)\rangle\cdots\langle\phi_i(q_i)\phi_j(q_j)|g_{ij}|\phi_j(q_i)\phi_i(q_j)\rangle\cdots\langle\phi_N(q_N)|\phi_N(q_N)\rangle\\ =\langle\phi_i(q_i)\phi_j(q_j)|g_{ij}|\phi_j(q_i)\phi_i(q_j)\rangle\\ =K_{ij} \]

则\(K_{ij}\)项被称为交换积分(注意由于仅交换了两个轨道,所以前面的系数为-1)

对于所有的电子对,无论轨道的自旋方向如何,库伦积分均非零;但对于交换积分,只有当\(i、j\)两个轨道自旋相同时,该项才不会零。不是很理解

综上,最终能量可表达为:

\[E=\sum_{i=1}^{N_{elec}}\langle\phi_i(q_i)|h_i|\phi_i(q_i)\rangle+\\ \sum_i^{N_{elec}}\sum_{j>i}^{N_{elec}}(\langle\phi_i(q_i)\phi_j(q_j)|g_{ij}|\phi_i(q_i)\phi_j(q_j)\rangle-\langle\phi_i(q_i)\phi_j(q_j)|g_{ij}|\phi_j(q_i)\phi_i(q_j)\rangle)+V_{nn}\\ =\sum_{i=1}^{N_{elec}}\bar{h}_i+\sum_{i=1}^{N_{elec}}\sum_{j>i}^{N_{elec}}(J_{ij}-K_{ij})+V_{nn}\\ =\sum_{i=1}^{N_{elec}}\bar{h}_i+\frac{1}{2}\sum_{i=1}^{N_{elec}}\sum_{j=1}^{N_{elec}}(J_{ij}-K_{ij})+V_{nn} \]

注意上面提出了系数\(1/2\),在没提出前,交换项仅包括来自不同轨道的电子对的贡献,而提出系数后,交换项也包括了相同轨道的自作用项。但是该自作用项会被库伦积分项抵消,因为此时\(J_{ii}\)和\(K_{ii}\)相等

进一步,我们引入库伦算符\(J\)和交换算符\(K\):

\[J_i|\phi_j(q_j)\rangle=\langle\phi_i(q_i)|g_{ij}|\phi_i(q_i)\rangle|\phi_j(q_j)\rangle\\ K_i|\phi_j(q_j)\rangle=\langle\phi_i(q_i)|g_{ij}|\phi_j(q_i)\rangle|\phi_i(q_j)\rangle \\ E=\sum_i^{N_{elec}}\langle\phi_i|h_i|\phi_i\rangle+\frac{1}{2}\sum_{ij}^{N_{elec}}(\langle\phi_j|J_i|\phi_j\rangle-\langle\phi_j|K_i|\phi_j\rangle)+V_{nn}\\ \]

注意上面定义交换算符\(K\)时,右矢并非\(j\)轨道,而是\(i\)轨道,只是为了形式统一用了上述记法,可以看右矢的轨道虽然变为了\(i\)轨道,但是电子是\(j\)电子,根据这个记忆。

此外,算符也可以写成右矢为\(i\)轨道,如\(J_i|\phi_j(q_j)\rangle\)右矢为\(j\)轨道,而实际上库伦积分中左矢和右矢中\(i、j\)轨道在乘积时是可以交换的,故也可定义\(J_j|\phi_i(q_i)\rangle\),则交换积分就变为\(\langle\phi_i|J_j|\phi_i\rangle\),与\(\langle\phi_j|J_i|\phi_j\rangle\)是相等的。

同样对交换积分项也是一样的,但要注意交换积分定义为了形式采用了上述记法,不要弄混了。\(K_j\)算符定义为:\(K_j|\phi_i(q_i)=\langle\phi_j(q_j)|g_{ij}|\phi_i(q_j)\rangle|\phi_j(q_i)\rangle\),对应的交换积分为:\(\langle\phi_i|K_j|\phi_i\rangle\)

以上几点在后面求双电子项能量微分时用到。

接下来的目标是用变分原理找到一组分子轨道使得能量最低。采用变分时,限制条件是分子轨道必须保证正交归一化,通过拉格朗日乘子法来处理。当泛函(拉格朗日函数)关于分子轨道微小变化不变时,取到极值,得到的解即确解。

拉格朗日函数:

\[L=E-\sum_{i}^{N_{elec}}\sum_{j}^{N_{elec}}\lambda_{ij}(\langle\phi_i|\phi_j\rangle-\delta_{ij}) \]

求极值即求:

\[\frac{\delta L}{\delta\phi_{i(j)}}=0\\ \Rightarrow \delta L=L(\phi_i+\delta\phi_i)-L(\phi_i)=0 \]

分几项来看:

(1)单电子项:(忽略微量项)

\[\delta E_1=\sum_i^{N_{elec}}\langle(\phi_i+\delta\phi_i)|h_i|(\phi_i+\delta\phi_i)\rangle-\sum_i^{N_{elec}}\langle\phi_i|h_i|\phi_i\rangle\\ =\sum_i^{N_{elec}}\bigg(\langle\delta\phi_i|h_i|\phi_i\rangle+\langle\phi_i|h_i|\delta\phi_i\rangle\bigg)\\ =\sum_i^{N_{elec}}\bigg(\langle\delta\phi_i|h_i|\phi_i\rangle+\langle\delta\phi_i|h_i|\phi_i\rangle^*\bigg) \]

其中第二项为共轭项

(2)双电子项:(这里要分别对\(i、j\)轨道,下面仅对\(i\)轨道求,\(j\)轨道直接把\(i\)的结果换成\(j\)就行)

\[\delta E_2=\frac{1}{2}\sum_i^{N_{elec}}\sum_j^{N_{elec}}\Bigg(\bigg\langle\bigg(\phi_i(q_i)+\delta\phi_i(q_i)\bigg)\phi_j(q_j)\bigg|g_{ij}\bigg|\bigg(\phi_i(q_i)+\delta\phi_i(q_i)\bigg)\phi_j(q_j)\bigg\rangle\\ -\bigg\langle\bigg(\phi_i(q_i)+\delta\phi_i(q_i)\bigg)\phi_j(q_j)\bigg|g_{ij}\bigg|\phi_j(q_i)\bigg(\phi_i(q_j)+\delta\phi_i(q_j)\bigg)\bigg\rangle \Bigg)\\ -\frac{1}{2}\sum_i^{N_{elec}}\sum_j^{N_{elec}}\bigg(\langle\phi_i(q_i)\phi_j(q_j)|g_{ij}|\phi_i(q_i)\phi_j(q_j)\rangle-\langle\phi_i(q_i)\phi_j(q_j)|g_{ij}|\phi_j(q_i)\phi_i(q_j)\rangle\bigg)\\ =\frac{1}{2}\sum_i^{N_{elec}}\sum_j^{N_{elec}}\Bigg(\bigg\langle\delta\phi_i(q_i)\phi_j(q_j)\bigg|g_{ij}\bigg|\phi_i(q_i)\phi_j(q_j)\bigg\rangle+\bigg\langle\phi_i(q_i)\phi_j(q_j)\bigg|g_{ij}\bigg|\delta\phi_i(q_i)\phi_j(q_j)\bigg\rangle\\ -\bigg\langle\delta\phi_i(q_i)\phi_j(q_j)\bigg|g_{ij}\bigg|\phi_j(q_i)\phi_i(q_j)\bigg \rangle-\bigg\langle\phi_i(q_i)\phi_j(q_j)\bigg|g_{ij}\bigg|\phi_j(q_i)\delta\phi_i(q_j)\bigg \rangle \Bigg)\\ =\frac{1}{2}\sum_{i}^{N_{elec}}\sum_j^{N_{elec}}\bigg(\langle\delta\phi_i|J_j|\phi_i\rangle+\langle\phi_i|J_j|\delta\phi_i\rangle-\langle\delta\phi_i|K_j|\phi_i\rangle-\langle\phi_i|K_j|\delta\phi_i\rangle\bigg)\\ =\frac{1}{2}\sum_{i}^{N_{elec}}\sum_j^{N_{elec}}\bigg(\langle\delta\phi_i|J_j|\phi_i\rangle+\langle\delta\phi_i|J_j|\phi_i\rangle^*-\langle\delta\phi_i|K_j|\phi_i\rangle-\langle\delta\phi_i|K_j|\phi_i\rangle^*\bigg)\\ =\frac{1}{2}\sum_{i}^{N_{elec}}\sum_j^{N_{elec}}\bigg(\langle\delta\phi_i|J_j-K_j|\phi_i\rangle+\langle\delta\phi_i|J_j-K_j|\phi_i\rangle^*\bigg) \]

同理我们再对\(j\)轨道求,合并最终可得:

\[\delta E_2=\frac{1}{2}\sum_{i}^{N_{elec}}\sum_j^{N_{elec}}\bigg(\langle\delta\phi_i|J_j-K_j|\phi_i\rangle+\langle\delta\phi_i|J_j-K_j|\phi_i\rangle^*+\\ \langle\delta\phi_j|J_i-K_i|\phi_j\rangle+\langle\delta\phi_j|J_i-K_i|\phi_j\rangle^*\bigg) \]

由于上面一、三项和二、四项其实是在整个\(i,j\)求和的结果是相等的,因此可以合并成一项,去掉系数\(1/2\):

\[\delta E_2=\sum_{i}^{N_{elec}}\sum_j^{N_{elec}}\bigg(\langle\delta\phi_i|J_j-K_j|\phi_i\rangle+\langle\delta\phi_i|J_j-K_j|\phi_i\rangle^*\bigg) \]

(3)原子核-原子核排斥项:

\[\delta E_3=0 \]

综上,能量微分可表示为:

\[\delta E=\sum_i^{N_{elec}}\bigg(\langle\delta\phi_i|h_i|\phi_i\rangle+\langle\delta\phi_i|h_i|\phi_i\rangle^*\bigg)+\sum_{i}^{N_{elec}}\sum_j^{N_{elec}}\bigg(\langle\delta\phi_i|J_j-K_j|\phi_i\rangle+\langle\delta\phi_i|J_j-K_j|\phi_i\rangle^*\bigg) \]

我们令:

\[F_i=h_i+\sum_{j}^{N_{elec}}(J_j-K_j) \]

上述算符称为Fock算符,为单电子算符。注意其与能量算符的区别,Fock算符是通过变分原理得到的,能量算符并不等于Fock算符的求和。

进一步有:

\[\delta E=\sum_{i}^{N_{elec}}(\langle\delta\phi_i|F_i|\phi_i\rangle+\langle\delta\phi_i|F_i|\phi_i\rangle^*) \]

(4)正交归一化条件项:

\[=-\sum_{i}^{N_{elec}}\sum_{j}^{N_{elec}}\lambda_{ij}(\langle\delta\phi_i|\phi_j\rangle+\langle\phi_i|\delta\phi_j\rangle)\\ =-\sum_{i}^{N_{elec}}\sum_{j}^{N_{elec}}\lambda_{ij}(\langle\delta\phi_i|\phi_j\rangle+\langle\delta\phi_j|\phi_i\rangle^*)\\ \]

因为是对\(i,j\)全部求和,所以\(\langle\delta\phi_j|\phi_i\rangle^*\)项其实\(i、j\)符号可换,即\(\langle\delta\phi_i|\phi_j\rangle^{*}\),因此该项即为第一项的复共轭项

综上,泛函微分可变为:

\[\delta L=\sum_{i}^{N_{elec}}(\langle\delta\phi_i|F_i|\phi_i\rangle+\langle\delta\phi_i|F_i|\phi_i\rangle^*)-\sum_{i}^{N_{elec}}\sum_{j}^{N_{elec}}\lambda_{ij}(\langle\delta\phi_i|\phi_j\rangle+\langle\delta\phi_i|\phi_j\rangle^*)=0 \]

上述等式中实数项和复共轭项应该均为零,有:

\[\sum_{i}^{N_{elec}}\langle\delta\phi_i|F_i|\phi_i\rangle-\sum_{i}^{N_{elec}}\sum_{j}^{N_{elec}}\lambda_{ij}\langle\delta\phi_i|\phi_j\rangle=0\\ \sum_{i}^{N_{elec}}\langle\delta\phi_i|F_i|\phi_i\rangle^*-\sum_{i}^{N_{elec}}\sum_{j}^{N_{elec}}\lambda_{ij}\langle\delta\phi_i|\phi_j\rangle^*=0 \]

对上述复共轭等式取复共轭有:(简单理解就是去掉*号,不然展开写一下也行)

\[\sum_{i}^{N_{elec}}\langle\delta\phi_i|F_i|\phi_i\rangle-\sum_{i}^{N_{elec}}\sum_{j}^{N_{elec}}\lambda_{ji}^*\langle\delta\phi_i|\phi_j\rangle=0 \]

用上式与前面的实数项等式相减,可得:

\[\sum_{i}^{N_{elec}}\sum_{j}^{N_{elec}}(\lambda_{ij}-\lambda_{ji}^*)\langle\delta\phi_i|\phi_j\rangle=0\\ \Rightarrow \lambda_{ij}=\lambda_{ji}^* \]

所以所有的拉格朗日乘子构成了一个厄米矩阵

观察上述泛函微分式的实数项或复共轭项,提出左矢和一次加和:

\[\sum_{i}^{N_{elec}}\langle\delta\phi_i|\bigg[F_i|\phi_i\rangle-\sum_{j}^{N_{elec}}\lambda_{ij}|\phi_j\rangle \bigg]=0 \]

因为等式对所有的\(\delta\phi_i\)均满足,因此方括号中的式子恒为0,即得HF方程:

\[F_i\phi_i=\sum_{j}^{N_{elec}}\lambda_{ij}\phi_j \]

下面对拉格朗日乘子矩阵进行对角化,在此之间补充知识:

(1)酉变换(幺正变换)

表象:在量子力学中,任何一个量子态都可以看作是抽象的希尔伯特空间的一个”矢量“,而体系的任何一组力学量完全集\(F\)的共同本征态\(\{\psi_k\}\)构成此态空间的一组正交归一完备的基矢。量子态可以用不同的基矢来展开,如坐标表象、动量表象

例:

动量的本征函数为:\(\psi_p(x)=\frac{1}{\sqrt{2\pi\hbar}}exp(ipx/\hbar)\)

我们将波函数用上述动量本征函数展开:

\[\Psi(x,t)=\sum_{p}C(p,t)\psi_p(x) \]

由于动量本征函数是连续的,改为积分形式:

\[\Psi(x,t)=\int_p C(p,t)\psi_p(x)dp \]

积分中第一项就是展开系数,第二项是展开的基底函数,整个积分就等价于一个向量(希尔伯特空间内)。进一步我们可以认为展开系数就是该向量在基底上的投影

那么由内积的定义,我们可以得到系数表达式即为:

\[C(p,t)=\int \psi_p^*(x)\Psi(x,t)dx=\langle\psi_p(x)|\Psi(x,t)\rangle \]

下面考虑将波函数用按照某一完备基展开:

\[|\psi\rangle=\sum_na_n|Q_n\rangle \]

则:

\[\langle Q_m|\psi\rangle = \sum_{n}a_n\langle Q_m|Q_n\rangle=\sum_n a_n\delta_{mn}=a_m \]

代回波函数可得:

\[|\psi\rangle=\sum_n |Q_n\rangle\langle Q_n|\psi\rangle \]

观察我们可知内积\(\langle Q_n|\psi\rangle\)的结果是一个常量,表示的波函数在\(Q_n\)为基底的空间的投影大小,则\(|Q_n\rangle\langle Q_n|\psi\rangle\)

则是波函数在\(Q_n\)为基矢的空间的波函数。同时因为有\(|Q_n\rangle\langle Q_n|\psi\rangle=\langle Q_n|\psi\rangle |Q_n\rangle\),可以理解为\(|Q_n\rangle\langle Q_n|\)是一个算符,将波函数\(\psi\)映射到了\(Q_n\)上

可知有:

\[\sum_n|Q_n\rangle\langle Q_n|=I \]

上述方程称为完备性关系或封闭性

考虑两套归一化的基矢,用\(B\)表象的本征函数展开\(A\)表象中的左矢\(\psi_A\):(如果基底函数是连续函数可写为积分)

\[|\psi_A\rangle=\sum_m|\phi_m\rangle\langle\phi_m|\psi_A\rangle \quad(\phi为B基矢的完备基) \]

令\(\langle\phi_m|\psi_A\rangle=S_{mA}\),则有:

\[|\psi_A\rangle=\sum_m|\phi_m\rangle S_{mA} \]

其中\(S_{mA}\)就是将\(\psi_A\)在\(B\)表象下展开得到的系数矩阵,\(S\)就是基矢变换的变换矩阵

上述等式用矩阵形式可以表达为:

\[\begin{pmatrix}\psi_1 \\\psi_2\\\vdots\\\psi_A\end{pmatrix}=\begin{pmatrix}S_{11} & S_{21} & \cdots & S_{m1} \\ S_{12} & S_{22} & \cdots & S_{m2}\\ \vdots & \vdots & \ddots & \vdots \\ S_{1A} & S_{2A} & \cdots & S_{mA}\end{pmatrix}_{A\times m}\begin{pmatrix}\phi_1 \\ \phi_2 \\ \vdots \\ \phi_m\end{pmatrix} \]

展开右矢有:(注意左矢和右矢互为转置共轭)

\[\langle\psi_A|=\sum_n\langle\psi_A|\phi_n\rangle\langle\phi_n| \]

令\(\langle\psi_A|\phi_n\rangle=S_{nA}^{\dagger}\),则有:

\[\langle\psi_A|=\sum_n S_{nA}^{\dagger}\langle\phi_n| \]

矩阵形式表达为:

\[\begin{pmatrix} \psi_1^* & \psi_2^* & \cdots & \psi_A^*\end{pmatrix}=\begin{pmatrix} \phi_1^* & \phi_2^* & \cdots & \phi_n^*\end{pmatrix}\begin{pmatrix} S_{11}^* & S_{12}^* & \cdots & S_{1A}^*\\S_{21}^* & S_{22}^* &\cdots & S_{2A}^* \\ \vdots & \vdots & \ddots & \vdots \\ S_{n1}^* & S_{n2}^* & \cdots & S_{nA}^* \end{pmatrix}_{n\times A} \]

则有:

\[\langle\psi_A|\psi_A\rangle=\sum_m\sum_n S_{nA}^{\dagger}\langle\phi_n|\phi_m\rangle S_{mA}\\ \Rightarrow 1=\sum_{mn}S_{nA}^{\dagger}\delta_{nm} S_{mA}\\ \Rightarrow \sum_n S_{nA}^{\dagger}S_{nA}=(S^{\dagger}S)_{AA}=I\\ (注意S_{nA}矩阵在左边,S_{nA}^{\dagger}在右边,最终矩阵大小为A\times A) \]

我们可以进一步计算得到:\(SS^{\dagger}=I\)

参考下面链接:

量子力学变换幺正变换 - 百度文库 (baidu.com)

因此我们的结论是两个表象之间的变化矩阵满足:\(S^{\dagger}=S^{-1}\)

满足上式的矩阵称为幺正矩阵,由幺正矩阵所表示的变化称为幺正变换,即由一个表象到另一个表象的变换为幺正变换

(2)矩阵对角化复习

对于同阶方阵\(A、B\),如果存在可逆矩阵\(U\),使得\(B=U^{-1}AU\),则称\(A\)与\(B\)是相似的,记为\(A\sim B\)。有如下性质: 相似矩阵的行列式的值相等 相似矩阵或者都可逆或者都不可逆,在可逆的情况下,逆矩阵也相似 若\(A\sim B\),则\(A^n\sim B^n\)

若\(A\)和对角阵\(B\)相似,如何求可逆矩阵?

\[\because B=U^{-1}AU\\ \therefore AU=UB\\ 假设B=\begin{pmatrix}\lambda_1 & & & 0 \\ & \lambda_2 & & \\ & & \ddots \\0 & & & \lambda_n\end{pmatrix}\\ 记可逆矩阵U的列向量为: U=(\alpha_1,\alpha_2,\cdots,\alpha_n)\\ \therefore A(\alpha_1,\alpha_2,\cdots,\alpha_n)=(\alpha_1,\alpha_2,\cdots,\alpha_n)\begin{pmatrix}\lambda_1 & & & 0 \\ & \lambda_2 & & \\ & & \ddots \\0 & & & \lambda_n\end{pmatrix}\\ \therefore (A\alpha_1,A\alpha_2,\cdots,A\alpha_n)=(\lambda_1\alpha_1,\lambda_2\alpha_2,\cdots,\lambda_n\alpha_n) \]

所以必须满足\(A\alpha_1=\lambda_1\alpha_1,A\alpha_2=\lambda_2\alpha_2,\cdots,A\alpha_n=\lambda_n\alpha_n\)

特征值和特征向量

设\(A\)是\(n\)阶方阵,如果数\(\lambda\)和非零的列向量满足\(A\alpha=\lambda\alpha\),则\(\lambda\)为\(A\)的特征值,\(\alpha\)为\(A\)属于\(\lambda\)的特征向量

几何意义:对于一个向量\(x\),我们将它乘上一个矩阵\(A\),相当于进行一次线性变换,变换后的向量\(Ax\)的方向和长度都发生了变化。而对于一个特定的矩阵,总存在一些特定方向的向量,使得变换后的向量只是改变长度而不改变方向,则该向量即称为矩阵的特征向量,对应的长度即特征向量的特征值

由上有齐次线性方程组:

\[(\lambda I-A)\alpha=0 \]

而要有非零解,行列式的值必须为0,否则满秩,只有唯一解(零解)

则得到特征方程:(\(\lambda\)的\(n\)次方程)

\[|\lambda I-A|=0 \]

展开特征多项式:

\[\left|\begin{array}{c} \lambda_1-a_{11} & -a_{12} & \cdots & -a_{1n} \\ -a_{21} & \lambda_2-a_{22} & \cdots & -a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ -a_{n1} & \cdots & \cdots & \lambda_n-a_{nn} \end{array}\right| \]

只有对角线有\(\lambda\)未知数,所以可知特征值就是方程组的解

对于\(n\)阶方阵,有如下性质:

\[\lambda_1+\lambda_2+\cdots+\lambda_n=a_{11}+a_{22}+\cdots+a_{nn}\\ \lambda_1\lambda_2\cdots\lambda_n=|A| \]

特征值和特征向量的解法:

第一步:计算行列式\(|\lambda I-A|\),求出特征方程的全部根,即得到\(A\)的全部特征值

第二步:对每个特征值\(\lambda_0\),求齐次线性方程组\((\lambda_0I-A)x=0\)的一个基础解系\(\alpha_1,\alpha_2,\cdots\alpha_s\),并写成列向量的形式,则\(A\)属于\(\lambda_0\)的全部特征向量为:\(k_1\alpha_1+k_2\alpha_2+\cdots+k_s\alpha_s\)

例:

\[A=\begin{pmatrix}1 & 2 & 2 \\ 2 & 1 & 2 \\ 2 & 2 & 1\end{pmatrix},\text{求A的全部特征值和特征向量}\\ |\lambda I-A|=\left| \begin{array}{c}\lambda-1 & -2 & -2 \\ -2 & \lambda-1 & -2 \\ -2 & -2 & \lambda-1\end{array}\right|=(\lambda+1)^2(\lambda-5)=0\\ \Rightarrow \lambda_1=-1,\lambda_2=5\\ 对于\lambda_1,代入方程组有:(有两个独立变量,线性无关)\\ \begin{cases}-2x_1-2x_2-2x_3=0 \\-2x_1-2x_2-2x_3=0 \\ -2x_1-2x_2-2x_3=0\end{cases}\\ 取x_2=1,x_3=0时,x_1=-1\\ 取x_2=0,x_3=1时,x_1=-1\\ \therefore \lambda_1的特征向量为:k_1\begin{pmatrix}-1\\1\\0\end{pmatrix}+k_2\begin{pmatrix}-1\\0\\1\end{pmatrix}\\ 对于\lambda_2=5,代入方程组有:\\ \begin{cases}4x_1-2x_2-2x_3=0 \\-2x_1+4x_2-2x_3=0 \\ -2x_1-2x_2+4x_3=0\end{cases}\Rightarrow\begin{cases}2x_1-x_2-x_3=0\\0x_1+ x_2-x_3=0\end{cases}\\ 取x_2=1时,则x_3=1,x_1=1\\ \therefore \lambda_2的特征向量为:k\begin{pmatrix}1\\1\\1\end{pmatrix},k\neq0 \]

补充性质:相似矩阵具有相同的特征多项式和特征值

矩阵对角化

定理1:\(n\)阶方阵\(A\)可对角化的充要条件是\(A\)有\(n\)个线性无关的特征向量

定理2:方阵\(A\)属于不同特征值的特征向量线性无关(故只要分别考虑每个特征值对应是否满足即可)

定理3:若\(n\)阶方阵\(A\)有\(n\)个互异的特征值,则\(A\)可对角化

对角化的流程:

1)求解\(|\lambda I-A|=0\)的解,得到全部特征值

2)代入特征值确定基础解系,将不同特征值的基础解系排列得到向量,若个数等于矩阵的阶数,则可对角化

3)最后可逆矩阵\(U=(\alpha_1,\cdots,\alpha_n)\),即所有特征向量(列向量)的组合

4)对角化得到的矩阵\(\Lambda=U^{-1}AU=\begin{pmatrix}\lambda_1 & & & 0 \\ & \lambda_2 & & \\ & & \ddots \\0 & & & \lambda_n\end{pmatrix}\)

或者看下面的:(\(n\)阶矩阵推广即可)

假设\(P=(v_1,v_2)\),其中\(v_1,v_2\)分别为矩阵\(A\)的特征向量(列向量)

则\(AP=A(v_1,v_2)=A\begin{pmatrix}v_{11} & v_{21} \\ v_{12} & v_{22}\end{pmatrix}=\begin{pmatrix}A\begin{bmatrix}v_{11}\\v_{12}\end{bmatrix} &A\begin{bmatrix}v_{21}\\v_{22}\end{bmatrix}\end{pmatrix}\)

因为有:\(Av_1=\lambda_1v_1,Av_2=\lambda_2v_2\)

则\(AP=\begin{pmatrix}\lambda_1v_{11} & \lambda_2v_{21} \\ \lambda_1v_{12} & \lambda_2v_{22}\end{pmatrix}=\begin{pmatrix}v_{11} & v_{21} \\ v_{12} & v_{22}\end{pmatrix}\begin{pmatrix}\lambda_1 & 0\\ 0 & \lambda_2\end{pmatrix}=P\Lambda\)

则\(P^{-1}AP=\Lambda\)

回顾HF方程为:

\[F_i\phi_i=\sum_{j}^{N_{elec}}\lambda_{ij}\phi_j \]

展开为矩阵:

\[\begin{pmatrix}F_1&&&\\&F_2&&\\&&\ddots&\\&&&F_N \end{pmatrix}\begin{pmatrix}\phi_1\\\phi_2\\\vdots\\\phi_N\end{pmatrix}=\begin{pmatrix}\lambda_{11}&\lambda_{12}&\cdots&\lambda_{1N}\\\lambda_{21}&\lambda_{22}&\cdots&\lambda_{2N}\\\vdots&\vdots&\ddots&\vdots\\\lambda_{N1}&\lambda_{N2}&\cdots&\lambda_{NN}\end{pmatrix}\begin{pmatrix}\phi_1\\\phi_2\\\vdots\\\phi_N\end{pmatrix} \]

即\(F\Phi=\lambda\Phi\)

取酉变换矩阵\(U\),则\(UU^{-1}=U^{-1}U=I\),有:

\[F\Phi=\lambda(U^{-1}U)\Phi\\ FU\Phi=U\lambda(U^{-1}U)\Phi\\ F(U\Phi)=(U\lambda U^{-1})(U\Phi) \]

因此可以确定一个酉矩阵\(U\),使得\(\lambda\)对角化:

\[U\lambda U^{-1}=\begin{pmatrix}\varepsilon_1&&&\\&\varepsilon_2&&\\&&\ddots&\\&&&\varepsilon_N \end{pmatrix} \]

此时原来的分子轨道\(\Phi\)经过酉变换得到了一个新的轨道\(\Phi^{'}=U\Phi\),称之为正则分子轨道

因此HF方程即变为正则HF方程:

\[F_i\phi_i^{'}=\varepsilon_i\phi_i^{'} \]

将上述方程积分:

\[\langle\phi_i^{'}|F_i|\phi_i^{'}\rangle=\varepsilon_i\langle\phi_i^{'}|\phi_i^{'}\rangle=\varepsilon_i \]

可以看到对角化后的拉格朗日乘子是Fock算符在正则轨道下的期望值,其物理意义可以看作是分子轨道的能量 由于HF算符的具体形式需要由所有的分子轨道确定,因而要求解该方程必须采用自洽迭代(SCF)的方法,得到解就是Self-Consistent Field orbitals

酉变换不会改变总波函数

当我们确定了正则分子轨道,就可以通过线性组合的方法得到其它的分子轨道

最后,我们计算一下看看总能量与正则HF方程确定的轨道能量的差别:

前面确定的近似波函数的总能量表达式为:

\[E=\sum_{i=1}^{N_{elec}}\bar{h}_i+\frac{1}{2}\sum_{i=1}^{N_{elec}}\sum_{j=1}^{N_{elec}}(J_{ij}-K_{ij})+V_{nn}\\ =\sum_i^{N_{elec}}\int\phi_i^*\bigg[h_i+\frac{1}{2}\sum_j^{N_{elec}}(J_j-K_j)\bigg]\phi_idx+V_{nn}\\ =\frac{1}{2}\sum_i^{N_{elec}}\int\phi_i^*\bigg[h_i+(h_i+\sum_j^{N_{elec}}(J_j-K_j))\bigg]\phi_idx+V_{nn}\\ =\frac{1}{2}\sum_i^{N_{elec}}\int\phi_i^*\bigg[h_i+F_i\bigg]\phi_idx+V_{nn}\\ =\frac{1}{2}\sum_i^{N_{elec}}(\bar{h}_i+\varepsilon_i)+V_{nn}\\ =\sum_i^{N_{elec}}(\frac{1}{2}\bar{h}_i+\frac{1}{2}\varepsilon_i)+V_{nn}\\ =\sum_i^{N_{elec}}(\varepsilon_i+\frac{1}{2}(\bar{h}_i-\varepsilon_i))+V_{nn}\\ \because \varepsilon_i=\langle\phi_i|F_i|\phi_i\rangle=\langle\phi_i|h_i+\sum_j^{N_{elec}}(J_j-K_j)|\phi_i\rangle=\bar{h}_i+\sum_j^{N_{elec}}(J_{ij}-K_{ij})\\ \therefore \bar{h}_i-\varepsilon_i=-\sum_j^{N_{elec}}(J_{ij}-K_{ij})\\ \quad= \sum_i^{N_{elec}}\bigg[\varepsilon_i-\frac{1}{2}\sum_j^{N_{elec}}(J_{ij}-K_{ij})\bigg]+V_{nn}\neq\sum_i^{N_{elec}}\varepsilon_i+V_{nn} \]

可以看到,总能量并不是Fock方程的轨道能量的加和,而是要去掉一个双电子项,这是因为定义Fock算符时,双电子排斥项重复计算了,故要去掉双电子排斥多的部分能量

3.5 Koopmans' Theorem

考虑第一解离能和电子亲合能(这里只考虑解离能,亲合能同理),假设N电子体系和N-1电子体系的分子轨道相同,有:

\[E_N=\sum_{i=1}^{N_{elec}}\bar{h}_i+\frac{1}{2}\sum_{i=1}^{N_{elec}}\sum_{j=1}^{N_{elec}}(J_{ij}-K_{ij})+V_{nn}\\ E_{N-1}=\sum_{i=1}^{N_{elec}-1}\bar{h}_i+\frac{1}{2}\sum_{i=1}^{N_{elec}-1}\sum_{j=1}^{N_{elec}-1}(J_{ij}-K_{ij})+V_{nn} \]

则:

\[E_N-E_{N-1}=\bar{h}_k+\frac{1}{2}\sum_{i=1}^{N_{eleci}}(J_{ik}-K_{ik})+\frac{1}{2}\sum_{j=1}^{N_{eleci}}(J_{kj}-K_{kj})\\ =\bar{h}_k+\sum_{i=1}^{N_{eleci}}(J_{ik}-K_{ik})=\varepsilon_k \]

同理:

\[E_{N+1}-E_N=\varepsilon_k \]

因此Koopmans提出假设在分子轨道固定(Frozen MO approximation)的前提下(即电子的得失不会改变体系的分子轨道,但实际上有电子占据和无电子占据的轨道的本征值是不一样的),第一解离能和电子亲和能的值为对应的分子轨道的能量。

3.6 Roothaan-Hall 方程

在一般的计算中,都是采用一组已知的函数作为基底用来展开分子轨道,这些函数可以是指数函数、高斯函数、多项式、立方函数、平面波等等,选取的原则是保证易收敛且积分计算方便。

采用如\(e^{-\zeta r}\)的指数函数作为基底函数,虽然其在描述电子结构上更为准确(实际上该函数为H原子的确解),但是在处理双电子积分项时计算复杂度为\(O(N^4)\);而采用高斯函数可以使得复杂度为\(O(N^4)\)的问题变为复杂度为\(O(N^2)\)的问题,同时其在较远处与指数函数较为接近,而在中心处则偏离比较大,可以通过多个高斯函数的线性组合逼近,后面详述。

image-20211018171421738 image-20211111224353755

假设基底函数已经确定,则可将分子轨道展开为基底函数的线性组合(MO=LCAO, Linear Combination of Atomic Orbitals):

\[\phi_i=\sum_{\alpha}^{M_{basis}}c_{\alpha i}\chi_\alpha \]

代入正则HF方程:(我们假设共有n个轨道,每个轨道由m个基底函数展开)

\[F_i\sum_{\alpha}^{M_{basis}}c_{\alpha i}\chi_\alpha=\varepsilon_i\sum_{\alpha}^{M_{basis}}c_{\alpha i}\chi_\alpha \quad(i=1,2,\cdots n) \]

在方程两边分别左乘一个基底函数并积分,注意加和\(\sum\)和系数可以提到前面:

\[\sum_{\alpha}^{M_{basis}}c_{\alpha i}\int \chi_\beta^*F_i\chi_\alpha dx=\varepsilon_i\sum_{\alpha}^{M_{basis}}c_{\alpha i}\int\chi_{\beta}^*\chi_\alpha dx \]

令:

\[F_{\beta\alpha(i)}=\int \chi_\beta^*F_i\chi_\alpha dx=\langle\chi_{\beta}|F_i|\chi_\alpha\rangle\\ S_{\beta\alpha(i)}=\int\chi_{\beta}^*\chi_\alpha dx=\langle\chi_\beta|\chi_{\alpha}\rangle \]

注意上面\(F\)和\(S\)定义其实分别定义了n个积分

则上述方程可变为:(即Roothaan-Hall方程)

\[\sum_{\alpha}^{M_{basis}}c_{\alpha i}F_{\beta\alpha(i)}=\varepsilon_i\sum_{\alpha}^{M_{basis}}c_{\alpha i}S_{\beta\alpha(i)} \]

我们分别将方程左右两边展开看:

方程左边:

\[\begin{pmatrix} \sum_{\alpha}^{M_{basis}}c_{\alpha i}F_{1\alpha(i)}\\ \sum_{\alpha}^{M_{basis}}c_{\alpha i}F_{2\alpha(i)}\\ \vdots\\ \sum_{\alpha}^{M_{basis}}c_{\alpha i}F_{m\alpha(i)} \end{pmatrix}= \begin{pmatrix} \sum_{\alpha}^{M_{basis}}c_{\alpha 1}F_{1\alpha(1)}&\sum_{\alpha}^{M_{basis}}c_{\alpha 2}F_{1\alpha(2)}&\cdots&\sum_{\alpha}^{M_{basis}}c_{\alpha n}F_{1\alpha(n)}\\ \sum_{\alpha}^{M_{basis}}c_{\alpha 1}F_{2\alpha(1)}&\sum_{\alpha}^{M_{basis}}c_{\alpha 2}F_{2\alpha(2)}&\cdots&\sum_{\alpha}^{M_{basis}}c_{\alpha n}F_{2\alpha(n)}\\ \vdots&\vdots&\ddots&\vdots\\ \sum_{\alpha}^{M_{basis}}c_{\alpha 1}F_{m\alpha(1)}&\sum_{\alpha}^{M_{basis}}c_{\alpha 2}F_{m\alpha(2)}&\cdots&\sum_{\alpha}^{M_{basis}}c_{\alpha n}F_{m\alpha(n)} \end{pmatrix} \]

其中从左向右是是分别取每个分子轨道的Fock算符,从上到下是分别取每个基底函数

接着我们再看右边的每一项:

\[=\begin{pmatrix} [c_{11}F_{11(1)}+c_{21}F_{12(1)}+\cdots] & [c_{12}F_{11(2)}+c_{22}F_{12(2)}+\cdots] & \cdots &[c_{1n}F_{11(n)}+c_{2n}F_{12(n)}+\cdots]\\ [c_{11}F_{21(1)}+c_{21}F_{22(1)}+\cdots] & [c_{12}F_{21(2)}+c_{22}F_{22(2)}+\cdots] & \cdots &[c_{1n}F_{21(n)}+c_{2n}F_{22(n)}+\cdots]\\ \vdots & \vdots & \ddots & \vdots\\ [c_{11}F_{m1(1)}+c_{21}F_{m2(1)}+\cdots] & [c_{12}F_{m1(2)}+c_{22}F_{m2(2)}+\cdots] & \cdots &[c_{1n}F_{m1(n)}+c_{2n}F_{m2(n)}+\cdots] \end{pmatrix} \]

\[=\begin{pmatrix} \mathbf{F_{11}} & \mathbf{F_{12}} & \cdots & \mathbf{F_{1m}}\\ \mathbf{F_{21}} & \mathbf{F_{22}} & \cdots & \mathbf{F_{2m}}\\ \vdots & \vdots & \ddots & \vdots\\ \mathbf{F_{m1}} & \mathbf{F_{m2}} & \cdots & \mathbf{F_{mm}} \end{pmatrix} \begin{pmatrix} c_{11}&c_{12}&\cdots&c_{1n}\\ c_{21}&c_{22}&\cdots&c_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ c_{m1}&c_{m2}&\cdots&c_{mn} \end{pmatrix} \]

令:

\[\mathbf{F}=\begin{pmatrix} \mathbf{F_{11}} & \mathbf{F_{12}} & \cdots & \mathbf{F_{1m}}\\ \mathbf{F_{21}} & \mathbf{F_{22}} & \cdots & \mathbf{F_{2m}}\\ \vdots & \vdots & \ddots & \vdots\\ \mathbf{F_{m1}} & \mathbf{F_{m2}} & \cdots & \mathbf{F_{mm}} \end{pmatrix}\\ \mathbf{C}= \begin{pmatrix} c_{11}&c_{12}&\cdots&c_{1n}\\ c_{21}&c_{22}&\cdots&c_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ c_{m1}&c_{m2}&\cdots&c_{mn} \end{pmatrix} \]

其中要注意\(F\)的每一项需要根据\(C\)的每一项的下标的第二位改变为第\(i\)个分子轨道的Fock算符与基底函数的积分

方程右边:

\[\begin{pmatrix} \varepsilon_i\sum_{\alpha}^{M_{basis}}c_{\alpha i}S_{1\alpha(i)}\\ \varepsilon_i\sum_{\alpha}^{M_{basis}}c_{\alpha i}S_{2\alpha(i)}\\ \vdots\\ \varepsilon_i\sum_{\alpha}^{M_{basis}}c_{\alpha i}S_{m\alpha(i)} \end{pmatrix} = \begin{pmatrix} \varepsilon_1\sum_{\alpha}^{M_{basis}}c_{\alpha 1}S_{1\alpha(1)}&\varepsilon_2\sum_{\alpha}^{M_{basis}}c_{\alpha 2}S_{1\alpha(2)}&\cdots&\varepsilon_n\sum_{\alpha}^{M_{basis}}c_{\alpha n}S_{1\alpha(n)}\\ \varepsilon_1\sum_{\alpha}^{M_{basis}}c_{\alpha 1}S_{2\alpha(1)}&\varepsilon_2\sum_{\alpha}^{M_{basis}}c_{\alpha 2}S_{2\alpha(2)}&\cdots&\varepsilon_n\sum_{\alpha}^{M_{basis}}c_{\alpha n}S_{2\alpha(n)}\\ \vdots&\vdots&\ddots&\vdots\\ \varepsilon_1\sum_{\alpha}^{M_{basis}}c_{\alpha 1}S_{m\alpha(1)}&\varepsilon_2\sum_{\alpha}^{M_{basis}}c_{\alpha 2}S_{m\alpha(2)}&\cdots&\varepsilon_n\sum_{\alpha}^{M_{basis}}c_{\alpha n}S_{m\alpha(n)} \end{pmatrix} \]

\[=\begin{pmatrix} \varepsilon_1(c_{11}S_{11}+c_{21}S_{12}+\cdots) & \varepsilon_2(c_{12}S_{11}+c_{22}S_{12}+\cdots) & \cdots & \varepsilon_n(c_{1n}S_{11}+c_{2n}S_{12}+\cdots) \\ \varepsilon_1(c_{11}S_{21}+c_{21}S_{22}+\cdots) & \varepsilon_2(c_{12}S_{21}+c_{22}S_{22}+\cdots) & \cdots & \varepsilon_n(c_{1n}S_{21}+c_{2n}S_{22}+\cdots) \\ \vdots & \vdots & \ddots & \vdots \\ \varepsilon_1(c_{11}S_{m1}+c_{21}S_{m2}+\cdots) & \varepsilon_2(c_{12}S_{m1}+c_{22}S_{m2}+\cdots) & \cdots & \varepsilon_n(c_{1n}S_{m1}+c_{2n}S_{m2}+\cdots) \end{pmatrix} \]

\[=\begin{pmatrix} S_{11} & S_{12} & \cdots & S_{1m} \\ S_{21} & S_{22} & \cdots & S_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ S_{m1} & S_{m2} & \cdots & S_{mm} \end{pmatrix} \begin{pmatrix} c_{11} & c_{12} & \cdots & c_{1n} \\ c_{21} & c_{22} & \cdots & c_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ c_{m1} & c_{m2} & \cdots & c_{mn} \end{pmatrix} \begin{pmatrix} \varepsilon_{11} & 0 & \cdots & 0 \\ 0 & \varepsilon_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \varepsilon_{nn} \end{pmatrix}\\ =\mathbf{SC}\varepsilon \]

综上我们可以得到Roothaan-Fall方程的矩阵形式:

\[\mathbf{FC}=\mathbf{SC}\varepsilon \]

我们先回顾下目前为止一共进行了哪些推导:

(1)总波函数可以表示为不同量子态(本征态)的波函数的线性组合,由于引入了BO近似,限制波函数处于一个本征态,即\(\Psi_k\)

(2)由绝热近似,可以得到某个本征态下原子核的薛定谔方程,而方程中需要电子的能量,即要求解电子的薛定谔方程

(3)全电子波函数可以用Slater行列式表示,其由一系列单电子的自旋轨道乘积组成

(4)将全电子波函数代入电子哈密顿算符作用在波函数上的全空间积分,可得期望值即能量

(5)利用变分法可得HF方程,对角化后得到正则HF方程

(6)用基底函数展开单电子波函数,即得Roothaan-Fall方程

下面我们对Roothaan-Fall方程的各项的具体形式逐个分析:

Fock矩阵:

\[F_{\beta\alpha(i)}=\int\chi_\beta^*F_i\chi_\alpha dx \]

注意也有书中用\(\hat{F}\)代替\(F_i\)统一表示第\(i=1,2,\cdots n\)电子的Fock算符,这里为了与前面统一,采用单电子形式

Fock算符为:

\[F_i=h_i+\sum_{j=1}^{N_{elec}}(J_{j}-K_j) \]

则:

\[F_{\beta\alpha(i)}=\int\chi_{\beta}^*\bigg[h_i+\sum_{j=1}^{N_{elec}}(J_j-K_j)\bigg]\chi_\alpha dx \]

其中:

\[h_i=-\frac{1}{2}\nabla_i^2-\sum_{A=1}^{N_{nuc}}\frac{Z_A}{|R_A-r_i|}\\ J_j=\int \phi_j^*(q_j)\frac{1}{|r_i-r_j|}\phi_j(q_j) dx\\ K_j=\int\phi_j^*(q_j)\frac{1}{|r_i-r_j|}\phi_j(q_i)dx \]

我们观察Roothaan-Fall方程的矩阵形式:

\[\mathbf{FC}=\mathbf{SC}\varepsilon \]

上述形式并不是标准的本征方程形式,因此我们需要对其进行变化

这里需要利用Lowdin正交化方法:包括对称正交化(Symmetric Orthogonalization)以及正则正交化(Canonical Orthogonalizaiton):

矩阵\(S\)的形式为:\(S=\langle\chi_k|\chi_l\rangle\)

(1)重叠积分矩阵\(S\)对角化:

\[S=PDP^{-1} \]

要求特征矩阵P,只需要求\(S\)的特征向量,即求解\(|AI-S|=0\),详细可参考之前的矩阵对角化补充知识

最终我们就可以得到对角化后的矩阵\(D=diag(\lambda_1,\cdots,\lambda_{M_{basis}})\)

(2)矩阵\(S\)为半正定矩阵:

假设\(S\)的归一化的特征向量为\(c_i\),有:(注意两个下标时,第一个表示行,第二个表示列)

\[Sc_i = \lambda_ic_i\\ c_i^{\dagger}c_i = \sum_j c_{ji}^{*}c_{ji}=\sum_j|c_{ji}|^2=1 \]

构造下面的函数:

\[\psi_i = \sum_kc_{ki}\chi_k \]

则:

\[\langle \psi_i|\psi_i\rangle =\langle \sum_jc_{ji}\chi_j|\sum_kc_{ki}\chi_k\rangle = \sum_j\sum_k c_{ji}^* S_{jk}c_{ki}\\= \sum_jc_{ji}^*(\sum_k S_{jk}c_{ki})=\sum_jc_{ji}^*\lambda_ic_{ji}=\lambda_i\sum_jc_{ji}^*c_{ji}\\ =\lambda_i\sum_j|c_{ji}|^2=\lambda_i \]

由于\(\langle\psi_i|\psi_i\rangle \geq 0\),所以\(\lambda_i\geq0\)

若矩阵\(S\)的所有特征向量\(c_i\)线性无关,当且仅当所有特征向量为零时,积分\(\langle\psi_i|\psi_i\rangle = 0\)才能成立,因此积分必定大于零;若矩阵\(S\)中的特征向量存在线性相关,则行列式\(|S|=0\),而\(|S|=\lambda_1\lambda_2\cdots\lambda_n\),因此必定存在特征值为零,故积分\(\langle\psi_i|\psi_i\rangle = 0\)能成立 因此有当矩阵\(S\)的特征向量线性无关时,矩阵\(S\)为正定阵;当矩阵\(S\)中存在特征向量线性相关时,矩阵\(S\)为半正定阵;因此我们知道矩阵\(S\)至少为半正定阵 我们可以通过矩阵\(S\)的本征值的最小值来衡量矩阵\(S\)对应的基底函数的相关性:

我们构造:

\[\psi=\sum_j^{m}d_j\chi_j \]

其中\(\chi_j\)即组成矩阵\(S\)的基底函数,则\(d\)为任意归一化(没有正交化)的展开系数(\(\sum_j|d_j|^2=1\))

由于重叠矩阵是两个函数的内积,故其为厄米矩阵。虽然厄米矩阵的属于不同本征值的特征向量相互正交,但是属于相同本征值的特征函数不一定正交。而根据谱定理,当我们已经确定了一个不变子空间invariant subspace,也就是这个空间里是由所有满足:\(Av=\lambda v\) 特征方程的相互正交的特征向量构成的,即\(V=\{v_1,v_2,\cdots,v_k\}\),而与该空间正交的空间也是厄米矩阵\(A\)的不变子空间。定理指出,在每个不变子空间中,都至少有一个特征向量。因此,我们要得到全部的特征向量,只要在已有不变子空间的基础上,求其正交空间,接着再求这个新的正交空间的正交空间,一直求就能确定全部正交归一化的特征向量。这个好像就是斯密特正交化的原理。因此最后全部的正交归一化后的特征向量可构成了一个完备集。

谱定理参考内容:

https://en.wikipedia.org/wiki/Spectral_theorem

https://www.youtube.com/watch?v=WAhM77lBkPk

厄米矩阵内容参考知识3:

量子化学知识点3-厄米矩阵 - miccoui - 博客园 (cnblogs.com)

因此上面的展开系数可表示为:(\(c\)为特征向量)

\[d=\sum_i\eta_ic_i\\ \sum_i|\eta_i|^2=1 \]

则:

\[\langle\psi|\psi\rangle = \langle\sum_jd_j\chi_j|\sum_kd_k\chi_k\rangle=\sum_{j,k}d_j^*S_{jk}d_k\\ =d^{\dagger}Sd=\sum_{i,j}\eta_i^*c_i^{\dagger}Sc_j\eta_j=\sum_{i,j}\eta_i^*\eta_j\lambda_jc_i^{\dagger}c_j\\=\sum_{i,j}\eta_i^*\eta_j\lambda_j\delta_{ij}=\sum_i|\eta_i|^2\lambda_i \]

设特征值中最大和最小的值分别为:\(\lambda_{max},\lambda_{min}\),因此有:

\[\langle\psi|\psi\rangle = \sum_i|\eta_i|^2\lambda_i \geq \sum_i|\eta_i|^2\lambda_{min}=\lambda_{min} \sum_i|\eta_i|^2=\lambda_{min}\\ \Rightarrow \langle\psi|\psi\rangle \geq \lambda_{min}\\ \langle\psi|\psi\rangle = \sum_i|\eta_i|^2\lambda_i \leq \sum_i|\eta_i|^2\lambda_{max}=\lambda_{max} \sum_i|\eta_i|^2=\lambda_{max}\\ \Rightarrow \langle\psi|\psi\rangle \leq \lambda_{max}\\ \Rightarrow \lambda_{min} \leq \langle\psi|\psi\rangle \leq \lambda_{max} \]

当组成矩阵\(S\)的基底函数线性相关或者接近线性相关,\(\lambda_{min}\)为零或接近零,这样会导致我们无法求矩阵的逆,因此我们可以用\(\lambda_{min}\)来判断,当其值过小时,则不能采用对称正交化方法,而要改用正则正交化方法

正定阵知识,之后有时间复习:

Definite matrix - Wikipedia

(3)求矩阵\(S^{-1/2}\)

Square root of a matrix - Wikipedia

image-20211110235140006

根据维基百科的解释,当矩阵\(A\)是正定阵或半正定阵时,若其对应的特征值为\(\lambda_1,\lambda_2,\cdots,\lambda_n\),均大于等于零,\(A\)对角化为:

\[A=Pdiag(\lambda_1,\lambda_2,\cdots,\lambda_n)P^{-1} \]

取矩阵\(B=Pdiag(\sqrt{\lambda_1},\sqrt{\lambda_2},\cdots,\sqrt{\lambda_n})P^{-1}\),则有:

\[BB=B^{T}B=\\ Pdiag(\sqrt{\lambda_1},\sqrt{\lambda_2},\cdots,\sqrt{\lambda_n})P^{-1}Pdiag(\sqrt{\lambda_1},\sqrt{\lambda_2},\cdots,\sqrt{\lambda_n})P^{-1}\\ =Pdiag(\lambda_1,\lambda_2,\cdots,\lambda_n)P^{-1}=A \]

如果有\(B^2=A\)成立,我们称\(B=A^{1/2}\),即矩阵的square root 的definition是建立在square的definition下的

重叠矩阵\(S\)的\(-1/2\)次幂即为:

\[S^{-\frac{1}{2}}=PD^{-\frac{1}{2}}P^{-1} \]

由于是在复数域,我们改写为:(注意在之前的知识1中定义的\(U^{\dagger}\)为\(S\)的特征向量组成的酉矩阵,这里我们令\(U=U^{\dagger},U^{\dagger}=U\)交换一下以便与实数域的统一。同时其中小写的\(s\)为\(S\)的对角阵)

\[S^{-\frac{1}{2}} = Us^{-\frac{1}{2}}U^{\dagger} \]

(4)对称正交化和正则正交化

Roothaan-Fall方程的矩阵形式为:

\[\mathbf{FC}=\mathbf{SC}\varepsilon \]

取酉变换矩阵\(X\):

\[XX^{\dagger}=X^{\dagger}X=I \]

注意酉变换矩阵的厄米共轭其实也是矩阵的逆:

\[X^{\dagger} = X^{-1} \]

则有:

\[FXX^{\dagger}C=SXX^{\dagger}C\varepsilon\\ \Rightarrow X^{\dagger}FXX^{\dagger}C = X^{\dagger}SXX^{\dagger}C\varepsilon\\ \Rightarrow F^{'}C^{'} = (X^{\dagger}SX)C^{'}\varepsilon \]

其中:

\[F^{'}=X^{\dagger}FX, C^{'}=X^{\dagger}C \]

对称正交化:

取酉变换矩阵为:(\(U\)也是酉矩阵)

\[X = S^{-\frac{1}{2}}=Us^{-\frac{1}{2}}U^{\dagger} \]

有:

\[X^{\dagger}=X^{-1}=(S^{-\frac{1}{2}})^{-1}=S^{\frac{1}{2}}\\ =(Us^{-\frac{1}{2}}U^{\dagger})^{-1}=(U^{\dagger})^{-1}s^{\frac{1}{2}}U^{-1}=Us^{\frac{1}{2}}U^{\dagger}\\ \Rightarrow XX^{\dagger}=Us^{-\frac{1}{2}}U^{\dagger}Us^{\frac{1}{2}}U^{\dagger}=I \]

注意上面我们把\(S^{1/2}\)即定义为\(S^{-1/2}\)的厄米共轭矩阵

则:

\[X^{\dagger}SX=(Us^{-\frac{1}{2}}U^{\dagger})^{\dagger}S(Us^{-\frac{1}{2}}U^{\dagger})\\ =Us^{-\frac{1}{2}}U^{\dagger}SUs^{-\frac{1}{2}}U^{\dagger}\\ =U^{\dagger}s^{-\frac{1}{2}}ss^{-\frac{1}{2}}U\\ =U^{\dagger}U=I \]

因此Roothaan-Fall方程变为:(变为了本征方程形式)

\[F^{'}C^{'} = C^{'}\varepsilon \]

上述的对称正交化其实是将原本不正交的基底函数\(\chi_k\)通过线性变换(酉变换)得到了一组相互正交归一化的函数,即:

\[\psi_j=\sum_{k=1}^mS_{kj}^{-1/2}\chi_k \quad (j=1,2,\cdots,m) \]

因此重叠积分的形式才发生变化:

\[S=\langle\chi_k|\chi_l\rangle\\ \Downarrow\\ S=\langle\psi_i|\psi_j\rangle=\langle \sum_kS_{ki}^{-1/2}\chi_k|\sum_lS_{lj}^{-1/2}\chi_l\rangle=\sum_{k,l}S_{ki}^{-1/2}\langle\chi_k|\chi_l\rangle S_{lj}^{-1/2}\\ =\sum_{k,l}S_{ki}^{-1/2}S_{kl} S_{lj}^{-1/2}=(S^{0})_{ij}=\delta_{ij} \]

可以看到新的函数确实是正交归一化的

正则正交化:

对称正交化要求重叠矩阵的特征值的最小值不能为极小值或为零(特征向量接近线性相关或直接就线性相关),如果遇到这种情况,则要采用正则正交化

此时我们取酉变换矩阵为:

\[X=Us^{-\frac{1}{2}} \]

同样我们可证明上述确实是酉变换矩阵:(还是注意酉矩阵的厄米共轭等于逆)

\[XX^{\dagger}=Us^{-\frac{1}{2}}(Us^{-\frac{1}{2}})^{\dagger}=Us^{-\frac{1}{2}}s^{\frac{1}{2}}U^{\dagger}=I \]

则:

\[X^{\dagger}SX=(Us^{-\frac{1}{2}})^{\dagger}SUs^{-\frac{1}{2}}=s^{\frac{1}{2}}U^{\dagger}SUs^{-\frac{1}{2}}=s^{-\frac{1}{2}}ss^{-\frac{1}{2}}=I \]

因此Roothaan-Fall方程也可以正交化

类比对称正交化,这里我们也是将原本的基底函数通过上述酉矩阵变换到了新的基底函数:(观察上面的酉矩阵)

\[\psi_i^{'}=\frac{1}{\sqrt{\lambda_i}}\sum_jc_{ji}\chi_j \]

我们同样可以证明新的基底函数也是正交归一化的:

\[\langle\psi_i^{'}|\psi_j^{'}\rangle=\frac{1}{\sqrt{\lambda_i\lambda_j}}\langle \sum_kc_{ki}\chi_k|\sum_lc_{lj}\chi_l\rangle\\ =\frac{1}{\sqrt{\lambda_i\lambda_j}}\sum_{k,l}c_{ki}^{*}S_{kl}c_{lj}=\frac{1}{\sqrt{\lambda_i\lambda_j}}\sum_kc_{ki}^*\sum_lS_{kl}c_{lj}\\ =\frac{\lambda_j}{\sqrt{\lambda_i\lambda_j}}\sum_kc_{ki}^*c_{kj}=\frac{\lambda_j}{\sqrt{\lambda_i\lambda_j}}\delta_{ij} \]

有点不明白正则的意思,上面狄拉克函数前面的就是正则吗?

由于可能存在特征值极小或者为零的情况,在正则正交化方法中,我们抛弃掉对角矩阵中的这些特征值,因此我们构建的酉矩阵并不为方阵:(假设去掉了k个,并注意我们用第二个下标表示列,因为我们假定特征向量写为列向量,同时基底函数共m维)

\[X_{ji}=U_{ji}s_{ii}^{-\frac{1}{2}}\quad(j=1,2,\cdots,m;i=1,2,\cdots,m-k) \]

通过上述方法,我们就可以实现Roothaan-Fall方程的正交化:

\[FC=SC\varepsilon \rightarrow F^{'}C^{'}=C^{'}\varepsilon\\ X^{\dagger}SX = I,F^{'}=X^{\dagger}FX, C^{'}=X^{\dagger}C \]

接下来再看SCF求解过程:

我们的目的是为了求解波函数,但是构造\(F\)算符需要波函数的积分,而解决这个矛盾的方法就是先假设初始值,然后进行SCF迭代直到收敛。就是说,我们先用已知的基底函数去展开初始的分子轨道(设置初始展开系数\(c's\)),然后用得到的分子轨道去计算\(F\)矩阵的每个元素\(F_{\beta\alpha}\)。注意这里需要计算几个积分:动能项积分、势能项积分以及重叠积分。根据分子轨道的重叠积分矩阵,我们可以得到正交归一化的酉变换矩阵,从而将\(F\)和\(C\)变换到新的\(F^{'}\)和\(C^{'}\)。将\(F^{'}\)矩阵对角化就可以求得\(\varepsilon\)和\(C^{'}\),再反求\(C\)即可求得新的\(c's\)。然后不断持续这个过程直到收敛。

这其中最重要就是\(F\)矩阵的计算,下面进行矩阵元素\(F_{\alpha\beta(i)}\)进一步推导:

\[F_{\beta\alpha(i)}=\langle\chi_{\beta}(q_i)|h_i|\chi_{\alpha}(q_i)\rangle+\sum_{j=1}^{N_{elec}}\bigg[\langle \chi_{\beta}(q_i)|J_j|\chi_{\alpha}(q_i)\rangle-\langle\chi_{\beta}(q_i)|K_j|\chi_{\alpha}(q_i)\rangle\bigg] \]

其中单电子项为:

\[h_{\beta\alpha}(q_i)=\langle\chi_{\beta}(q_i)|h_i|\chi_{\alpha}(q_i)\rangle\\ =\int\chi_{\beta}^*(q_i)(-\frac{1}{2}\nabla^2)\chi_{\alpha}(q_i)dx+\int\chi_{\beta}^*(q_i)(-\sum_{A=1}^{N_{nuc}}\frac{Z_A}{|R_A-r_i|})\chi_{\alpha}(q_i)dx = h_{\beta\alpha}(q_i) \]

我们重点关注双电子项:

库伦积分:

\[\because J_j=\int\phi_j^*(q_j)\frac{1}{r_{ij}}\phi_j(q_j)dx^{'}\\ \therefore \text{J integral} =\sum_{j=1}^{N_{elec}}\langle \chi_{\beta}(q_i)|J_j|\chi_{\alpha}(q_i)\rangle\\ =\sum_{j=1}^{N_{elec}}\int\int\chi_{\beta}^*(q_i)\phi_j^*(q_j)\frac{1}{r_{ij}}\phi_j(q_j)\chi_{\alpha}dxdx^{'}\\ \text{Let:}\quad \phi_j(q_j)^*=\sum_{\gamma}^{M_{basis}}c_{\gamma j}^*\chi_{\gamma}^*(q_j),\quad\phi_j(q_j) = \sum_{\delta}^{M_{basis}}c_{\delta j}\chi_{\delta}(q_j)\\ \therefore \text{J integral} = \sum_{j=1}^{N_{elec}}\sum_{\gamma}^{M_{basis}}\sum_{\delta}^{M_{basis}}c_{\gamma j}^*c_{\delta j}\int\int\frac{\chi_{\beta}^*\chi_{\alpha}\chi_{\gamma}^*\chi_{\delta}}{r_{ij}}dxdx^{'}\\ \text{Simplified to:}\quad = \sum_{j=1}^{N_{elec}}\sum_{\gamma}^{M_{basis}}\sum_{\delta}^{M_{basis}}c_{\gamma j}^*c_{\delta j}(\beta\alpha|\gamma\delta)\\ \text{where:} \quad (\beta\alpha|\gamma\delta)=\int\int\frac{\chi_{\beta}^*\chi_{\alpha}\chi_{\gamma}^*\chi_{\delta}}{r_{ij}}dxdx^{'} \]

交换积分:(同理,但是注意交换积分算符定义的记法稍微有点不同,下面的\(\chi_{\alpha}\)函数对应的电子变成了\(q_j\))

\[\because K_j=\int\phi_j^*(q_j)\frac{1}{r_{ij}}\phi_j(q_i)dx^{'}\\ \therefore \text{K integral} = \sum_{j=1}^{N_{elec}}\langle\chi_{\beta}(q_i)|K_j|\chi_{\alpha}(q_i)\rangle \\ =\sum_{j=1}^{N_{elec}}\int\int\chi_{\beta}^*(q_i)\phi_j^*(q_j)\frac{1}{r_{ij}}\phi_j(q_i)\chi_\alpha(q_j)dxdx^{'}\\ \text{Let:} \quad \phi_j^*(q_j)=\sum_{\gamma}^{M_{basis}}c_{\gamma j}^*\chi^*_{\gamma}(q_j),\quad \phi_j(q_i)=\sum_{\delta}^{M_{basis}}c_{\delta j}\chi_{\delta}(q_i)\\ \therefore \text{K integral} = \sum_{j=1}^{N_{elec}}\sum_{\gamma}^{M_{basis}}\sum_{\delta}^{M_{basis}}c_{\gamma j}^*c_{\delta j}\int\int\frac{\chi_{\beta}^*\chi_{\delta}\chi_{\gamma}^*\chi_{\alpha}}{r_{ij}}dxdx^{'}\\ \text{Simplied to: }\quad =\sum_{j=1}^{N_{elec}}\sum_{\gamma}^{M_{basis}}\sum_{\delta}^{M_{basis}}c_{\gamma j}^*c_{\delta j}(\beta\delta|\gamma\alpha)\\ \text{where: }\quad (\beta\delta|\gamma\alpha)=\int\int\frac{\chi_{\beta}^*\chi_{\delta}\chi_{\gamma}^*\chi_{\alpha}}{r_{ij}}dxdx^{'} \]

最终有:

\[F_{\beta\alpha(i)} = h_{\beta\alpha}(q_i) + \sum_{j=1}^{N_{elec}}\sum_{\gamma}^{M_{basis}}\sum_{\delta}^{M_{basis}}c_{\gamma j}^*c_{\delta j}[(\beta\alpha|\gamma\delta)-(\beta\delta|\gamma\alpha)] \]

进一步,我们引入密度矩阵\(P\):

\[\text{Let: }\quad P_{\gamma\delta}=\sum_{j=1}^{N_{elec}}c_{\gamma j}^*c_{\delta j}\\ \therefore F_{\beta\alpha(i)} = h_{\beta\alpha}(q_i) + \sum_{\gamma}^{M_{basis}}\sum_{\delta}^{M_{basis}}P_{\gamma\delta}[(\beta\alpha|\gamma\delta)-(\beta\delta|\gamma\alpha)]\\ \text{Here:}\quad P = \begin{pmatrix}P_{11} & P_{12} & \cdots & P_{1m}\\ P_{21} & P_{22} & \cdots & P_{2m}\\ \vdots & \vdots & \ddots & \vdots\\ P_{m1} & P_{m2} & \cdots & P_{mm} \end{pmatrix} \]

这里先跳过对\(E_{HF}\)的讨论,直接进去计算环节

3.7 HF程序

首先,正如前文提到的,由于采用高斯函数展开分子轨道能够减少积分的计算量,因此程序采用高斯函数作为基底函数,这里先对高斯函数的知识先进行学习:

关于高斯函数,这里先确定一下定义,假设高斯函数为\(\chi=Ae^{-\alpha (r-R_A)^2}\),即关于中心\(R_A=(x_A,y_A,z_A)\)对称分布,这里就理解成原子的位置。首先求解归一化系数\(A\):(复习\(\Gamma\)函数)

\[\int_{-\infty}^{+\infty}\chi^2dr=1\\ A^2\int_{-\infty}^{+\infty}e^{-2\alpha (r-R_A)^2}dr=1\\ A^2\int_{-\infty}^{+\infty}e^{-2\alpha(x-x_A)^2}dx\int_{-\infty}^{+\infty}e^{-2\alpha(y-y_A)^2}dy\int_{-\infty}^{+\infty}e^{-2\alpha(z-z_A)^2}dz=1\\ A^2\bigg(\frac{1}{\sqrt{2\alpha}}\int^{+\infty}_{-\infty}e^{-(\sqrt{2\alpha}(x-x_A))^2}d(\sqrt{2\alpha}(x-x_A)) \bigg)^3=1\\ A^2(\frac{\pi}{2\alpha})^{\frac{3}{2}}=1\\ \Rightarrow A=(\frac{2\alpha}{\pi})^{\frac{3}{4}} \]

补充知识:

定义:(\(\Gamma\)函数实现了阶乘在整个实数和复数域的扩充泛化)

\[\Gamma(z)=\int_0^{\infty}x^{z-1}\cdot e^{-x}dx \]

对定义式进行分部积分有:

\[令f=x^{z-1},g^{'}=e^{-x}\Rightarrow g=-e^{-x}\\ 由\int fg^{'}=fg-\int f^{'}g\\ \Gamma(z)=-x^{z-1}e^{-x}|_{0}^{\infty}+\int_0^{\infty}(z-1)x^{z-2}e^{-x}dx=(z-1)\int_0^{\infty}x^{z-2}e^{-x}dx=(z-1)\Gamma(z) \]

假设\(z\)为正整数,则有:

\[\Gamma(n)=(n-1)\Gamma(n-1)=(n-1)(n-2)\Gamma(n-2)=(n-1)(n-2)\cdots 3\cdot2\cdot1\cdot\Gamma(1)=(n-1)!\Gamma(1) \]

而:

\[\Gamma(1)=\int_0^{\infty}e^{-x}dx=-e^{-x}|_0^{\infty}=1 \]

故:

\[\Gamma(n)=(n-1)! \]

伽马函数也可以扩充到复数域

除了上述定义式之外,还可以写成以下公式:

\[令 x=t^2,则dx=2tdt\\ \Gamma(z)=\int_0^{\infty}x^{z-1}\cdot e^{-x}dx\\ \Rightarrow \Gamma(z)=\int_0^{\infty}t^{2z-2}\cdot e^{-t^2}\cdot2tdt\\ \Rightarrow \Gamma(z)=2\int_0^{\infty}x^{2z-1}\cdot e^{-x^2}dx\\ \]

其中当\(z=1/2\)时,有:

\[\Gamma(\frac{1}{2})=2\int_0^{\infty}e^{-x^2}dx\\ =2\sqrt{\int_0^{\infty}e^{-x^2}dx\int_0^{\infty}e^{-y^2}dy}\\ =2\sqrt{\int_0^{\infty}\int_0^{\infty}e^{-(x^2+y^2)}dxdy}\\ =2\sqrt{\int_0^{\frac{\pi}{2}}d\theta\int_0^{\infty}\rho e^{-\rho^2}d\rho}\\ =\sqrt{2\pi}\sqrt{-\frac{1}{2}\int_0^{\infty}-2\rho e^{-\rho^2}d\rho}\\ =\sqrt{2\pi}\sqrt{-\frac{1}{2} e^{-\rho^2}|_0^{\infty}}\\ =\sqrt{2\pi}\sqrt{\frac{1}{2}}\\ =\sqrt{\pi} \]

所以,归一化的高斯函数为:

\[\chi=(\frac{2\alpha}{\pi})^{\frac{3}{4}}e^{-\alpha (r-R_A)^2} \]

再看两个高斯函数的乘积:(双电子积分项中需要)

\[\chi_1\chi_2=(\frac{2\alpha}{\pi})^{\frac{3}{4}}(\frac{2\beta}{\pi})^{\frac{3}{4}}e^{-\alpha (r-R_A)^2}e^{-\beta (r-R_B)^2} \]

忽略常数项,有:

\[=e^{-\alpha (r-R_A)^2-\beta (r-R_B)^2}=e^{-\alpha r^2+2\alpha R_Ar-\alpha R_A^2-\beta r^2+2\beta R_Br-\beta R_B^2}\\ =e^{(-\alpha-\beta)r^2+2(\alpha R_A+\beta R_B)r-\alpha R_A^2-\beta R_B^2}\\ =e^{-(\alpha+\beta)[r^2-\frac{2(\alpha R_A+\beta R_B)}{\alpha+\beta}r+(\frac{\alpha R_A+\beta R_B}{\alpha+\beta})^2-(\frac{\alpha R_A+\beta R_B}{\alpha+\beta})^2+\frac{\alpha R_A^2+\beta R_B^2}{\alpha+\beta}]}\\ =e^{-(\alpha+\beta)(r-\frac{\alpha R_A+\beta R_B}{\alpha+\beta})^2}\cdot e^{\frac{(\alpha R_A+\beta R_B)^2}{\alpha+\beta}-\alpha R_A^2-\beta R_B^2}\\ =e^{-(\alpha+\beta)(r-\frac{\alpha R_A+\beta R_B}{\alpha+\beta})^2}\cdot e^{\frac{-\alpha\beta(R_A-R_B)^2}{\alpha+\beta}} \]

最终有:

\[\chi_1\chi_2 = (\frac{4\alpha\beta}{\pi^2})^{\frac{3}{4}}e^{\frac{-\alpha\beta(R_A-R_B)^2}{\alpha+\beta}} \cdot e^{-(\alpha+\beta)(r-\frac{\alpha R_A+\beta R_B}{\alpha+\beta})^2} \]

先求重叠积分矩阵:

\[S_{\alpha\beta} = \langle\chi_\alpha|\chi_{\beta}\rangle=(\frac{4\alpha\beta}{\pi^2})^{\frac{3}{4}}e^{\frac{-\alpha\beta(R_A-R_B)^2}{\alpha+\beta}}\int_{-\infty}^{\infty} e^{-p(r-R_P)^2}dr\\ \text{Here: } p = \alpha+\beta,R_P=\frac{\alpha R_A+\beta R_B}{\alpha+\beta}\\ \because \quad \int_{-\infty}^{\infty} e^{-p(r-R_P)^2}dr=\bigg[\frac{1}{\sqrt{p}}\int_{-\infty}^{\infty} e^{-[\sqrt{p}(x-X_P)]^2}d[\sqrt{p}(x-X_P)]\bigg]^3=\bigg[\frac{1}{\sqrt{p}}\Gamma(1/2)\bigg]^3=(\frac{\pi}{p})^{3/2}\\ \therefore S_{\alpha\beta} = (\frac{4\alpha\beta}{\pi^2})^{\frac{3}{4}}e^{\frac{-\alpha\beta(R_A-R_B)^2}{\alpha+\beta}}(\frac{\pi}{p})^{3/2} \]

注意分子轨道是由基底函数线性组合而成的,因此上面矩阵中还得乘以展开系数

再求动能积分矩阵:

\[T_{\alpha\beta} = \langle \chi_\alpha|-\frac{1}{2}\nabla^2|\chi_\beta\rangle = -\frac{1}{2}(\frac{4\alpha\beta}{\pi^2})^{3/4}\int e^{-\alpha|\vec{r}-\vec{R}_A|^2}\nabla^2e^{-\beta|\vec{r}-\vec{R}_B|^2}d\vec{r} \]

我们关注上面关于高斯函数的二阶导,这里使用了一个trick。虽然二阶导数是关于\(r\)的,但是我们如果我们令二阶导数关于\(R_B\)也能得到一样的结果:

\[G = e^{-\beta(x-X_B)^2}\cdot e^{-\beta(y-Y_B)^2}\cdot e^{-\beta(z-Z_B)^2} \\ \frac{\partial G}{\partial x}=-2\beta(x-X_B)e^{-\beta(x-X_B)^2}\cdot e^{-\beta(y-Y_B)^2}\cdot e^{-\beta(z-Z_B)^2} \\ \frac{\partial^2G}{\partial x^2} = 4\beta^2(x-X_B)^2e^{-\beta(x-X_B)^2}\cdot e^{-\beta(y-Y_B)^2}\cdot e^{-\beta(z-Z_B)^2} \\ \therefore \nabla_r^2G = 4\beta^2\bigg[(x-X_B)^2+(y-Y_B)^2+(z-Z_B)^2\bigg]e^{-\beta(x-X_B)^2}\cdot e^{-\beta(y-Y_B)^2}\cdot e^{-\beta(z-Z_B)^2}\\=4\beta^2(\vec{r}-\vec{R}_B)^2e^{-\beta(\vec{r}-\vec{R}_B)^2} \\ \frac{\partial G}{\partial X_B} = 2\beta(x-X_B)e^{-\beta(x-X_B)^2}\cdot e^{-\beta(y-Y_B)^2}\cdot e^{-\beta(z-Z_B)^2}\\ \frac{\partial^2 G}{\partial X_B^2} = 4\beta^2(x-X_B)^2e^{-\beta(x-X_B)^2}\cdot e^{-\beta(y-Y_B)^2}\cdot e^{-\beta(z-Z_B)^2}=4\beta^2(x-X_B)^2e^{-\beta(\vec{r}-\vec{R}_B)^2}\\ \therefore \nabla^2_{R_B} G =4\beta^2(\vec{r}-\vec{R}_B)^2e^{-\beta(\vec{r}-\vec{R}_B)^2} \]

因此我们这里很tricky的发现:

\[\nabla^2_rG = \nabla^2_{R_B}G \]

而由于右边的算符是与\(r\)无关的,因此可以提到积分外面,这样动能积分就变为先求一个重叠积分,再求导:

\[T_{\alpha\beta} = -\frac{1}{2}(\frac{4\alpha\beta}{\pi^2})^{3/4}\nabla^2_{R_B}\int e^{-\alpha|\vec{r}-\vec{R}_A|^2}e^{-\beta|\vec{r}-\vec{R}_B|^2}d\vec{r} \\ = -\frac{1}{2}(\frac{4\alpha\beta}{\pi^2})^{3/4}(\frac{\pi}{p})^{3/2}(\frac{\partial^2}{\partial X_B^2}+\frac{\partial^2}{\partial Y_B^2}+\frac{\partial^2}{\partial Z_B^2}) \bigg[e^{\frac{-\alpha\beta(R_A-R_B)^2}{\alpha+\beta}}\bigg] \]

我们看对\(X_B\)的偏导:

\[H = e^{-\frac{\alpha\beta}{\alpha+\beta}[(X_A-X_B)^2+(Y_A-Y_B)^2+(Z_A-Z_B)^2]} \\ \frac{\partial H}{\partial X_B}= \frac{2\alpha\beta}{\alpha+\beta}(X_A-X_B)e^{-\frac{\alpha\beta}{\alpha+\beta}[(X_A-X_B)^2+(Y_A-Y_B)^2+(Z_A-Z_B)^2]}\\ \frac{\partial^2 H}{\partial X_B^2} = \frac{2\alpha\beta}{\alpha+\beta}\bigg[-1+(X_A-X_B)^2\frac{2\alpha\beta}{\alpha+\beta}\bigg]e^{-\frac{\alpha\beta}{\alpha+\beta}[(X_A-X_B)^2+(Y_A-Y_B)^2+(Z_A-Z_B)^2]} \]

所以:

\[\nabla^2_{R_B}H = \frac{2\alpha\beta}{\alpha+\beta}\bigg[-3+\frac{2\alpha\beta}{\alpha+\beta}(\vec{R}_A-\vec{R}_B)^2\bigg]e^{-\frac{\alpha\beta}{\alpha+\beta}(\vec{R}_A-\vec{R}_B)^2}\\ T_{\alpha\beta} = (\frac{4\alpha\beta}{\pi^2})^{3/4}(\frac{\pi}{\alpha+\beta})^{3/2}(\frac{\alpha\beta}{\alpha+\beta})\bigg[3-\frac{2\alpha\beta}{\alpha+\beta}(\vec{R}_A-\vec{R}_B)^2\bigg]\exp[-\frac{\alpha\beta}{\alpha+\beta}(\vec{R}_A-\vec{R}_B)^2] \]

在求电子和原子核作用以及电子电子作用的积分时,可通过傅里叶变换将函数变换到复数域进行积分更加方便,为此我们先复习傅里叶变换的知识:

参考下面DR_CANup的视频,讲的特别好

纯干货数学推导_傅里叶级数与傅里叶变换_Part1_三角函数的正交性_哔哩哔哩_bilibili

Part 1. 三角函数的正交性

三角函数系:\(\{0,sin(x),cos(x),sin(2x),cos(2x),\cdots,sin(nx),cos(nx),\cdots\}\)

正交:

\[\int_{-\pi}^{\pi}sin(nx)cos(mx)dx=0 \quad (n\neq m)\\ \int_{-\pi}^{\pi}cos(nx)cos(mx)dx=0 \quad (n\neq m)\\ \int_{-\pi}^{\pi}sin(nx)sin(mx)dx=0 \quad (n\neq m)\tag{1}\\ \]

我们知道两个向量正交(垂直),则内积为零。推广到函数,若两个函数在某区间内正交,则有:

\[\int_{x_1}^{x_2}f(x)g(x)dx=0 \]

我们取\((1)\)中的一项进行证明:

\[\int_{-\pi}^{\pi}cos(nx)cos(mx)dx=0 \]

其中根据积化和差公式有:

\[cos(nx)cos(mx)=\frac{1}{2}[cos(n+m)x+cos(n-m)x] \]

则:

\[\int_{-\pi}^{\pi}cos(nx)cos(mx)dx=\frac{1}{2}[\int_{-\pi}^{\pi}cos(n+m)xdx+\int_{-\pi}^{\pi}cos(n-m)xdx]\\ =\frac{1}{2}[\frac{1}{n+m}sin(n+m)x|_{-\pi}^{\pi}+\frac{1}{n-m}sin(n-m)x|_{-\pi}^{\pi}]\\ =0 \]

如果\(m=n\),则有:

\[\int_{-\pi}^{\pi}cos(mx)cos(mx)dx=\int_{-\pi}^{\pi}cos^2(mx)dx=\int_{-\pi}^{\pi}\frac{1}{2}[1+cos(2mx)]dx\\ =\frac{1}{2}x|_{-\pi}^{\pi}+\frac{1}{2}\int_{-\pi}^{\pi}cos(0x)cos(2mx)dx=\pi \]

Part 2. 周期为\(2\pi\)的函数展开为傅里叶级数

考虑一个周期为\(2\pi\)的函数,有\(f(x)=f(x+2\pi)\)

进行傅里叶级数展开:

\[f(x)=\sum_{n=0}^{\infty}a_ncos(nx)+\sum_{n=0}^{\infty}b_nsin(nx)\tag{2} \]

教科书中给出的式子为:

\[f(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}[a_ncos(nx)+b_nsin(nx)] \]

对\((2)\)进行变换:

\[f(x)=a_0cos(0x)+\sum_{n=1}^{\infty}a_ncos(nx)+b_0sin(0x)+\sum_{n=1}^{\infty}b_nsin(nx)\\ =a_0+\sum_{n=1}^{\infty}a_ncos(nx)+\sum_{n=1}^{\infty}b_nsin(nx) \]

观察可知与教科书中的差别在\(a_0\)的系数

求\(a_0\):

\[\int_{-\pi}^{\pi}f(x)dx=\int_{-\pi}^{\pi}a_0dx+\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}a_ncos(nx)dx+\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}b_nsin(nx)dx\\\ =2a_0\pi+\sum_{n=1}^{\infty}a_n\int_{-\pi}^{\pi}cos(nx)dx+\sum_{n=1}^{\infty}b_n\int_{-\pi}^{\pi}sin(nx)dx\\ =2a_0\pi+\sum_{n=1}^{\infty}a_n\int_{-\pi}^{\pi}cos(0x)cos(nx)dx+\sum_{n=1}^{\infty}b_n\int_{-\pi}^{\pi}cos(0x)sin(nx)dx\\ =2a_0\pi\\ \Rightarrow a_0=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x)dx \]

假设令\(a_0=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)dx\),则有教科书中的形式

求\(a_n\):

\[\int_{-\pi}^{\pi}f(x)cos(mx)dx= \\\int_{-\pi}^{\pi}\frac{a_0}{2}cos(mx)dx+\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}a_ncos(nx)cos(mx)dx+\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}b_nsin(nx)cos(mx)dx\\ =0+\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}a_ncos(nx)cos(mx)dx+0\\ =\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}a_ncos(nx)cos(mx)dx \]

由于三角函数的正交性,仅当\(n=m\)时,积分不为零:

\[\int_{-\pi}^{\pi}f(x)cos(nx)dx=\int_{-\pi}^{\pi}a_ncos^2(nx)dx=a_n\pi\\ \Rightarrow a_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)cos(nx)dx \]

求\(b_n\):

\[\int_{-\pi}^{\pi}f(x)sin(mx)dx=\\ \int_{-\pi}^{\pi}\frac{a_0}{2}sin(mx)dx+\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}a_ncos(nx)sin(mx)dx+\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}b_nsin(nx)sin(mx)dx\\ =0+0+\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}b_nsin(nx)sin(mx)dx\\ =\int_{-\pi}^{\pi}\sum_{n=1}^{\infty}b_nsin(nx)sin(mx)dx \]

同理由于三角函数的正交性,仅当\(n=m\)时,积分不为零:

\[\int_{-\pi}^{\pi}f(x)sin(nx)dx=\int_{-\pi}^{\pi}b_nsin^2(nx)dx=b_n\pi\\ \Rightarrow b_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)sin(nx)dx \]

综上,我们可知周期为\(2\pi\)的函数展开为傅里叶级数:

\[f(x)=f(x+2\pi)\\ f(x)=\frac{a_0}{2}+\sum_{n=1}^{\infty}[a_ncos(nx)+b_nsin(nx)]\\ a_0=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)dx\\ a_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)cos(nx)dx\\ b_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)sin(nx)dx \]

Part 3. 周期为\(2L\)的函数展开为傅里叶级数

已知:\(f(t)=f(t+2L)\)

换元:

\[x=\frac{\pi}{L}t\\ t=\frac{L}{\pi}x \]

则:

\[f(t)=f(\frac{L}{\pi}x)=g(x) \]

其中\(g(x)\)是以\(2\pi\)为周期的函数

有:

\[x=\frac{\pi}{L}t\\ cos(nx)=cos(\frac{n\pi}{L}t)\\ sin(nx)=sin(\frac{n\pi}{L}t)\\ g(x)=f(t)\\ \int_{-\pi}^{\pi}dx=\int_{-L}^Ld(\frac{\pi}{L}t)=\frac{\pi}{L}\int_{-L}^Ldt\\ \frac{1}{\pi}\int_{-\pi}^{\pi}dx=\frac{1}{L}\int_{-L}^Ldt \]

对照part 2中的傅里叶级数,有:

\[f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}[a_ncos(\frac{n\pi}{L}t)+b_nsin(\frac{n\pi}{L}t)]\\ a_0=\frac{1}{L}\int_{-L}^{L}f(t)dt\\ a_n=\frac{1}{L}\int_{-L}^{L}f(t)cos(\frac{n\pi}{L}t)dt\\ b_n=\frac{1}{L}\int_{-L}^{L}f(t)cos(\frac{n\pi}{L}t)dt \]

令\(\omega=\frac{2\pi}{T},T=2L\) ,则有:(\(f(t)=f(t+T)\))

\[\int_{-L}^Ldt=\int_0^{2L}dt=\int_0^Tdt\\ f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}[a_ncos(n\omega t)+b_nsin(n\omega t)]\\ a_0=\frac{2}{T}\int_{0}^{T}f(t)dt\\ a_n=\frac{2}{T}\int_{0}^{T}f(t)cos(n\omega t)dt\\ b_n=\frac{2}{T}\int_{0}^{T}f(t)sin(n\omega t)dt \]

Part 4. 傅里叶级数的复数形式

由欧拉公式有:

\[cos(\theta)=\frac{1}{2}(e^{i\theta}+e^{-i\theta})\\ sin(\theta)=-\frac{1}{2}i(e^{i\theta}-e^{-i\theta}) \]

代入得:

\[f(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty}[a_n\cdot\frac{1}{2}(e^{in\omega t}+e^{-in\omega t})-b_n\cdot\frac{1}{2}i(e^{in\omega t}-e^{-in\omega t})]\\ =\frac{a_0}{2}+\sum_{n=1}^{\infty}[\frac{a_n-ib_n}{2}e^{in\omega t}+\frac{a_n+ib_n}{2}e^{-in\omega t}]\\ =\frac{a_0}{2}+\sum_{n=1}^{\infty}\frac{a_n-ib_n}{2}e^{in\omega t}+\sum_{n=1}^{\infty}\frac{a_n+ib_n}{2}e^{-in\omega t}\\ =\frac{a_0}{2}+\sum_{n=1}^{\infty}\frac{a_n-ib_n}{2}e^{in\omega t}+\sum_{n=-\infty}^{-1}\frac{a_n+ib_n}{2}e^{in\omega t}(令n\rightarrow -n )\\ =\sum_{n=0}^0\frac{a_0}{2}e^{in\omega t}+\sum_{n=1}^{\infty}\frac{a_n-ib_n}{2}e^{in\omega t}+\sum_{n=-\infty}^{-1}\frac{a_{-n}+ib_{-n}}{2}e^{in\omega t}\\ =\sum_{-\infty}^{\infty}C_ne^{in\omega t} \]

故:

\[C_n=\begin{cases}\frac{a_0}{2}, & n=0\\ \frac{a_n-ib_n}{2}, & n=1,2,3,\cdots \\ \frac{a_{-n}+ib_{-n}}{2}, & n=-1,-2,-3,\cdots \end{cases} \]

接下来代入前述的周期为\(T\)的傅里叶级数的系数表达式:

\[a_0=\frac{2}{T}\int_{0}^{T}f(t)dt\\ a_n=\frac{2}{T}\int_{0}^{T}f(t)cos(n\omega t)dt\\ b_n=\frac{2}{T}\int_{0}^{T}f(t)sin(n\omega t)dt \]

则有:

\[C_n=\frac{a_0}{2}=\frac{1}{T}\int_{0}^{T}f(t)dt=\frac{1}{T}\int_0^Tf(t)e^{-in\omega t}dt \quad (n=0)\\ C_n=\frac{1}{2}[\frac{2}{T}\int_{0}^{T}f(t)cos(n\omega t)dt-i\frac{2}{T}\int_{0}^{T}f(t)sin(n\omega t)dt]\\ =\frac{1}{T}\int_0^Tf(t)[cos(n\omega t)-isin(n\omega t)]dt\\ =\frac{1}{T}\int_0^Tf(t)e^{-in\omega t}dt \quad (n=1,2,3,\cdots)\\ C_n=\frac{1}{2}[\frac{2}{T}\int_{0}^{T}f(t)cos(-n\omega t)dt+i\frac{2}{T}\int_{0}^{T}f(t)sin(-n\omega t)dt]\\ =\frac{1}{T}\int_0^Tf(t)e^{-in\omega t}dt \quad (n=-1,-2,-3,\cdots) \]

综上有:

\[f(t)=f(t+T)\\ f(t)=\sum_{-\infty}^{\infty}C_ne^{in\omega t}\\ C_n=\frac{1}{T}\int_0^Tf(t)e^{-in\omega t}dt \]

Part 5. 傅里叶变换

前面的周期性函数有下标:

\[f_T(t)=f(t+T)\\ f_T(t)=\sum_{-\infty}^{\infty}C_ne^{in\omega_0 t}\\ C_n=\frac{1}{T}\int_0^Tf_T(t)e^{-in\omega_0 t}dt \]

观察可知,在时域空间的函数\(f_T(t)\)可以转化为频域空间不同频率值的加和形式,即:

\[f_T(t)=\cdots +C_{-1}e^{-i\omega_0t}+C_0e^0+C_1e^{i\omega_0t}+C_1e^{2i\omega_0t}+\cdots \]

考虑当周期\(T\rightarrow \infty\)时的一般形式,则\(\lim_{T\rightarrow\infty}f_T(t)=f(t)\)

取频域横坐标的间隔为\(\Delta \omega=(n+1)\omega_0-n\omega_0=\omega_0=2\pi/T\)

当\(T\rightarrow \infty\)时,\(\Delta \omega \rightarrow 0\),此时离散的加和就变为连续,在频域空间得到对应的曲线,此时横坐标由\(n\omega_0\)变为\(\omega\)

代入系数得:

\[f_T(t)=\sum_{-\infty}^{\infty}C_ne^{in\omega_0 t}=\sum_{-\infty}^{\infty}\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f_T(t)e^{-in\omega_0 t}dte^{in\omega_0 t}\\ \because \frac{1}{T}=\frac{\Delta \omega}{2\pi}\\ f_T(t)=\sum_{-\infty}^{\infty}\frac{\Delta \omega}{2\pi}\int_{-\frac{T}{2}}^{\frac{T}{2}}f_T(t)e^{-in\omega_0 t}dte^{in\omega_0 t} \]

令\(T\rightarrow \infty\),则有:

\[\int_{-\frac{T}{2}}^{\frac{T}{2}}dt\rightarrow\int_{-\infty}^{\infty}dt\\ n\omega_0\rightarrow\omega\\ \Rightarrow\sum_{-\infty}^{\infty}\Delta \omega=\int_{-\infty}^{\infty}d\omega \]

则:

\[f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(t)e^{-i\omega t}dte^{i\omega t}d\omega \]

其中:

\[F(\omega)=\int_{-\infty}^{\infty}f(t)e^{-i\omega t}dt \]

称为傅里叶变换

而:

\[f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{i\omega t}d\omega \]

称为傅里叶变换的逆变换

在电子积分中\(t\rightarrow \vec{r},\omega \rightarrow \vec{k}\),则有:

\[F(\vec{k}) = \int_{-\infty}^{\infty}f(\vec{r})e^{-i\vec{k}\cdot\vec{r}}d\vec{r}\\ f(\vec{r}) = (\frac{1}{2\pi})^{3}\int F(\vec{k})e^{i\vec{k}\cdot\vec{r}}d\vec{k} \]

电子积分中涉及多个函数,下面给出各种函数的傅里叶变换:

\(\frac{1}{r}\):

我们一般对可积函数做FT变换,但是这里由于\(\frac{1}{r}\)是不可积的,我们需要将库伦势先改成汤川势,即:

\[\frac{1}{r} \rightarrow \frac{e^{-\alpha r}}{r} \]

先求解汤川势的傅里叶变换,然后令\(\alpha \rightarrow 0\)即可得库伦势的傅里叶变换

注意由于被积函数是球对称函数,只和\(r\)有关,所以旋转不会改变积分值,因此可以将坐标系旋转,令z轴指向\(\vec{k}\)的方向,因而有:(转成极坐标时,雅可比行列式的值为\(r^2\sin\theta\))

\[\mathcal{F}(\frac{e^{-\alpha r}}{r})=\int \frac{e^{-\alpha r}}{r} e^{-i\vec{k}\cdot \vec{r}}d\vec{r}^3=\int_0^{2\pi}d\phi\int_{0}^{\infty}\frac{e^{-\alpha r}}{r}r^2dr\int_0^{\pi}e^{-ikrcos\theta}sin \theta d\theta\\ =2\pi\int_0^{\infty}dr\frac{e^{-\alpha r}}{r}r^2(\frac{1}{ikr}e^{-ikrcos\theta}|_0^\pi)\\ =\frac{2\pi}{ik}\int_0^{\infty}e^{-\alpha r}(e^{ikr}-e^{-ikr})dr\\ =\frac{2\pi}{ik}\bigg[\int_0^{\infty}e^{(-\alpha+ik)r}dr - \int_0^{\infty}e^{(-\alpha-ik)r}dr\bigg]\\ =\frac{2\pi}{ik}\bigg[\frac{1}{-\alpha+ik}e^{(-\alpha+ik)r}\bigg|_0^\infty - \frac{1}{-\alpha-ik}e^{(-\alpha-ik)r}\bigg|_0^\infty\bigg]\\ =\frac{2\pi}{ik}\bigg[\frac{-(\alpha+ik)}{\alpha^2+k^2}e^{(-\alpha+ik)r}+\frac{\alpha-ik}{\alpha^2+k^2}e^{(-\alpha-ik)r}\bigg]_0^{\infty}\\ =\frac{2\pi}{ik(\alpha^2+k^2)}\bigg[-\alpha e^{(-\alpha+ik)r}-ike^{(-\alpha+ik)r}+\alpha e^{(-\alpha-ik)r}-ike^{(-\alpha-ik)r}\bigg]_0^{\infty}\\ \]

因为当\(r\rightarrow\infty\)时,若是\(e^{ik r}|_{x\rightarrow \infty}\),由于欧拉公式展开后是sin和cos函数,因此结果为有限值,因此我们只看\(e^{-\alpha r}|_{x\rightarrow \infty}\)项即可,结果为零,故:

\[=\frac{2\pi}{ik(\alpha^2+k^2)}\bigg[\alpha+ik-\alpha+ik\bigg]=\frac{4\pi}{\alpha^2+k^2} \]

当\(\alpha\rightarrow0\)时,得:

\[\mathcal{F}(\frac{1}{r})=\frac{4\pi}{k^2} \]

\(e^{-\alpha r^2}\):

\[\mathcal{F}(e^{-\alpha r^2})= \int e^{-\alpha r^2}e^{-i\vec{k}\cdot\vec{r}}d\vec{r}^3\\ =\int e^{-\alpha(x^2+y^2+z^2)}e^{-i(k_xx+k_yy+k_zz)}dxdydz\\ =\int e^{-\alpha x^2}e^{-ik_xx}dx \int e^{-\alpha y^2}e^{-ik_yy}dy \int e^{-\alpha z^2}e^{-ik_zz}dz\\ =\int_{-\infty}^{\infty} e^{-\alpha(x+\frac{ik_x}{2\alpha})^2-\frac{k_x^2}{4\alpha}}dx \int_{-\infty}^{\infty} e^{-\alpha(y+\frac{ik_y}{2\alpha})^2-\frac{k_y^2}{4\alpha}}dy \int_{-\infty}^{\infty} e^{-\alpha(z+\frac{ik_z}{2\alpha})^2-\frac{k_z^2}{4\alpha}}dz\\ \text{x integral} =e^{-\frac{k_x^2}{4\alpha}}\int_{-\infty}^{\infty}\frac{1}{\sqrt{\alpha}}e^{-[\sqrt{\alpha}(x+\frac{ik_x}{2\alpha})]^2}d(\sqrt{\alpha}(x+\frac{ik_x}{2\alpha})) =e^{-\frac{k_x^2}{4\alpha}}\frac{1}{\sqrt{\alpha}}\Gamma(\frac{1}{2})\\ \because \Gamma(\frac{1}{2})=\sqrt{\pi} \\ \therefore \text{x integral}=(\frac{\pi}{\alpha})^{\frac{1}{2}}e^{-\frac{k_x^2}{4\alpha}}\\ \therefore \text{whole integral:}\quad\mathcal{F}(e^{-\alpha r^2})=(\frac{\pi}{\alpha})^{\frac{3}{2}}e^{-\frac{k_x^2+k_y^2+k_z^2}{4\alpha}}=(\frac{\pi}{\alpha})^{\frac{3}{2}}e^{-\frac{k^2}{4\alpha}} \]

\(\delta(r)\):

狄拉克delta函数定义为:(三维同理)

\[\delta(x)= \begin{cases}+\infty, & x=0 \\ 0,&x\neq0 \end{cases}\\ \int_{-\infty}^{\infty}\delta(x)dx=1 \]

所以有:

\[\mathcal{F}(\delta(r))=\int \delta(r)e^{-i\vec{k}\cdot\vec{r}}d\vec{r}^3\\ \because \text{when }r\neq 0,\delta(r)=0\\ \therefore \mathcal{F}(\delta(r))=\int\delta(r)dr\bigg|_{r\neq 0}+\int\delta(r)dr\bigg|_{r=0}=1 \]

下面我们考虑电子和原子核库伦作用的积分:

假设高斯函数为:(不进行归一化)

\[\chi_1 = e^{-\alpha (\vec{r}-\vec{R}_A)^2}\\ \chi_2 = e^{-\beta(\vec{r}-\vec{R}_B)^2} \]

则:(考虑某个电子的位置\(r_1\))

\[\langle\chi_1|-\frac{Z_A}{r_{1A}}|\chi_2\rangle = -Z_A\int d\vec{r}_1e^{-\alpha|\vec{r}_1-\vec{R}_A|^2}|\vec{r}_1-\vec{R}_A|^{-1}e^{-\beta|\vec{r}_1-\vec{R}_B|^2}\\ \text{Combine two Gaussians to a single one:}\\ =-Z_A\tilde{K}\int d\vec{r}_1e^{-p|\vec{r}_1-\vec{R}_P|^2}|\vec{r}_1-\vec{R}_A|^{-1}\\ \text{Here: }\tilde{K} = e^{\frac{-\alpha\beta(\vec{R}_A-\vec{R}_B)^2}{\alpha+\beta}},\vec{R}_p=\frac{\alpha \vec{R}_A+\beta \vec{R}_B}{\alpha+\beta},p=\alpha+\beta\\ \text{Then substitude the FT of Gaussian and 1/r into the integral:}\\ =-Z_A\tilde{K}\int d\vec{r}_1\bigg[(\frac{1}{2\pi})^3\int d\vec{k}_1(\frac{\pi}{p})^{3/2}e^{-\vec{k}_1^2/4p}\cdot e^{i\vec{k}_1\cdot(\vec{r}_1-\vec{R}_P)}\bigg]\bigg[(\frac{1}{2\pi})^3\int d\vec{k}_2\frac{4\pi}{\vec{k}_2^2}\cdot e^{i\vec{k}_2\cdot(\vec{r}_1-\vec{R}_A)}\bigg]\\ =-Z_A\tilde{K}(2\pi)^{-6}\cdot 4\pi\cdot(\frac{\pi}{p})^{3/2}\int\int\int d\vec{r}_1 d\vec{k}_1 d\vec{k}_2\cdot \bigg[e^{-\vec{k}_1^2/4p}\vec{k}_2^{-2}\bigg]\bigg[e^{i\vec{k}_1\cdot(\vec{r}_1-\vec{R}_P)}e^{i\vec{k}_2\cdot(\vec{r}_1-\vec{R}_A)}\bigg]\\ =-Z_A\tilde{K}(2\pi)^{-6}\cdot 4\pi\cdot(\frac{\pi}{p})^{3/2}\int\int\int d\vec{r}_1 d\vec{k}_1 d\vec{k}_2\cdot \bigg[e^{-\vec{k}_1^2/4p}\vec{k}_2^{-2}\bigg]\bigg[e^{-i\vec{k}_1\vec{R}_P}e^{-i\vec{k}_2\vec{R}_A}e^{i\vec{r}_1(\vec{k}_1+\vec{k}_2)}\bigg]\\ \]

由于delta函数有:

\[\delta(\vec{r}_1-\vec{r}_2)=(\frac{1}{2\pi})^3\int d\vec{k}e^{i\vec{k}(\vec{r}_1-\vec{r}_2)} \]

对应\(r\)为复数情况有:

\[\delta(\vec{k}_1-\vec{k}_2)=(\frac{1}{2\pi})^3\int d\vec{r}e^{i\vec{r}(\vec{k}_1-\vec{k}_2)} \]

因此上述积分中:

\[(\frac{1}{2\pi})^3\int d\vec{r}_1e^{i\vec{r}_1(\vec{k}_1+\vec{k}_2)}=\delta(\vec{k}_1+\vec{k}_2) \]

则积分变为:

\[=-Z_A\tilde{K}(\frac{1}{2\pi^2})\cdot(\frac{\pi}{p})^{3/2}\int\int d\vec{k}_1 d\vec{k}_2\cdot \bigg[e^{-\vec{k}_1^2/4p}\vec{k}_2^{-2}\bigg]\bigg[e^{-i\vec{k}_1\vec{R}_P}e^{-i\vec{k}_2\vec{R}_A}\bigg]\delta(\vec{k}_1+\vec{k}_2) \]

令\(\vec{k}_1=\vec{k},\vec{k}_2=-\vec{k}_1=-\vec{k}\),则\(d\vec{k}_1=d\vec{k},d\vec{k}_2=-d\vec{k}\)

有:(注意积分上下限)

\[\int_{-\infty}^{\infty}\delta(\vec{k}_1+\vec{k}_2)d\vec{k}_2=-\int_{\infty}^{-\infty}\delta(0)d\vec{k}=\int_{-\infty}^{\infty}\delta(0)d\vec{k}=1 \]

则积分变为:

\[=-Z_A\tilde{K}(\frac{1}{2\pi^2})\cdot(\frac{\pi}{p})^{3/2}\int d\vec{k} \cdot \bigg[e^{-\vec{k}^2/4p}\vec{k}^{-2}\bigg]\bigg[e^{-i\vec{k}\vec{R}_P}e^{i\vec{k}\vec{R}_A}\bigg]\\ =-Z_A\tilde{K}(\frac{1}{2\pi^2})\cdot(\frac{\pi}{p})^{3/2}\int d\vec{k} \cdot e^{-\vec{k}^2/4p}\vec{k}^{-2}e^{-i\vec{k}(\vec{R}_P-\vec{R}_A)} \]

由于被积函数只是\(\vec{k}\)的函数,可旋转坐标系,使得\(z\)轴指向\(\vec{k}\)方向,则有\(\vec{k}\cdot(\vec{R}_P-\vec{R}_A)=k|\vec{R}_P-\vec{R}_A|\cos\theta\)

球坐标积分为:

若\(|\vec{R}_P-\vec{R}_A| \neq 0\):

\[=-Z_A\tilde{K}(\frac{1}{2\pi^2})\cdot(\frac{\pi}{p})^{3/2}\int_0^{2\pi}d\phi\int_0^{\infty}dk\cdot e^{-k^2/4p}k^{-2}\cdot k^2\int_0^{\pi}d\theta e^{-ik|\vec{R}_P-\vec{R}_A|\cos\theta}\sin\theta \\ =-Z_A\tilde{K}\frac{1}{\pi}(\frac{\pi}{p})^{3/2}\int_0^{\infty}dk\cdot e^{-k^2/4p}\cdot\frac{1}{ik|\vec{R}_P-\vec{R}_A|}e^{-ik|\vec{R}_P-\vec{R}_A|\cos\theta}\bigg|_0^{\pi}\\ =-Z_A\tilde{K}\frac{1}{\pi}(\frac{\pi}{p})^{3/2}\int_0^{\infty}dk\cdot e^{-k^2/4p}\cdot\frac{1}{ik|\vec{R}_P-\vec{R}_A|}(e^{ik|\vec{R}_P-\vec{R}_A|}-e^{-ik|\vec{R}_P-\vec{R}_A|})\\ =-Z_A\tilde{K}\frac{1}{\pi}(\frac{\pi}{p})^{3/2}\int_0^{\infty}dk\cdot e^{-k^2/4p}\cdot\frac{1}{ik|\vec{R}_P-\vec{R}_A|}2i\sin(k|\vec{R}_P-\vec{R}_A|)\\ =-2Z_A\tilde{K}(\pi|\vec{R}_P-\vec{R}_A|)^{-1}(\frac{\pi}{p})^{3/2}\int_0^{\infty}dk\cdot e^{-k^2/4p}\cdot k^{-1}\cdot\sin(k|\vec{R}_P-\vec{R}_A|) \]

若\(|\vec{R}_P-\vec{R}_A| = 0\):

\[=-Z_A\tilde{K}(\frac{1}{2\pi^2})\cdot(\frac{\pi}{p})^{3/2}\int_0^{2\pi}d\phi\int_0^{\infty}dk\cdot e^{-k^2/4p}k^{-2}\cdot k^2\int_0^{\pi}d\theta e^{-ik|\vec{R}_P-\vec{R}_A|\cos\theta}\sin\theta\\ =-Z_A\tilde{K}(\frac{1}{2\pi^2})\cdot(\frac{\pi}{p})^{3/2}\cdot 2\pi\cdot\int_0^{\pi}\sin\theta d\theta\cdot \int_0^{\infty}e^{-k^2/4p}dk\\ =-2Z_A\tilde{K}\frac{1}{\pi}(\frac{\pi}{p})^{3/2}\sqrt{4p}\int_0^{\infty} e^{-(k/\sqrt{4p})^2}d(k/\sqrt{4p})\\ =-4Z_A\tilde{K}\frac{1}{\pi}(\frac{\pi}{p})^{3/2}p^{1/2}\frac{1}{2}\pi^{1/2}\\ =-2Z_A\tilde{K}(\frac{\pi}{p}) \]

令\(N= -2Z_A\tilde{K}(\pi|\vec{R}_P-\vec{R}_A|)^{-1}(\frac{\pi}{p})^{3/2}\),则:

\[=N\int_0^{\infty}dk\cdot e^{-k^2/4p}\cdot k^{-1}\cdot\sin(k|\vec{R}_P-\vec{R}_A|) \]

为了计算上面的积分,我们先考虑普遍的形式:(注意这里\(x=|\vec{R}_P-\vec{R}_A|\))

\[I(x)=\int_0^{\infty}dk\cdot e^{-ak^2}\cdot k^{-1}\cdot \sin kx\\ I(0)=0 \]

对应的导数为:(是个偶函数)

\[I^{'}(x)=\int_0^{\infty}dk\cdot e^{-ak^2}\cos kx \]

则有:

\[I(x)-I(0)=I(x)=\int_0^x dyI^{'}(y) \]

因为:

\[I^{'}(x)=\frac{1}{2}\int_{-\infty}^{\infty}dk\cdot e^{-ak^2}\cos kx=\frac{1}{2}\mathcal{Re}\bigg[\int_{-\infty}^{\infty}dk\cdot e^{-ak^2}e^{ikx}\bigg]\\ =\frac{1}{2}\mathcal{Re}\bigg[\int_{-\infty}^{\infty}dk\cdot e^{-[(\sqrt{a}k)^2-ixk+(\frac{ix}{2\sqrt{a}})^2-(\frac{ix}{2\sqrt{a}})^2]}\bigg]\\ =\frac{1}{2}e^{-x^2/(4a)}\mathcal{Re}\bigg[ \int_{-\infty}^{\infty}dk\cdot e^{-(\sqrt{a}k-\frac{ix}{2\sqrt{a}})^2} \bigg] \]

令\(u=\sqrt{a}k-\frac{ix}{2\sqrt{a}}\),则\(dk=\frac{1}{\sqrt{a}}du\)

\[I^{'}(x)=\frac{1}{2}e^{-x^2/(4a)}\frac{1}{\sqrt{a}}\int_{-\infty}^{\infty}e^{-u^2}du=\frac{1}{2}e^{-x^2/(4a)}\frac{1}{\sqrt{a}}\Gamma(\frac{1}{2})=\frac{1}{2}(\pi/a)^{1/2}e^{-x^2/(4a)} \]

所以有:

\[I(x)=\frac{1}{2}(\pi/a)^{1/2}\int_0^x e^{-y^2/(4a)}dy \]

因此积分为:(\(\frac{1}{4p}=a,|\vec{R}_P-\vec{R}_A|=x\))

\[=-2Z_A\tilde{K}(\pi|\vec{R}_P-\vec{R}_A|)^{-1}(\frac{\pi}{p})^{3/2}\frac{1}{2}(4p\pi)^{1/2}\int_0^{|\vec{R}_P-\vec{R}_A|}e^{-py^2}dy\\ =-2\pi Z_A\tilde{K}(p|\vec{R}_P-\vec{R}_A|)^{-1}\int_0^{|\vec{R}_P-\vec{R}_A|}e^{-py^2}dy\\ \text{Let: } \sqrt{p}y\rightarrow t\Rightarrow dy=\frac{1}{\sqrt{p}}dt\\ \therefore \quad = -2\pi Z_A\tilde{K}(p|\vec{R}_P-\vec{R}_A|)^{-1}\frac{1}{\sqrt{p}}\int_0^{\sqrt{p}|\vec{R}_P-\vec{R}_A|}e^{-t^2}dt\\ =-2\pi Z_A\tilde{K}p^{-1}(p^{1/2}|\vec{R}_P-\vec{R}_A|)^{-1}\int_0^{p^{1/2}|\vec{R}_P-\vec{R}_A|}e^{-y^2}dy \quad (t\rightarrow y ) \]

定义\(F_0(t)\)函数:

\[F_0(t)=t^{-1/2}\int_0^{t^{1/2}}e^{-y^2}dy \]

而误差函数\(erf(x)\)定义为:

\[erf(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-\eta^2}d\eta \]

则有:

\[F_0(t)=\frac{1}{2}(\frac{\pi}{t})^{1/2}erf(t^{1/2}) \]

令\(t=p|\vec{R}_p-\vec{R}_A|^2\),则\(t^{1/2}=p^{1/2}|\vec{R}_p-\vec{R}_A|\)

积分可变为:

\[= -2\pi Z_A\tilde{K}p^{-1}t^{-1/2}\int_0^{t^{1/2}}e^{-y^2}dy\\ \]

最终积分为:

\[\langle\chi_1|-\frac{Z_A}{r_{1A}}|\chi_2\rangle = -\frac{2\pi Z_A}{\alpha+\beta} \exp(-\frac{\alpha\beta |\vec{R}_A-\vec{R}_B|^2}{\alpha+\beta})F_0(t)\\ \text{Here we let: }\\ t=p|\vec{R}_P-\vec{R}_A|^2=(\alpha+\beta)|\vec{R}_P-\vec{R}_A|^2\\ F_0(t)=\frac{1}{2}(\frac{\pi}{t})^{1/2}erf(t^{1/2})\\ \text{when }t=0: \text{let }F_0(t) = 1 \]

误差函数可以调库计算得到,因此积分得解

再计算双电子积分:

\[\langle\chi_1\chi_2|\frac{1}{r_{12}}|\chi_3\chi_4\rangle = \int\int d\vec{r}_1d\vec{r}_2e^{-\alpha|\vec{r}_1-\vec{R}_A|^2}e^{-\beta|\vec{r}_1-\vec{R}_B|^2}\frac{1}{r_{12}}e^{-\gamma|\vec{r}_2-\vec{R}_C|^2}e^{-\delta|\vec{r}_2-\vec{R}_D|}\\ \text{Combine two gaussians to a single one:}\\ =\exp[-\frac{\alpha\beta|\vec{R}_A-\vec{R}_B|^2}{\alpha+\beta}-\frac{\gamma\delta|\vec{R}_C-\vec{R}_D|^2}{\gamma+\delta}]\times \int\int d\vec{r}_1d\vec{r}_2e^{-p|\vec{r}_1-\vec{R}_P|^2}\frac{1}{r_{12}}e^{-q|\vec{r}_2-\vec{R}_Q|^2}\\ \text{Let: } M = \exp[-\frac{\alpha\beta|\vec{R}_A-\vec{R}_B|^2}{\alpha+\beta}-\frac{\gamma\delta|\vec{R}_C-\vec{R}_D|^2}{\gamma+\delta}]\\ \text{substitute the FT into the integral:}\\ =M(\frac{1}{2\pi})^9\int\int\int\int\int d\vec{r}_1d\vec{r}_2d\vec{k}_1d\vec{k}_2d\vec{k}_3(\frac{\pi}{p})^{3/2}e^{-k_1^2/4p}e^{i\vec{k}_1\cdot(\vec{r}_1-\vec{R}_P)}\times \\ \frac{4\pi}{k_2^2}e^{i\vec{k}_2\cdot(\vec{r}_1-\vec{r}_2)}\times (\frac{\pi}{q})^{3/2}e^{-k_3^2/4q}e^{i\vec{k}_3\cdot(\vec{r}_2-\vec{R}_Q)}\\ \text{The integral signs are abbreviated to one sign for simplicity:}\\ =4\pi M(2\pi)^{-9}(\frac{\pi^2}{pq})^{3/2}\int d\vec{r}_1d\vec{r}_2d\vec{k}_1d\vec{k}_2d\vec{k}_3 e^{-k_1^2/4p}e^{-k_3^2/4q}\times k_2^{-2}\times \\e^{-i\vec{k}_1 \vec{R_P}}e^{-i\vec{k}_3\vec{R}_Q}e^{i\vec{r}_1(\vec{k}_1+\vec{k}_2)}e^{-i\vec{r}_2(\vec{k}_3-\vec{k}_2)}\\ \text{use the definition of delta function, we have:}\\ =4\pi M(2\pi)^{-3}(\frac{\pi^2}{pq})^{3/2}\int d\vec{k}_1d\vec{k}_2d\vec{k}_3 e^{-k_1^2/4p}e^{-k_3^2/4q}\times k_2^{-2}\times \\ e^{-i\vec{k}_1 \vec{R_P}}e^{-i\vec{k}_3\vec{R}_Q}\delta(\vec{k}_1+\vec{k}_2)\delta(\vec{k}_3-\vec{k}_2)\\ \text{Let: } \quad \vec{k}_1=-\vec{k}_2=-\vec{k}, \vec{k}_3=\vec{k}_2=\vec{k}\\ \therefore \quad = \frac{M}{2\pi^2}(\frac{\pi^2}{pq})^{3/2}\int d\vec{k}e^{-(p+q)k^2/(4pq)}k^{-2}e^{i\vec{k}\cdot(\vec{R}_P-\vec{R}_Q)} \]

由于被积函数只是\(\vec{k}\)的函数,可旋转坐标系,使得\(z\)轴指向\(\vec{k}\)方向,则有\(\vec{k}\cdot(\vec{R}_P-\vec{R}_Q)=k|\vec{R}_P-\vec{R}_Q|\cos\theta\)

球坐标积分为:

若\(|\vec{R}_P-\vec{R}_Q| \neq 0\):

\[= \frac{M}{2\pi^2}(\frac{\pi^2}{pq})^{3/2}\int_0^{2\pi}d\phi\int_0^{\infty}dk\cdot e^{-(p+q)k^2/(4pq)}k^{-2}\cdot k^2 \int_0^{\pi}e^{ik|\vec{R}_P-\vec{R}_Q|\cos\theta}\sin\theta d\theta\\ =\frac{M}{\pi}(\frac{\pi^2}{pq})^{3/2}\int_0^{\infty}dk\cdot e^{-(p+q)k^2/(4pq)}\cdot \frac{-1}{ik|\vec{R}_P-\vec{R}_Q|}e^{ik|\vec{R}_P-\vec{R}_Q|\cos\theta}\bigg|_0^{\pi}\\ =M\pi^2(pq)^{-3/2}\int_0^{\infty}dk\cdot e^{-(p+q)k^2/(4pq)}\cdot \frac{-1}{ik|\vec{R}_P-\vec{R}_Q|}\cdot[-2i\sin(k|\vec{R}_P-\vec{R}_Q|)]\\ =2M\pi^2(pq)^{-3/2}|\vec{R}_P-\vec{R}_Q|^{-1}\int_0^{\infty} e^{-(p+q)k^2/(4pq)}\cdot k^{-1}\cdot \sin(k|\vec{R}_P-\vec{R}_Q|)dk \]

若\(|\vec{R}_P-\vec{R}_Q| = 0\):

\[=\frac{M}{2\pi^2}(\frac{\pi^2}{pq})^{3/2}\int_0^{2\pi}d\phi\int_0^{\infty}dk\cdot e^{-(p+q)k^2/(4pq)}k^{-2}\cdot k^2 \int_0^{\pi}e^{ik|\vec{R}_P-\vec{R}_Q|\cos\theta}\sin\theta d\theta\\ =\frac{M}{2\pi^2}(\frac{\pi^2}{pq})^{3/2}\cdot 2\pi\cdot\int_0^{\pi}\sin\theta d\theta \cdot\frac{2\sqrt{pq}}{\sqrt{p+q}}\int_0^{\infty}e^{-(\sqrt{(p+q)/(4pq)}k)^2}d(\sqrt{(p+q)/(4pq)}k)\\ =\frac{4M}{\pi}(\frac{\pi^2}{pq})^{3/2}\frac{(pq)^{1/2}}{(p+q)^{1/2}}\frac{1}{2}\pi^{1/2}\\ =2M\pi^{5/2}(p+q)^{-1/2}(pq)^{-1} \]

令:(\(a=\frac{p+q}{4pq},x=|\vec{R}_P-\vec{R}_Q|\))

\[\int_0^{\infty} e^{-(p+q)k^2/(4pq)}\cdot k^{-1}\cdot \sin(k|\vec{R}_P-\vec{R}_Q|)dk=\int_0^{\infty}e^{-ak^2}k^{-1}\sin(kx)dk=I(x) \]

可以看到,上述积分\(I(x)\)的形式和库伦势积分的形式是一样的,因此直接套用即可,即:

\[I(x)=\frac{1}{2}(\pi/a)^{1/2}\int_0^x e^{-y^2/(4a)}dy=\frac{1}{2}\pi^{1/2}(\frac{4pq}{p+q})^{1/2}\int_0^{|\vec{R}_P-\vec{R}_Q|}e^{-pqy^2/(p+q)}dy \]

积分即为:

\[=2M\pi^2(pq)^{-3/2}|\vec{R}_P-\vec{R}_Q|^{-1}I(x)\\ =2M\pi^2(pq)^{-3/2}|\vec{R}_P-\vec{R}_Q|^{-1}\frac{1}{2}\pi^{1/2}(\frac{4pq}{p+q})^{1/2}\int_0^{|\vec{R}_P-\vec{R}_Q|}e^{-pqy^2/(p+q)}dy\\ =2M\pi^{5/2}|\vec{R}_P-\vec{R}_Q|^{-1}(pq)^{-1}(p+q)^{-1/2}\int_0^{|\vec{R}_P-\vec{R}_Q|}e^{-pqy^2/(p+q)}dy\\ \text{Let: } \sqrt{\frac{pq}{p+q}}y\rightarrow t \Rightarrow dy=\sqrt{\frac{p+q}{pq}}dt\\ \therefore \quad = 2M\pi^{5/2}|\vec{R}_P-\vec{R}_Q|^{-1}(pq)^{-1}(p+q)^{-1/2}(p+q)^{1/2}(pq)^{-1/2}\int_0^{(\frac{pq}{p+q})^{1/2}|\vec{R}_P-\vec{R}_Q|}e^{-t^2}dt\\=2M\pi^{5/2}|\vec{R}_P-\vec{R}_Q|^{-1}(pq)^{-3/2}\int_0^{(\frac{pq}{p+q})^{1/2}|\vec{R}_P-\vec{R}_Q|}e^{-y^2}dy \quad\quad(t\rightarrow y) \]

同样定义\(F_0(t)\):

\[F_0(t)=t^{-1/2}\int_0^{t^{1/2}}e^{-y^2}dy \]

误差函数为:

\[erf(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-\eta^2}d\eta \]

则有:

\[F_0(t)=\frac{1}{2}(\frac{\pi}{t})^{1/2}erf(t^{1/2}) \]

令\(t=\frac{pq}{p+q}|\vec{R}_p-\vec{R}_Q|^2\),则\(t^{1/2}=(\frac{pq}{p+q})^{1/2}|\vec{R}_p-\vec{R}_Q|\)

积分可变为:

\[= 2M\pi^{5/2}(p+q)^{-1/2}(pq)^{-1}\bigg[(\frac{p+q}{pq})^{1/2}|\vec{R}_p-\vec{R}_Q|^{-1}\bigg] \int_0^{(\frac{pq}{p+q})^{1/2}|\vec{R}_P-\vec{R}_Q|}e^{-y^2}dy\\ =2M\pi^{5/2}(p+q)^{-1/2}(pq)^{-1}t^{-1/2}\int_0^{t^{1/2}}e^{-y^2}dy\\ =2M\pi^{5/2}(p+q)^{-1/2}(pq)^{-1}F_0(t) \]

最终积分为:

\[\langle\chi_1\chi_2|\frac{1}{r_{12}}|\chi_3\chi_4\rangle = 2M\pi^{5/2}(p+q)^{-1/2}(pq)^{-1}F_0(t)\\=2\pi^{5/2}/[(\alpha+\beta)(\gamma+\delta)(\alpha+\beta+\gamma+\delta)^{1/2}]\times \exp[-\frac{\alpha\beta|\vec{R}_A-\vec{R}_B|^2}{\alpha+\beta}-\frac{\gamma\delta|\vec{R}_C-\vec{R}_D|^2}{\gamma+\delta}]F_0(t)\\ (M = \exp[-\frac{\alpha\beta|\vec{R}_A-\vec{R}_B|^2}{\alpha+\beta}-\frac{\gamma\delta|\vec{R}_C-\vec{R}_D|^2}{\gamma+\delta}])\\ \text{Here we let: }\\ t=\frac{pq}{p+q}|\vec{R}_P-\vec{R}_Q|^2=\frac{(\alpha+\beta)(\gamma+\delta)}{\alpha+\beta+\gamma+\delta}|\vec{R}_P-\vec{R}_Q|^2\\ F_0(t)=F_0(t)=\frac{1}{2}(\frac{\pi}{t})^{1/2}erf(t^{1/2})\\ \text{when }t=0: \text{let }F_0(t) = 1\\ \]

误差函数可以调库计算得到,因此积分得解



【本文地址】

公司简介

联系我们

今日新闻

    推荐新闻

    专题文章
      CopyRight 2018-2019 实验室设备网 版权所有