简单利用Mathematica计算氦合质子的分子轨道
—Edward Calhoun
1. 引言
量子化学是理论化学的一个重要分支学科，自20世纪60年代计算机技术突飞猛进，计算化学也同样成为量子化学研究的有力交叉学科，基于非相对论近似，绝热近似和单电子近似的全电子的非经验计算方法被迅速提出，其包括不借助任何经验参数求解全电子体系的Schrödinger方程方法，自洽场水平上的超Hartree–Fock-Roothaan方程计算法（HFR法）。原则上讲，有了HFR方程，就可以计算任何多原子体系的电子结构。
1928年，D．R．Hartree提出了解决分子和离子波函数和能量问题的假设，Hartree将n个电子体系中的任意一个电子都看成是其余电子所提供的平均势场中运动的单电子，即可导出可表示任意一个电子的单电子方程，再用自洽场法进行迭代，所得方程成为Hartree方程，但由于Hartree并未考虑到Pauli原理，2年后，由B.A.Fock和J.C.Slater进行了Pauli修正，加入了单Slater行列式波函数和多Slater行列式波函数，使得HF方程正式被提出，而C.C.J.Roothann提出的LCAO-MO法进一步简化了HF方程的表达形式，利用LCAO法可以将分子轨道看成是原子轨道的线性组合，即将HF方程转变为易于求解的代数方程组，使得HFR方程（1-1式）得以成型。其中
$$f_{(1)}$$为Fock算子， $$\chi_{a(1)}$$为体系的分子轨道。
（1-1式）
$$f_{(1)}\chi_{a(1)}=\in_a\chi_{a(1)}$$
20世纪70年代，基于HF方法的密度泛函理论（DFT）被提出，尽管DFT的概念起源于Thomas-Fermi模型，但仍然离不开HF方程的理论基础，DFT最大的优点在于利用电子密度替代波函数作为研究的基本量，因为以每个电子作为质点，N个电子有3N个变量，而用整体密度代替，则缩减到3个变量，大大减少了计算量。
氦合质子，化学式HeH+，是热力学理论计算最强的Bronsted酸，其pKa≈-63，HeH+无法稳定存在于凝聚相中，其pKa数据根据《无机化学》上Hess定律设计Born-Harbor循环可以计算得出。HeH+与同核分子H2有所不同，HeH+有永久的偶极矩，使得其有强烈的光谱特征，这也是HeH+发现的来源，通过观察光谱特征判断该物种的存在与否。
本文尝试利用Mathematica软件可视化HeH+的分子轨道，利用HFR法首先要进行电子积分的求解，Mathematica无疑是强有力的工具。
2. 正文
分子积分包括重叠积分，动能积分，势能积分和双电子积分，这个部分是量子化学计算量最大的地方，对于本弱的高数基础，不能完全胜任，需要借助Mathematica语言解决问题。由于其计算量和计算基组有关，这里本弱用经典的STO-3G基组进行分子积分，首先输入两个基组 $$\varphi_1$$$$\varphi_2$$
为计算方便，确定He原子坐标为原点（0,0,0），H原子坐标为（1.4632,0,0）。计算重叠积分$$S_{12}$$ ，输入语句如下：Integrate[φ1*φ2,{x,-∞,+∞},{y,-∞,+∞},{z,-∞,+∞}]
计算动能积分$$T_{12}$$Integrate[(-0.5)*φ1*(D[φ2,{x,2}]+ D[φ2,{y,2}]+ D[φ2,{z,2}],{x,-∞,+∞},{y,-∞,+∞},{z,-∞,+∞}]
对于势能积分-NIntegrate[φ1*2/Sqrt[x^2+y^2+z^2]* φ2,{x,-∞,+∞},{y,-∞,+∞},{z,-∞,+∞}]
对于双电子积分，经过Laplace变换后可得解析式，但Mathematica不支持该解析式，只能寻求以下做法。
HeH+的双电子积分是4个Gauss函数的双电子积分式，根据LCAO，一个基函数由3个Gauss函数构成，基本组合原理指出，基函数的双电子积分由3^4=81组Gauss函数的双电子积分加和而得，那么首先在Mathematica内定义Gauss函数f[],进行一次双电子积分，再用Do命令循环81次再予以加和，过程本弱不贴示（←太懒）
有了分子积分的结果就好说了，接下来只需要分别写成矩阵形式，按照HFR方程求解，很轻易的就可以解出HeH+的分子轨道了。
HFR方程的求解过程其实就是迭代过程，方程反复迭代，直到体系能量和密度矩阵达到自洽的结果为止。Mathematica在这段计算中起到了方便的功效，HFR的迭代需要的数学运算有矩阵的加法乘法、转置和特征值的计算，徒手解矩阵那是dalao做的事情，本弱的线代基础只够借助计算机算。
不废话，根据基本数学性质，Fock算子$$f_{(1)}$$F={[-1.50621,-0.363043],[-0.363043,-0.152917]} Eigenvalues[F] Out[5]={-1.59745,-0.0616766} 这两个值就是HeH+的两个分子轨道的能量。
下面是关键的可视化环节。
根据LCAO，结合计算得到的系数矩阵，可以得到分子轨道数学表达式。 Mathematica 软件提供了丰富的作图命令，可以绘制轨道的径向图、网格图和等值图等。
Fig-1就是HeH+的ψ1-x占据轨道径向图，和经典教材上的分子轨道图像一致，首先是Bonding Orbit上两原子间电子云密度增大，符合基本分子轨道原理，其次He的电子云密度比H的要高，在图像上也很容易体现出来。
3. 后记总结
Mathematica在进行HeH+的量子化学计算时，极少的运用编程知识，是对于本语言渣的一个巨大福音。
了解量子化学，可视化的图像会比枯燥的行列式或者微分方程更好接受，但手动绘图的过程及其艰辛，尤其是要手解矩阵，手动Laplace变换的环节，这些是我这种弱渣很难做到的，只有julao才能徒手解量子力学的方程。
HeH+作为一个竞赛模拟题里出现的物质引起了我很大的兴趣，于是有感而发运用浅陋的Mathematica语言写下了这篇水文，
其作为最强的Bronsted酸，比高中化学常见的硫酸要强上60多个数量级，这些是难以想象的。He的化合物还有很多，例如Na2He，在Nature Chemistry期刊上发表的新物种，31届的国初题正好考察到了它的相关晶体结构，尽管来讲该物种可能只归类于电子盐，但是也是化学界的一个重大发现，然而本弱依然在此题上翻了车，这是重大的遗憾。
不扯太多，量子化学计算一般有化学软件代替完成，其过程犹如黑箱，是不可知的，市面上有许多从头计算法和运用DFT计算的软件，虽然应用方便，但是对于学习该门课程来讲是于事无补的，然而手工计算量子化学几乎是登天的事情，这时就需要数学软件的辅助，才能透彻的了解经典模型，才能真真正正的把《结构化学基础》《高等无机结构化学》的量子化学部分弄明白。可以这样说，踏入量子世界，要么脑壳里要装着Bohr，Heisenberg的大脑，要么兜里要兜着一台计算机，不然就只能乖乖离开。