关于HF、H2O的电荷计算

Zhou2009

在寻找与电负性相关性好的σ诱导电荷时， Mulliken、NBO等电荷都不合用。必须找到它们为什么不合用，到底问题在何处？这样才能试着找到合理、合用的电荷数值。在对分子中电荷计算的深入探讨过程中，HF、H2O是首先要观察、计算的有代表性典型分子，也许说例子比全盘分析还清楚一些，今后也只以例子分组贴出。显然这些还只是计算、思考的笔记，只是前进的足迹。

在多次改稿的过程中Sobereva都参加了讨论并提出了宝贵的意见，特此致谢！

1．        从Mulliken、NBO电荷与电负性相关性不好谈起：
前帖“原子电负性与其诱导电荷”，用原子电负性及其诱导电荷（用局部密度差HDD求得）曾经得到图1。
而用另外原子电负性及其诱导电荷（分别由Mulliken atomic charges和NBO atomic charges求得）曾经得到图2。




从图1看，原子的σ诱导电荷的值与电负性有良好的相关性。原子的σ诱导电荷的值，是随原子的电负性增大而增大的。这种横向的比较表明，在各分子体系中，量化计算的电荷值是统一的、相互之间是可以定量比较的。

从图2看，原子的σ诱导电荷Mulliken charges和NBO charges与电负性的相关性出现了严重的问题。而且Mulliken charges和NBO charges在图形上何其一致，表明它们有着共同的误差根源，NBO也并没有根本性的改进。
我们平日只是从实践经验中（如电负性）个别地感到Mulliken、NBO的电荷尚有问题，然而图2这两条线却从横向展示了问题的严重性，但这也揭示了产生这些问题本源的线索。从图2看，C、N、O的诱导电荷竟比F还要大，Cl又一同掉了下来！从图形看肯定是有问题的，它是怎样导致的？

在探讨这方面的问题时，入手分析Mulliken charges都采用HF/STO-3G计算，这样指认起来比较简单明了。进一步取能量、电荷数据时，则会采用较高基组HF/6-311++G**，使数据更准确一些，到时会作说明。

电荷的本源无非是MO-LCAO，由此就有了原子上的布居电荷和重叠布居电荷，利用得到的wfn文件，进入Multiwfn对分子的MO和MO片段进行计算，并与Mulliken population analysis数据作对照，可以深入了解Mulliken population的构成和来源，特别是要找到什么地方引起了差错导致电荷计算的失误。
同时用HDD方法已经得到的数据数量级来作参考依据，因为HDD得到的σ诱导电荷的值毕竟是与电负性有良好相关性的（图1）。

本来经过量化计算得到的系列MO，进一步得到的电子密度分布，只是在分子的空间取格点，求每一个格点的电子密度ρ。这里并没有原子上密度、也没有难分难解的重叠密度问题，电子密度ρ就是一个空间分布，它可以是整个分子的，也可以是某一个或某几个MO的，还可以是一个MO的部分片段。
但是为了实际应用，我们还是需要得到分子中原子上电荷的实际数量，比如σ键诱导电荷对相关原子来说，有多少电荷发生了转移呢？而且这种电荷转移又确实是存在于空间分布的电子密度ρ之中的。

怎样在量化计算的实际结果中确定我们所需要的电荷数值呢？
具体说在一个分子H-R中（-R为原子或基团），怎样从它的电子密度ρ，给H划定一个空间去得到H上的ρ或者ρ的变化呢？这样就可以进一步得到具体的电荷数值。
人为地划分空间并不容易，也不那么合用，历来许多学者都作过这样的探索了。

这时不免会想到：如果就拿一个孤立H原子的电子密度ρ和孤立R的电子密度ρ（或者是R中各原子的电子密度ρ），分别在H-R中H的位置和R的位置，与H-R的电子密度ρ来作电子密度差Δρ，不就可以比较出在成键前后这个H和R空间上的ρ的净变化了吗？用图形展示这种Δρ，研究电子密度ρ在空间的变化，这已是国际量子化学研究的通用的方法了。我们还试探了对Δρ或增或减的各部分作进一步定量。
但在实际操作时，有时R不易计算确定时，特别是R中各原子的电子密度ρ不易计算确定时，又想到：如果就拿一个孤立H原子的电子密度ρ，在H-R中H的位置，与H-R的电子密度ρ来作部分电子密度差Δρ，即HDD方法，不就可以比较出在成键前后仅H这个空间上的ρ的净变化了吗？

无论什么情况，现在只计较H这个空间的实际具有的电荷变化总应是合理、方便、合用的。
当R的电负性大于H时，H失去电荷，HDD的这个H处的电荷差值就为负值，HDD为三维空间格点文件cube，HDD进入GsGrid程序，就会计算出H处负值电荷的具体数值，这个值主要就是R从H上诱导走的电荷。前面图1就是这样得到的电荷，它与作为R的原子电负性有较好的相关性。其实，电荷与电负性的这种相关性，已经开始成为寻找合用电荷数量级的判据。用这种方法，可以从H上诱导走的电荷反过来推知-R的电负性数值，无论这-R是原子或基团。这时H就好像一个电负性的测量器，从H上被诱导走的电荷数值（HDD）可以知道此时分子中H实际具有的电负性，并用来测量-R的电负性数值。
当然，本帖虽然用HDD作为从H上被诱导走的电荷数值的数量级参照，但基本计算还是根据从MO得到的原子上的布居电荷和重叠布居电荷来进行分析。

当然，此时要排除非σ键诱导电荷对H电荷的叠加影响，如孤对电子反馈、超共轭等。

2．        HF 分子Mulliken charges的解析和深入探讨：
G03对HF分子作的HF/STO-3G计算输出Full Mulliken population analysis：
                             1         2         3          4            5           6
1 1   F  1S      2.11665
2        2S     -0.11740   2.13199
3        2PX    0.00000   0.00000   1.00200
4        2PY    0.00000   0.00000   0.00000   2.00000
5        2PZ    0.00000   0.00000   0.00000   0.00000   2.00000
6 2   H  1S     -0.00022  -0.07034   0.26556   0.00000   0.00000   0.59415
由此可以得到Mulliken atomic charges（净电荷）：
1  F    -0.210846      
2  H     0.210846

2.1  Mulliken population analysis中对角元H上的电荷计算：
在上面的population analysis中，对角元为相应原子上基函数所形成的电荷，非对角元为相应两个原子基函数重叠布居电荷。
运行multiwfn读进HF.wfn，用功能-3，选1（片段只包含一个原子），然后输入2（片段包含2号原子，即H），然后计算电子密度格点文件，得到得到H上的电荷0.5937。
可以从上面的population analysis矩阵对角元H看到相应的值：0.5942。它们是吻合的。
现在更深入地看看这个Mulliken矩阵对角元H电荷0.5942的构成，即它包含些什么？怎样得来的？

从计算结果看，HF分子共有5个MO，MO1为F的1s、MO4、5为F的孤对电子轨道，它们与H的诱导电荷都没有关系，其中只有2σ（MO3）是名副其实的分子轨道，而MO2也包含H基函数，也必须考查一下。

还是老办法，利用得到HF.wfn文件。先只看MO3这个成键轨道。
运行multiwfn读进HF.wfn，用功能-1，选择MO3，再用功能-3，选1，然后输入2（即H），计算电子密度格点文件，得到得到H上的电荷为0.5435。
运行multiwfn读进HF.wfn，用功能-1，选择MO2，再用功能-3，选1，然后输入2（即H），计算电子密度格点文件，得到得到H上的电荷为0.05412。
0.5435+0.05412=0.5936
这样就得到了Mulliken矩阵对角元H电荷0.5942。于是我们就知道Mulliken这个矩阵对角元H电荷数值的来源、构成了，它是MO3、MO2中H上电荷的加和。

2.2  Mulliken重叠布居电荷的计算：
同上面一样，运行multiwfn读进HF.wfn，用功能-3，选1，然后输入1（片段包含1号原子，即F），最后计算电子密度格点文件，得到电荷8.9865。这样再计算H的，得到电荷0.5937。算HF整体得到9.9703。于是9.9703-0.5937-8.9865=0.3901，这就是H、F基函数叠加形成的重叠布居电荷。
这也可从输出的population analysis矩阵计算到同样的值：F的2s与H的重叠布居电荷为-0.1407（矩阵列出的值只是重叠布居电荷的一半）。F的2p与H的1s的重叠布居电荷为0.5311。这两个重叠布居电荷之和为0.3904。

现在更深入地看看这个Mulliken重叠布居电荷0.3904的构成，即它包含些什么？怎样得来的？
还是老办法，利用得到HF.wfn文件。先只看MO3这个成键轨道。
前面已经作过，运行multiwfn读进HF.wfn，用功能-1，选择MO3，再用功能-3，选1，然后输入2（即H），计算电子密度格点文件，得到得到H的电荷为0.5435。照此办理，求得F的电荷为1.3426。
而对MO3全部基函数求得的电荷为2.0026。它与分别单独求得的F的电荷加上H的1s的电荷相差0.1165，这个新增的差值就是MO3中F与H重叠布居电荷。

再算一下MO2。
运行multiwfn读进HF.wfn，用功能-1，选择MO2，再用功能-3，选1，然后输入2（即H），计算电子密度格点文件，得到得到H的电荷为0.0501。照此办理，求得F的电荷为1.6702。再将HF.wfn读入Multiwfn，求得MO2总的电荷为1.9985。。
1.9985-1.6702-0.0501=0.2782
此即是MO2中F与H之间形成的重叠布居电荷。将MO3、MO2中F与H之间形成的重叠布居电荷相加，
0.1165+0.2782=0.3947
这样也就得到了Mulliken重叠布居电荷之和0.3904。于是我们就知道Mulliken这个重叠布居电荷数值的来源、构成了，它是MO3、MO2中H与F上重叠布居电荷的加和。
Mulliken将H上布居电荷加上重叠布居电荷的一半，就得到H原子的电荷了：
0.5942+0.3904/2=0.7894，
并由此得到H的Mulliken atomic charges（净电荷）：0.2106。

然而，化学实践的经验告诉我们，这个H的净电荷值太小了！

2.3  MO2是什么？
其实MO2应该只是F的2s轨道。它的能级为-1.60382（用6-311++G**计算，孤立F原子的2s为-1.67382，HF成键的MO3能级为-0.77001）。
MO2中的基函数，只是F的2s轨道在H的摄动下极化向空间展开它的形象、或说象黄油那样抹开成电子云时的轨迹，这时的基函数不是指真实的AO。即此时F2s轨道还是它自己，一个完整的2s轨道，基函数中的Hs并不是一组真实的AO，只是F的2s轨道的一个踩点而已。

此时MO2中Fsp与Hs的关系完全不同于MO3中的。
在MO3中，Hs既是一组基函数，它同时又是真实的AO，Hs有电荷转移到Fsp上了，它们之间还形成了重叠布居电荷0.1165。（此时，MO3中的2s也只是基函数而已，也不是一个真实的AO，它体现了sp杂化。）
可见在MO2中，这时所谓Hs电荷0.05412和基函数Fsp与Hs的重叠布居电荷0.2782这些都还只是F2s轨道的组成部分。当我们对整个MO2的电子进行加和时，它严格的是2！

我们还可以作出MO3的电子密度图（图3）和MO2电子密度图（图4）来作比较。
MO3除了F核为中心的2sp外，在以H为中心处还可见电子云涡旋，表明此处存在以H核为中心电子聚集。


MO2（图4）表明整个图形是以F核为中心的2s轨道，但是延伸到H处受到极化变形。此时基函数Hs安排了电荷0.05412，这在Mulliken电荷中也是要记入H的。


F2s这种极化、延伸，是由于H在与F成键时，H向F转移了较多电荷，在H与F之间还提供了所谓重叠布居电荷，使H核外电荷减少，吸引电子的能力增强所引起的。与此同时，F核外电荷增多，吸引电子的能力减弱，但由于2s毕竟是内层轨道，受影响较小，能级只是由孤立F原子的2s的-1.67382上升至-1.60382，并且因H的影响表现为受静电引力而极化变形。

其实F核外电荷增多使得其吸引电子的能力降低，对其上的孤对电子更是有直接影响的。孤立F原子的孤对电子为二重简并，能级为-0.73347，而在HF分子中孤对电子仍为二重简并，能级升高至-0.65110。这孤对电子也会受H核影响的，对比HF分子中孤对电子与孤立F原子的孤对电子的电子云分布，可以发现前者的电荷分布因静电引力向H也有前倾趋势，然而由于对称性不同、正交，不能形成所谓孤对电子反馈（此图易想象未列出）。电子云是多么灵动啊！

至此，Mulliken矩阵对角元H的电荷和重叠布居电荷的构成就找到了，它把MO2中，H上的电荷和重叠布居电荷也都等同于MO3中的了，并将它们加和起来作为H上的电荷和重叠布居电荷，这样H上的电荷就由0.5435算成了0.5937，重叠布居电荷就由0.1165算成了0.3904，这里存在着严重多算！这种均分MO2中过大的重叠布居电荷一个直接的影响是它缩小了原本的电荷因转移形成的差距。这样就使得H的电荷计算过多，使H的净电荷过于偏小！（注意这里的“电荷”是指H上实际所具有的“电荷值”，而“净电荷”则是核电荷（1）减去这个“电荷值”，表明H上电荷的净变化。）

这种将2s也纳入的严重多算，就是Mulliken布居分析的计算失误之一！

这不仅使用重叠布居来表征的Mulliken键级失误，把这原子之间的重叠布居的电荷均分给相应的两个原子，又造成Mulliken原子净电荷的失误。

2.4 问题就出在Mulliken布居分析的基本思想：
Mulliken布居分析的基本思想是：因为体系的整个电子都是用各个原子基函数表示的，只要把这些原子基函数再分给这些原子的轨道就可以了。

如今基函数已经不能完全等同于原子的轨道了，基函数常常是一个MO展开的态势，如在HF中MO2中基函数F2s、F2p与H1s，不能将HF中MO2中的电荷真正分归到F2s、F2p与H1s，也不能把它们之间的重叠布居电荷看作是H与F之间的，因为它只是F2s自己而已。

我们在考察分子中σ键诱导电荷时，曾经在CH2O中，需要将孤对电子反馈在H上的电荷叠加挑剔出去不于考虑。现在在HF中，我们在考虑H的电荷时，又需要将F的2s轨道挑剔出去不于考虑。因为他们都不直接是σ键诱导电荷的内容。

2.5  对σ键诱导电荷的计算：
在Molecular Orbital Coefficients矩阵中，我们不仅从纵向可以看到每个MO都是归一的，而且从横向看每个基函数的电荷又是可以加和，它又确实表征着该基函数与它联系的原子轨道所占据的电荷数值的。该矩阵是一个归一化的严格自洽系统，是Mulliken布居分析的出发点。

由于考虑到诱导电荷只是σ键的内容，我们试着只从MO3来推算诱导电荷的数值。为了使数值更精确一点，转而采用HF/6-311++G**计算数据。

从上面算例可以看到，无论是在Mulliken布居中还是用Multiwfn程序来计算电荷，原子电荷（对角元）因直接从MO系数来，它们是一致的。现在重点是考察MO层次的重叠布居电荷。

拿这里的HF分子来说，重叠布居电荷是HF某一个MO的电子密度ρ与组成这个MO的H的ρ和F的ρ之差，它是新的一种Δρ。通常这种相减原子部分的ρ抵消了，只剩下新形成的ρ（由此可计算的重叠布居电荷）。如图5就是HF分子MO2的这种Δρ，重叠布居分布在二原子之间。以至于可以用MO的ρ的数值，减去组成这个MO各片段的ρ的数值，就可得进一步计算得到重叠布居电荷。


在Mulliken布居分析中，或者用Multiwfn程序来计算，对HF分子的 MO3的重叠布居电荷也就是这样计算的。
这样求得MO3的电荷值为1.9995，H的电荷值为0.3534，F的电荷值为1.4509，重叠布居电荷为：1.9995-0.3534-1.4509=0.1952
这重叠布居电荷的一半0.0976划归H，加上H原有电荷0.3534，共有电荷0.4510，净电荷为+0.5490，与HDD得到的净电荷0.3990相比较，感到净电荷过于大了！

原来，MO3的所谓重叠布居电荷不能简单地这样用上面的Mulliken布居分析电荷数值差来计算。这可以作出MO3的Δρ（图6）来进行分析，它表达了MO与组成它的原子具有了组合系数的基函数电子密度的净变化。图6所展示的信息，与图5完全不一样！图6表明，MO3中独立的F 基函数片段，是sp型杂化，并且是以杂化后的小瓣朝向H，但一旦与H的基函数组合之后，却发生F的sp向H转移电荷的现象，转移电荷为0.2564，要完全算在H上，而不是象图5那样分布在F、H之间需要均分。当然这并不表明只是F向H转移了电荷，而是MO3中独立的F基函数片段分配了过多的电荷，一旦与H组合，电荷又进一步作了调整而已。


转移电荷0.2564，加上H上原有的电荷0.3534，共有电荷0.6098，净电荷为+0.3902，与HDD得到的净电荷0.3990相比较，是符合的！
而Mulliken布居分析将这里的所谓重叠布居电荷的一半0.0976划归H，这个数值与转移电荷0.2564相比较，差得太远。这是Mulliken布居分析的计算失误之二！

Mulliken布居分析一方面将HF分子的MO2（F2s）这种极化、延伸所形成的所谓重叠布居电荷的一半划归了H（0.1494，由Δρ求得），还给出了H上的电荷（0.05012），使得H的电荷算得过多，多算了0.1995；另一方面却又将MO3中的电荷调整转移（0.2564）只以所谓重叠布居电荷的一半（0.0976）计入H，电荷算得又过少了，少算了0.1588。而且这一个“过多”、一个“过少”并没有什么相关性而作为系统误差相抵消，相反二者各有千秋。多算了0.1995，少算了0.1588，总地来说还是多算了0.0407，这使得HF分子的电荷转移量（即净电荷）算得过于小了。

上面这种问题是普遍存在的，它造成了Mulliken布居分析给出的电荷值存在很大的失误。

现在电子密度差Δρ，已经由分子片作到MO的基函数片的层级了。当然它与用孤立原子分子片作的电子密度差Δρ，在概念、数值意义上又是有显著差别的。
MO的基函数片的Δρ，集中表达了各具有了组合系数的原子基函数与它们组成MO之后电子密度的净变化，定量地给出了重叠布居电荷和新产生的电荷在原子之间的调整。这不是“把这些原子基函数再分给这些原子的轨道”，只对原子的轨道进行加减计算就可以解决的问题，因为MO的基函数线性组合产生了许多新的质！而且一个MO的Δρ不一定总是在展示重叠布居的信息，需要具体分析Δρ的图象来分析判断。

    这里还只是双原子分子，在多键原子的分子中，问题又是怎样存在着呢？且算一个有代表性的H2O的例子考察、检验一下。

[ Last edited by yjcmwgk on 2010-7-16 at 11:47 ] 返回小木虫查看更多