量子化学计算原理及其应用

陈奕驰

(西北工业大学附属中学 西安 710000)

一、量子化学的发展

量子化学是利用量子力学研究化学问题，其核心是精确求解薛定谔方程。在量子力学中，粒子的状态用波函数（概率分布函数）来描述。为了描述粒子状态变化的规律，需要找出波函数所满足的运动方程，薛定谔方程够满足上述要求。

其中第一项表示电子动能，第二项是电子互斥势，第三项是电子与原子核的吸引势，第四项是原子核之间的互斥势，第五项是原子核的动能项。

上述公式都是采用原子单位，原子单位可以使得各种计算公式中不需要涉及单位转换，一方面使得代码简洁，另一个方面也可以避免额外的计算带来的精度的损失。上述公式中电量的原子单位就是电子所带电量的绝对值。能量的原子单位是hartree 相当于2625.500 kJ/mol。

量子化学的核心问题是求解分子体系的薛定谔方程，但是由于没有办法精确求解，在实际求解中会用到各种近似，其中最常用的近似是Born-Oppenheimer 近似。[1]考虑到原子核的质量是电子质量的1780倍。原子核的运动相对于电子的运动非常慢，因此，把原子核的运动认为是静止的。此时，对于公式2中的哈密顿量。原子核的动能项为0，原子核之间的相互作用就是常数了。这样只要考虑电子的运动，原子核的位置作为参数处理，能够减少计算量，使得关注点在电子的运动。除了Born-Oppenheimer近似外，还有有限基组近似。也就是说，我们在描述一个体系的时候，需要无数的基（向量）才能准确描述，但是如果做泰勒展开的话，我们只保留展开式的前几项，后面的项都忽略，这样也能极大减少计算量。

Hartree-Fock方法，又称为自洽场迭代方法，是量子化学求解薛定谔方程的基础方法，其他的方法都是在HF方法的基础上发展而来的。我们观察公式1可以发现。目前我们只知道H的表达式，而需要求解能量E，波函数 。一个公式中两个未知数，而且是等式两边都有需要求解的波函数。自洽场迭代的思想如下。首先，对波函数 进行一个猜测，这样，公式的左边H和波函数都已知了，此时我们把波函数记为，通过薛定谔方程，我们可以求出公式右边的E和波函数。我们把此次求得的结果记为E1和，然后我们把带入到公式的左边，继续求解，得到E1和，如此循环下去，直到En和En+1之间的差足够小。这个差值是我们自己设置的。一般是要小于1kJ/mol。此时可以认为，我们求出了合适的能量E和波函数。这个过程是不断循环迭代的，所以称为自洽场迭代方法。实际计算的时候是矩阵计算，需要对矩阵进行对角化计算。

量子化学计算精度的关键在于理论方法和基组两个部分。如图所示：

图1 量子化学计算精度与理论方法、基组的关系

考虑到计算的体系大小和计算的资源等，我们需要做的是理论方法能够和基组大小达到匹配。高精度的理论方法配合高质量的基组，否则就是一种浪费。

除了Hartree-Fock计算方法外，还有一种常见的计算方法是密度泛函方法。密度泛函方法是目前实际计算中广泛使用的方法，能够在可以接受的时间范围内，给出较为精确的计算结果。密度指电子数密度；泛函是说能量是电子密度的函数，而电子密度又是空间坐标的函数；函数的函数，是为泛函（Functional）。[2]密度泛函理论是一种通过电子密度研究多电子体系电子结构的方法。具体到操作中，密度泛函理论通过各种各样的近似，把难以解决的包含电子-电子相互作用的问题简化成无相互作用的问题，再将所有误差单独放进一项中（XC Potential），之后再对这个误差进行分析。下图给出了计算精度和时间花费的关系。

图2 计算精度与计算时间的关系

二、Gaussian

Gaussian是最早的一批量子化学计算程序，由1998年诺贝尔奖获得者John A.Pople从1970年起在卡内基梅隆大学主导开发。该软件功能全面，配合GaussView图形界面。输入文件的编写是所有主流量化程序中最简单的，能够实现很多的功能。包括对于基态、激发态、反应过渡态的几何结构优化；基态、激发态能量计算，化学键能计算，电子解离能和亲合能的计算；对于红外光谱，拉曼光谱，核磁共振谱的光谱计算；偶极矩，电荷分布，电荷密度，热力学参数等的计算。

三、基组

具体计算的时候，波函数是用原子轨道线性展开的。原子轨道是由基组构成的。考虑到不同的计算要求，需要对基组添加极化函数和弥散函数。极化函数是用来描述分子轨道发生很大形变时，更准确地描述。例如，我们知道s轨道是个球形的，但是如果收到周围电荷的影响，球形可能形变成椭球形，此时，用极化函数可以描述这种变形。另外，在电子迁移的过程中，电子会弥散到整个空间，这个时候用弥散函数描述这种大的范围的轨道性质，才能更真实反应体系的性质。同时需要注意的是，基组的变大会使得计算花费的时间增大很多。

（一）实验过程

通过计算振动频率判断乙酸、乙醛、丙酮中的羰基强弱。

在GView软件中搭建乙酸分子，点击Calculate，Job Type设置为Opt+Freq，Method设置为DFT下的B3LYP法和3-21G基组，打开Gaussian03W进行计算，再用GView打开输出文件，在Result一栏点击Vibrations，在弹出的对话框中设置Show Displacement Vectors，点击各个振动频率，观察选中哪一个频率时箭头显示为羰基的振动频率，将其记录下来。完成后对乙醛、丙酮进行相同操作。以上全部完成后，将基组改为6-31G，方法不变，重复全部操作。

（二）实验结果

图3 B3LYP/6-31G 优化后的（a）乙酸；（b）乙醛；（c）丙酮

表1 B3LYP/6-21g//6-31g乙酸、乙醛、丙酮羰基的振动伸缩频率的计算结果

羰基由强到弱排列顺序：乙酸、乙醛、丙酮。

四、结果讨论

由于化学键振动伸缩频率越高强度越大，我们从计算结果可以得知在使用B3LYP方法的前提下，使用3-21G基组计算时得出的结论是羰基强度从强到弱顺序为乙酸、丙酮、乙醛，使用6-31G基组时从强到弱顺序则为乙酸、乙醛、丙酮。要试图解释此结果可以考虑从羰基氧和羰基碳原子的电荷分布入手，而这正是从输出文件中可以得到的，具体数据见表2和表3。

表2 B3lyp/3-21G计算后的乙酸、乙醛和丙酮中羰基碳和羰基氧的电荷分布

由表2可见，用3-21G基组算出的羰基碳氧电荷差值的大小顺序与羰基的相对强弱顺序是吻合的，从静电作用角度理解，碳原子和氧原子间更大的电荷差距使得正负电荷之间的库伦相互作用增强，对于键能的增益更加明显。

表3 B3LYP/6-31G 计算后的乙酸、乙醛和丙酮中羰基碳和羰基氧的电荷分布

由表3可知，从羰基碳氧电荷差值大小顺序来看，6-31G基组算得的结果与3-21G基组相同，皆为乙酸＞丙酮＞乙醛，但是从振动伸缩频率上看，6-31G基组算得的羰基强弱顺序为乙酸＞乙醛＞丙酮，与3-21G的结果以及6-31G的电荷差值结果不符。因为尽管3-21G与6-31G属于不同的基组，但是它们都属于劈裂价键基组，对这两个基组而言，每个价层电子轨道都会被劈裂为两个基函数，3-21G基组中价层电子轨道由3个和1个高斯型函数线性组合而成，每个内层电子轨道由3个高斯型函数线性组合而成，而6-31G基组中后者为6个，价层电子轨道由3个和1个高斯型函数线性组合而成，每个内层电子轨道由6个高斯型函数线性组合而成，二者的结果为何不同还需要进一步探索。

五、结束语

计算化学的发展帮助我们更好地理解实验现象并且给出定量的解释。但是同时，我们要意识到很多时候，计算的结果和实验结果有出入的时候，要明确我们计算的物理量是否就是实验上的。比如，我们研究是羰基的键的活泼性，而不是羰基的键和键共同的。我们计算的频率是双键的频率。因此，要做出区分。Gaussian软件结合Gaussview软件的作图能帮助快速计算，并给出定量的结果去讨论，而计算结果需要结合有机化学的知识去进一步分析。