利用布居分析判断基函数与原子轨道的对应关系
利用布居分析判断基函数与原子轨道的对应关系
Determining correspondence between basis functions and atomic orbitals via population analysis
文/Sobereva @北京科音 2018-May-26
0 前言
波函数分析程序Multiwfn (http://sobereva.com/multiwfn)的很多分析，比如计算分子轨道中某些原子轨道所占成分（见《谈谈轨道成份的计算方法》http://sobereva.com/131）、绘制PDOS图（见《使用Multiwfn绘制态密度(DOS)图考察电子结构》http://sobereva.com/482）等等都需要在理解基函数与原子轨道的对应关系的基础上定义片段，把需要讨论的原子轨道所对应的基函数纳入片段。有些基组的这种对应关系根据基组名就很容易判断，比如6-31G*，对于碳原子，一看就知道有3个S基函数，第一个收缩度为6的S基函数对应1s原子轨道，而另外两个，即一个收缩度为3和一个未收缩的S基函数一起对应2s原子轨道。但是有些基组不好判断，或者可能有时还会被基组的形式上的名称误导，比如def2系列、Dunning相关一致性基组、UGBS等。实际上利用Multiwfn的布居分析查看基函数壳层的布居数和自旋布居，多数情况可判断出大致对应关系，本文就举一些例子。元素和基组种类繁多，本文不可能都讨论全，故在搞懂本文例子基础上一定要举一反三。如果对自旋布居概念不懂的话，看此文了解常识《谈谈自旋密度、自旋布居以及在Multiwfn中的绘制和计算》（http://sobereva.com/353）。不了解Multiwfn的话看《Multiwfn入门tips》（http://sobereva.com/167）和《Multiwfn波函数分析程序的意义、功能与用途》（http://sobereva.com/184）。本文使用Multiwfn 3.6(dev)版和Gaussian 16 A.03。注意下文中大写的字母比如S,P,D代表基函数壳层，小写字母比如s,p,d代表实际原子轨道壳层。
1 例1：判断6-31G*对于硫原子的情况
硫原子的基态是三重态，其电子结构为1s2 2s2 2p6 3s2 3p4。用Gaussian通过# b3lyp/6-31G*关键词算一个三重态硫原子，然后把.fch载入Multiwfn，之后依次输入7 //布居分析
5 //Mulliken分析
1 //输出Mulliken布居分析结果
从输出信息中我们可以看到各个基函数壳层的分析结果，其中Alpha_pop.、Beta_pop.、Total_pop.、Spin_pop.分别是壳层上的Alpha布居数、Beta布居数、总布居数、自旋布居数
Population of shells:
Shell Type Atom Alpha_pop. Beta_pop. Total_pop. Spin_pop.
1 S 1(S ) 0.99933 0.99934 1.99867 -0.00000
2 S 1(S ) 0.99463 0.99441 1.98904 0.00022
3 P 1(S ) 2.98942 2.97832 5.96773 0.01110
4 S 1(S ) 0.76405 0.76894 1.53299 -0.00490
5 P 1(S ) 2.05660 0.66106 2.71766 1.39553
6 S 1(S ) 0.25905 0.28410 0.54315 -0.02505
7 P 1(S ) 0.95399 0.36062 1.31461 0.59337
8 D 1(S ) -0.01706 -0.04679 -0.06385 0.02973
习俗上，基组在定义时都是按照基函数的主体分布由内（接近原子核）到外（远离原子核）排序的，因此与原子轨道壳层的对应关系在指认时应当也是从内到外的顺序。
要判断哪些P壳层对应硫的2p、3p可以看自旋布居。此原子的两个单电子都在3p上，而5P、7P上的自旋布居加和（1.39553+0.59337）几乎精确为2.0，因此可知5P与7P一起描述3p。而2p至少得有一个P壳层描述，故肯定对应于剩下的3P。
再看S壳层的情况。此体系中s轨道上没有单电子，因此不可能利用自旋布居来判断对应关系，但可以通过总布居数判断。我们可以发现1S和2S上的电子数都几乎精确为2.0，因此分别对应1s和2s，而且4S和6S上的总电子数1.53299+0.54315也几乎为2.0，因此肯定对应3s。
通过以上方式判断出的对应关系，和6-31G*的定义完全相符。下面再看更复杂的情况。
2 例2：判断def2-QZVP对于硫原子的情况
以下是对B3LYP/def2-QZVP计算的硫原子的fch文件按照上一节的做法做布居分析得到的结果Shell Type Atom Alpha_pop. Beta_pop. Total_pop. Spin_pop
1 S 1(S ) 0.82869 0.82873 1.65743 -0.00004
2 S 1(S ) -0.00037 -0.00032 -0.00069 -0.00004
3 S 1(S ) 0.14125 0.14127 0.28252 -0.00002
4 S 1(S ) 0.16419 0.16270 0.32689 0.00149
5 S 1(S ) 0.52143 0.52249 1.04392 -0.00106
6 S 1(S ) 0.30816 0.31310 0.62125 -0.00494
7 S 1(S ) 0.24370 0.24143 0.48513 0.00227
8 S 1(S ) 0.58979 0.55149 1.14128 0.03829
9 S 1(S ) 0.20191 0.23879 0.44069 -0.03688
10 P 1(S ) 2.92983 2.89745 5.82728 0.03238
11 P 1(S ) -0.03547 -0.03299 -0.06847 -0.00248
12 P 1(S ) 0.34015 0.21380 0.55395 0.12635
13 P 1(S ) 1.04061 0.31596 1.35656 0.72465
14 P 1(S ) 1.30799 0.41974 1.72772 0.88825
15 P 1(S ) 0.41602 0.18605 0.60206 0.22997
16 D 1(S ) 0.00002 0.00000 0.00002 0.00001
17 D 1(S ) 0.00028 0.00006 0.00034 0.00022
18 D 1(S ) 0.00070 0.00015 0.00085 0.00055
19 D 1(S ) 0.00027 0.00010 0.00036 0.00017
20 F 1(S ) 0.00063 0.00000 0.00063 0.00063
21 F 1(S ) 0.00025 0.00000 0.00025 0.00025
22 G 1(S ) 0.00000 0.00000 0.00000 -0.00000
先看P的情况。我们发现12P~15P的自旋布居数加和为1.969，和3p上有俩单电子的实际情况一致，因此可以判定12P~15P一起描述3p。虽然这4个P壳层电子数加和为4.24，看似比3p上本来有的4个电子数目要多，但是如果不计入12P，会发现和实际情况偏离明显更大。之所以4.24比本应有的4.0大一些，一方面是Mulliken布居分析本来就不严格（本来布居分析也没有绝对严格的），另一方面基函数本身径向特征就和原子轨道存在差异，而且搞def2基组的人的思路是为了令能量计算误差较小，而并未去刻意保持基函数与原子轨道的严格对应关系。
其它两个P壳层，即10P和11P，应当认为对应的是2p，确实二者的电子数之和接近实际2p上的电子数（6）。值得一提的是，11P上的电子数很接近于0，看似把它视为描述2p还是描述3p都可以，但考虑到从此基组的名字上看，此基组是4-zeta基组，因此应认为只有12P~15P这4个P才对应3p。
接下来看S的情况，还是从总布居数上判断对应关系。由于7S,8S,9S的电子数加和为2.0671，正好和3s上本来有的两个电子相对应，而如果再把6S算进去就会严重高估，因此应当认为此基组是用这三个S基函数描述的3s。看似这和此基组名义上的4-zeta有异，但这就是事实。然后我们看到4S,5S,6S的电子数加和也将近2.0，1S,2S,3S的电子数加和也接近2.0，因此1s、2s和3s在这个基组中都是通过3个S基函数来描述的。
3 例3：判断UGBS对于碳原子的情况
UGBS基组比较凶残，是完全未收缩的基组，基函数数目极多。用B3LYP/UGBS计算基态的碳原子（是三重态），布居分析结果如下Shell Type Atom Alpha_pop. Beta_pop. Total_pop. Spin_pop.
1 S 1(C ) 0.00000 0.00000 0.00000 0.00000
2 S 1(C ) -0.00000 -0.00000 -0.00000 -0.00000
3 S 1(C ) 0.00000 0.00000 0.00000 0.00000
4 S 1(C ) 0.00000 0.00000 0.00000 -0.00000
5 S 1(C ) 0.00000 0.00000 0.00001 0.00000
6 S 1(C ) 0.00001 0.00001 0.00001 -0.00000
7 S 1(C ) 0.00003 0.00003 0.00006 0.00000
8 S 1(C ) 0.00010 0.00010 0.00020 -0.00000
9 S 1(C ) 0.00038 0.00037 0.00075 0.00000
10 S 1(C ) 0.00129 0.00130 0.00259 -0.00000
11 S 1(C ) 0.00437 0.00435 0.00872 0.00001
12 S 1(C ) 0.01374 0.01375 0.02749 -0.00001
13 S 1(C ) 0.03943 0.03941 0.07884 0.00003
14 S 1(C ) 0.09792 0.09762 0.19555 0.00030
15 S 1(C ) 0.19619 0.19752 0.39372 -0.00133
16 S 1(C ) 0.28956 0.28668 0.57625 0.00288
17 S 1(C ) 0.24540 0.25141 0.49682 -0.00601
18 S 1(C ) 0.08091 0.07676 0.15767 0.00416
19 S 1(C ) 0.06111 0.06208 0.12319 -0.00096
20 S 1(C ) 0.25654 0.22257 0.47911 0.03397
21 S 1(C ) 0.39966 0.37481 0.77447 0.02485
22 S 1(C ) 0.25586 0.27426 0.53012 -0.01839
23 S 1(C ) 0.05747 0.09696 0.15444 -0.03949
...
P壳层的情况不用管，因为这体系p壳层只有2p一个，所以所有P壳层都可以当做描述2p，我们主要得需要区分哪些S描述1s哪些描述2s。如果从后往前对总电子数累加，会发现从S19~S23的电子数加和为2.06，看起来可以被视为是用来描述2s的样子。
当前例子比较复杂，S基函数极多，而且电子数分布比较分散。如果想更进一步确认上述归属是否合理，我们可以换个组态再计算试试。我们改成计算碳的五重态，此时电子组态是1s2 2s1 2p3，分析结果如下
Shell Type Atom Alpha_pop. Beta_pop. Total_pop. Spin_pop.
1 S 1(C ) 0.00000 0.00000 0.00000 0.00000
2 S 1(C ) -0.00000 -0.00000 -0.00000 -0.00000
3 S 1(C ) 0.00000 0.00000 0.00000 0.00000
4 S 1(C ) 0.00000 0.00000 0.00000 0.00000
5 S 1(C ) 0.00000 0.00000 0.00001 0.00000
6 S 1(C ) 0.00001 0.00001 0.00001 0.00000
7 S 1(C ) 0.00003 0.00003 0.00006 0.00000
8 S 1(C ) 0.00010 0.00010 0.00020 0.00001
9 S 1(C ) 0.00038 0.00036 0.00074 0.00002
10 S 1(C ) 0.00129 0.00122 0.00251 0.00007
11 S 1(C ) 0.00440 0.00422 0.00862 0.00019
12 S 1(C ) 0.01371 0.01298 0.02669 0.00072
13 S 1(C ) 0.03972 0.03813 0.07785 0.00159
14 S 1(C ) 0.09793 0.09265 0.19058 0.00528
15 S 1(C ) 0.19706 0.19205 0.38912 0.00501
16 S 1(C ) 0.28975 0.27928 0.56903 0.01047
17 S 1(C ) 0.24237 0.25911 0.50148 -0.01674
18 S 1(C ) 0.08282 0.10805 0.19087 -0.02523
19 S 1(C ) 0.05933 0.01147 0.07081 0.04786
20 S 1(C ) 0.26827 0.00043 0.26870 0.26784
21 S 1(C ) 0.41362 -0.00013 0.41349 0.41375
22 S 1(C ) 0.24255 0.00004 0.24259 0.24251
23 S 1(C ) 0.04666 -0.00001 0.04665 0.04666
...
对S19~S23电子数加和，结果为1.04224，对自旋布居数加和，结果为0.97665，都很接近于1.0，这和此时碳原子2s只有一个电子的事实完全相符，而且要把S18算进去，偏离1.0会较明显。因此有很强理由认为S19~S23描述2s，而S1~S18描述1s。
4 例4：判断def2-TZVP对于金原子的情况
def2-TZVP从第五周期开始是赝势基组，搭配的是Stuttgart小核赝势，对Au来说60个内核电子被赝势代替，只有价电子5s,5p,5d,6s被赝势基组表达出来。如果对赝势和赝势基组不熟悉，参看《谈谈赝势基组的选用》（http://sobereva.com/373）及其中的引文。用B3LYP/def2-TZVP计算金原子的基态（二重态），布居分析结果如下Shell Type Atom Alpha_pop. Beta_pop. Total_pop. Spin_pop.
1 S 1(Au) 0.01764 0.01588 0.03352 0.00175
2 S 1(Au) -0.25037 -0.22638 -0.47675 -0.02399
3 S 1(Au) 0.90648 0.84814 1.75462 0.05834
4 S 1(Au) 0.33047 0.35849 0.68895 -0.02802
5 S 1(Au) 0.60949 0.00436 0.61385 0.60513
6 S 1(Au) 0.38630 -0.00049 0.38581 0.38679
7 P 1(Au) 1.32302 1.32562 2.64864 -0.00260
8 P 1(Au) 1.43496 1.44146 2.87642 -0.00650
9 P 1(Au) 0.24145 0.23250 0.47395 0.00896
10 P 1(Au) 0.00057 0.00042 0.00099 0.00015
11 D 1(Au) 3.11664 3.17522 6.29186 -0.05859
12 D 1(Au) 1.48961 1.45237 2.94198 0.03724
13 D 1(Au) 0.39375 0.37241 0.76616 0.02135
14 F 1(Au) 0.00000 0.00000 0.00000 0.00000
金原子的组态是5s2 5p6 5d10 6s1，显然7P~10P对应5p，11D~13D对应5d。5S和6S的电子数及自旋布居数加和几乎都精确为1.0，正好对应6s有一个电子。而把1S~4S的电子数加和几乎精确为2.0，对应5s的两个电子。显然1S~4S描述的是5s，而5S和6S描述的是6s。