谈谈NBO中用bndidx关键词输出的三种键级

注：一般没必要用NBO算键级。Multiwfn（可在//www.umsyar.com/multiwfn免费下载）支持非常全面的键级计算方法。怎么用Multiwfn算各种键级看Multiwfn手册4.9节的例子，各种键级原理介绍见Multiwfn手册3.11节，对于Gaussian用户只要提供fch文件即可，敲两下键盘键级就输出出来了。更多和键级有关的博文看《谈谈原子间是否成键的判断问题》（//www.umsyar.com/414）、《通过键级曲线和ELF/LOL/RDG等值面动画研究化学反应过程》（//www.umsyar.com/200）。

On the three bond orders outputted by bndidx keyword in NBO

文/Sobereva@北京科音2010-Nov-19

1. 前言

用了bndidx关键词后，NBO程序会进行键级分析，输出三种键级矩阵，包括：
1 Wiberg bond index matrix in the NAO basis
2 Atom-atom overlap-weighted NAO bond order
3 MO bond order
矩阵中(i,j)元素对应i、j原子间的键级。对这三种不同键级，有些人对应当选取哪个感到疑惑，此文将谈谈这三种键级的含义及其合理性。值得一提的是，经常有人说“NBO键级”，这种说法是很不严格的，在NBO程序中，上面已经看到有三种键级了，NBO还能产生LMO键级、NRT键级，另外，NBO程序搜索出的NBO们的占据数本身也可以用于衡量键级。所以，不应当使用含糊不清的“NBO键级”这个词。

为了简明且便于验证，本文例子使用HF/STO-3G的基态水分子波函数。下面的是用DMNAO关键词输出的密度矩阵，将在后文用到。后文直接将此矩阵称为DMNAO。
NAO 1 2 3 4 5 6 7
---------- ------- ------- ------- ------- ------- ------- -------
1. O 1 (S) 2.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
2. O 1 (S) 0.0000 1.7325 0.0000 0.0000 0.3538 0.4113 0.4113
3. O 1 (px) 0.0000 0.0000 2.0000 0.0000 0.0000 0.0000 0.0000
4. O 1 (py) 0.0000 0.0000 0.0000 1.1399 0.0000 0.7001 -0.7001
5. O 1 (pz) 0.0000 0.3538 0.0000 0.0000 1.5321 -0.5439 -0.5439
6. H 2 (S) 0.0000 0.4113 0.0000 0.7001 -0.5439 0.7977 -0.0623
7. H 3 (S) 0.0000 0.4113 0.0000 -0.7001 -0.5439 -0.0623 0.7977

2. 第一种键级(Wiberg bond index matrix in the NAO basis)

说起第一种键级，先要回顾一下Wiberg键级和Mayer键级的原始定义。Wiberg键级最早来自Tetrahedron 24,1083，键级W_AB=∑[i∈A]∑[j∈B]P(i,j)^2，其中P为密度矩阵，这个定义只适合正交基。Mayer在CPL 97,270里将其扩展到非正交基，A、B间的键级定义为∑[i∈A]∑[j∈B]PS(i,j)*PS(j,i)，这里PS(i,j)是指P和重叠矩阵S做矩阵乘所得矩阵的(i,j)元素。当使用正交基后，变为∑[i∈A]∑[j∈B]P(i,j)*P(j,i)，由于P是对称的，所以又还原到了Wiberg键级。Mayer键级数值比较符合一般化学观念，比如单重键约为1，三重键约为3，苯的C-C键约为1.5。

NBO的各种分析主要基于由原始非正交基函数转化成的自然原子轨道(NAO)，NAO是正交基，所以可以直接使用Wiberg键级进行分析，也就是NBO输出的Wiberg bond index matrix in the NAO basis。这是三种键级中唯一推荐的键级。NBO给出的结果：
Atom 1 2 3
---- ------ ------ ------
1. O 0.0000 0.9552 0.9552
2. H 0.9552 0.0000 0.0039
3. H 0.9552 0.0039 0.0000
可以很容易地通过NAO的密度矩阵DMNAO来计算此矩阵。此矩阵的(2,3)代表H-H的键级，由于6、7号NAO分别是两个氢上的唯一基函数，键级也就是DMNAO(6,7)^2=(-0.0623)^2=0.00388。再比如(1,3)代表的O-H键级，由于1~5号基函数对应O，也就是∑[i=1,2,3,4,5]DMNAO(i,7)^2=0.95513。

在NBO输出这个矩阵之后还会输出Wiberg bond index, Totals by atom，指的就是每个原子形成键级的总和，也称为原子价。对于第二种、第三种键级也是如此。

3. 第二种键级(Atom-atom overlap-weighted NAO bond order)

这种键级号称是考虑了原子轨道间重叠，引入距离因素，并使得结果能看出成键反键，但实际结果不好，而且其定义我认为没有任何根据。A、B间键级定义为∑[i∈A]∑[j∈B]P(i,j)*S'(i,j)，其中S'是PNAO重叠矩阵（注：PNAO是原始基函数向NAO转换的中间基函数，PNAO在原子内正交而在原子间未必正交）。这种定义很莫名其妙，密度矩阵用NAO的，但重叠矩阵却用PNAO的。从其形式上，可以想起Mulliken键级。Mulliken键级就是Mulliken布居分析得到的两原子间交叉项的总和，∑[imo∈占据]∑[i∈A]∑[j∈B]2*占据数(imo)*C(i,imo)*C(j,imo)*S(i,j)，写成密度矩阵的形式就是∑[i∈A]∑[j∈B]2*P(i,j)*S(i,j)。所以，第二种键级实际上是一种既不完全是NAO、也不完全是PNAO下的不伦不类的Mulliken键级，而且还少乘了2。这种古怪的键级缺乏依据，实际结果又差，所以不要使用。NBO输出的结果为：
Atom 1 2 3
---- ------ ------ ------
1. O 0.0000 0.5399 0.5399
2. H 0.5399 0.0000 -0.0145
3. H 0.5399 -0.0145 0.0000

这也容易验证，使用SPNAO关键词可以输出PNAO重叠矩阵S'
NAO 1 2 3 4 5 6 7
---------- ------- ------- ------- ------- ------- ------- -------
1. O 1 (S) 1.0000 0.0000 0.0000 0.0000 0.0000 0.0935 0.0935
2. O 1 (S) 0.0000 1.0000 0.0000 0.0000 0.0000 0.4675 0.4675
3. O 1 (px) 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 0.0000
4. O 1 (py) 0.0000 0.0000 0.0000 1.0000 0.0000 0.3205 -0.3205
5. O 1 (pz) 0.0000 0.0000 0.0000 0.0000 1.0000 -0.2265 -0.2265
6. H 2 (S) 0.0935 0.4675 0.0000 0.3205 -0.2265 1.0000 0.2328
7. H 3 (S) 0.0935 0.4675 0.0000 -0.3205 -0.2265 0.2328 1.0000
设T(i,j)=DMNAO(i,j)*S'(i,j)，则T为：
1 2 3 4 5 6 7
1 2.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
2 0.00000 1.73250 0.00000 0.00000 0.00000 0.19228 0.19228
3 0.00000 0.00000 2.00000 0.00000 0.00000 0.00000 0.00000
4 0.00000 0.00000 0.00000 1.13990 0.00000 0.22438 0.22438
5 0.00000 0.00000 0.00000 0.00000 1.53210 0.12319 0.12319
6 0.00000 0.19228 0.00000 0.22438 0.12319 0.79770 -0.01450
7 0.00000 0.19228 0.00000 0.22438 0.12319 -0.01450 0.79770
求O-H键的这种键级也就是∑[i=1,2,3,4,5]T(i,7)=0+0.19228+0+0.22438+0.12319=0.53985，和NBO输出的一致。

另外值得一提的是，所谓NPA电荷（无数人谬称为NBO电荷），就是以NAO为基获得的Mulliken电荷。A原子布居数=∑[i∈A]∑[j]P(i,j)S(i,j)=∑[i∈A]∑[j]P(i,j)δ(i,j)=∑[i∈A]P(i,i)，所以氧的布居数∑[i=1,2,3,4,5]DMNAO(i,i)=2.0000+1.7325+2.0000+1.1399+1.5321=8.4045，净电荷为8-8.4045=-0.4045，和NBO程序输出的一致。

4. 第三种键级(MO bond order)

这个键级的名字十分令人疑惑，而且手册上也没说明。实际上它的定义极其简单，A、B间的键级为∑[i∈A]∑[j∈B]P(i,j)。NBO输出为
Atom 1 2 3
---- ------ ------ ------
1. O 0.0000 0.5675 -0.8328
2. H 0.5675 0.0000 -0.0623
3. H -0.8328 -0.0623 0.0000
求O-H2键级就是∑[i=1,2,3,4,5]DMNAO(i,7)=0.4113-0.7001-0.5439=-0.8327。结果十分荒诞，明明是典型sigma键，键级却为负。而另一个等价的O-H3键键级却是正值0.4113+0.7001-0.5439=0.5675！可见，这种键级从实际应用上看毫无意义。从理论上同样没有根据，从形式上看和Mulliken键级相仿佛，却少乘了重叠矩阵元，另外NAO又是正交基，不存在交叉项，也不可能求Mulliken键级。而与Wiberg键级相比，又少求了平方。我认为NBO程序中应该去掉这种毫无意义的键级。

5. 用PNAO计算Mayer键级的尝试

原理上，以非正交基PNAO做Mayer键级分析也是可以的，可以代替正交基NAO下的Wiberg键级分析，另外，做Mulliken分析通过交叉项还可以研究成键反键和各分子轨道贡献，但是NBO程序里没有提供。笔者尝试利用PNAO做Mayer键级分析，这首先需要获得PNAO密度矩阵，NBO3.1里不支持输出它的DMPNAO关键词，故通过.47文件在GENNBO5.0下输出，所得矩阵如下：
NAO 1 2 3 4 5 6 7
---------- ------- ------- ------- ------- ------- ------- -------
1. O 1( 1s ) 2.0000 0.0000 0.0000 0.0000 0.0000 0.1870 0.1870
2. O 1( 2s ) 0.0000 1.9767 0.0000 0.0000 0.0759 1.0236 1.0236
3. O 1( 2px) 0.0000 0.0000 2.0000 0.0000 0.0000 0.0000 0.0000
4. O 1( 2py) 0.0000 0.0000 0.0000 1.5602 0.0000 0.9391 -0.9391
5. O 1( 2pz) 0.0000 0.0759 0.0000 0.0000 1.7530 -0.7415 -0.7415
6. H 2( 1s ) 0.1870 1.0236 0.0000 0.9391 -0.7415 1.4612 0.3307
7. H 3( 1s ) 0.1870 1.0236 0.0000 -0.9391 -0.7415 0.3307 1.4612
然而，不幸的是，所得Mayer键级矩阵很糟
1 2 3
1 0.00000000 7.23697119 7.23697119
2 7.23697119 0.00000000 1.06885215
3 7.23697119 1.06885215 0.00000000
键级大得出奇。NBO5.0输出的PNAO密度矩阵估计错误的，理论上Tr(PS)应当等于分子总电子数，然而用PNAO的密度矩阵和重叠矩阵竟得到16.22606912。虽然测试这个方法由于NBO程序的bug失败了，但原理上是合理的。