利用布居分析判断基函数与原子轨道的对应关系

sobereva

利用布居分析判断基函数与原子轨道的对应关系

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。

flying

谢谢分享，学习了。之前一直不明白基函数与原子轨道关系，这下了解了。

Uus/pMeC6H4-/キ

对6-31G*算硫的例子的疑问：
1. 那个单独的D壳层基函数是怎么来的呢？如果对应3d价层，按理2-zeta基组应该有2个基函数；如果对应d极化函数，那一层d极化函数应该有5个（球谐型）或6个（笛卡尔型）基函数才对。
2. 在获得“5P与7P一起描述3p，4S和6S一起描述3s”的结论后，有可能进一步区分基函数中哪个收缩度为3、哪个收缩度为1么？
3. 所谓“SP壳层”（也有标记为L的）可能导致基组定义文件里写的一组“指数与收缩系数”实际对应基组描述的两个壳层的基函数，那么这类壳层是Pople系列基组特有的，还是也有其他主流的非广义收缩的基组会用？

sobereva

Uus/pMeC6H4-/キ 发表于 2025-5-17 11:00
对6-31G*算硫的例子的疑问：
1. 那个单独的D壳层基函数是怎么来的呢？如果对应3d价层，按理2-zeta基组应该 ...

1 看不懂。D壳层要带就是带一层，不可能只有2个D基函数
6-31G*的D壳层对硫是极化函数

2 5P、4S基函数壳层收缩度为3。实际上按照基组定义就知道（体现了基函数的顺序）

3 Pople系列特有的

Uus/pMeC6H4-/キ

sobereva 发表于 2025-5-17 22:22
1 看不懂。D壳层要带就是带一层，不可能只有2个D基函数
6-31G*的D壳层对硫是极化函数

谢谢。第1点说的“单独的D”是指列表中的8D，总布居-0.06385的那个壳层，我昨天的困惑可能是把基函数和壳层弄混了。定义硫6-31G*基组的文件只给出一套D壳层的GTH指数和收缩系数，文中布居分析结果的列表里也只有8D一个是D壳层；但是作为极化函数的D壳层有5或6个基函数，所以文中8D一个极化函数的壳层就同时对应了5或6个基函数，而不像文中5P、7P、4S、6S等描述价层的壳层那样一个壳层对应一个基函数，这么总结对吗？

另外还想再加几个问题：
4. 之前http://sobereva.com/235计算单原子波函数时，考虑了B3LYP、HF等理论方法能否正确描述基态原子组态，那么此文中做布居分析用的理论方法是否也需要考虑呢？对同一原子的同一基组，用不同理论方法按此文中操作判断的对应关系能不能通用，比如借助B3LYP的结果解释HF、BLYP、CAM-B3LYP、B2PLYP甚至MP2的计算结果？
5. 文中判断的对应关系仅精确到1s、2p、3d等特定主量子数n与角量子数l的原子轨道，如果还需要精确到2pz、3dz^2等特定磁量子数m_l的原子轨道，具体如何操作呢？比如不同体系的轨道成分分析横向对比，计算前已考虑http://sobereva.com/52提到的问题统一到贡献最大化的朝向，后续也按http://sobereva.com/482定义了特定基函数类型标识来绘制态密度图，此时就有这个需要。

sobereva

Uus/pMeC6H4-/キ 发表于 2025-5-18 12:05
谢谢。第1点说的“单独的D”是指列表中的8D，总布居-0.06385的那个壳层，我昨天的困惑可能是把基函数和壳 ...

此文给出的都是壳层。D是你说的那样。一个P基函数壳层也有X,Y,Z共三个基函数。仅有S才是一个壳层一个基函数

4 用一个对当前体系合适的泛函做的布居分析，就可以判断基函数与原子轨道的大致对应关系，这个对应关系也可以套用到其它计算级别产生的波函数并得到相应的原子轨道布居数

5 得到了基函数壳层与原子轨道壳层对应的关系，自然也就得到了相应壳层里特定类型的基函数（如XY）与原子轨道（如dxy）的对应关系

Uus/pMeC6H4-/キ

今天看http://bbs.keinsci.com/thread-766-1-1.html一帖，提到multiwfn的隐藏选项1000有个子选项15 Make orbitals equivalent to basis functions可以把轨道和基函数一一对应，再进主功能0可以直接看基函数图形。仍用例1的三重态硫原子，用#p B3LYP/6-31G(d)关键词（不得不说这个基组默认6D 7F的笛卡尔型/球谐型混用实在奇葩），这样处理得到19个alpha轨道（提示Now each alpha orbital corresponds to a basis function, index of alpha orbitals is identical to index of basis functions）

1.          1      E(au/eV):   -88.93372   -2420.0097 Occ: 1.000000 Typ: A
   
2.          2      E(au/eV):    -7.99324    -217.5070 Occ: 1.000000 Typ: A
   
3.          3      E(au/eV):    -5.96283    -162.2570 Occ: 1.000000 Typ: A
   
4.          4      E(au/eV):    -5.96283    -162.2570 Occ: 1.000000 Typ: A
   
5.          5      E(au/eV):    -5.94344    -161.7292 Occ: 1.000000 Typ: A
   
6.          6      E(au/eV):    -0.72314     -19.6776 Occ: 1.000000 Typ: A
   
7.          7      E(au/eV):    -0.33262      -9.0509 Occ: 1.000000 Typ: A
   
8.          8      E(au/eV):    -0.33262      -9.0509 Occ: 1.000000 Typ: A
   
9.          9      E(au/eV):    -0.28443      -7.7398 Occ: 1.000000 Typ: A
   
10.         10      E(au/eV):     0.25230       6.8655 Occ: 0.000000 Typ: A
    
11.         11      E(au/eV):     0.32670       8.8899 Occ: 0.000000 Typ: A
    
12.         12      E(au/eV):     0.32670       8.8899 Occ: 0.000000 Typ: A
    
13.         13      E(au/eV):     0.35291       9.6033 Occ: 0.000000 Typ: A
    
14.         14      E(au/eV):     0.68195      18.5568 Occ: 0.000000 Typ: A
    
15.         15      E(au/eV):     0.68195      18.5568 Occ: 0.000000 Typ: A
    
16.         16      E(au/eV):     0.72247      19.6595 Occ: 0.000000 Typ: A
    
17.         17      E(au/eV):     0.72247      19.6595 Occ: 0.000000 Typ: A
    
18.         18      E(au/eV):     0.73738      20.0653 Occ: 0.000000 Typ: A
    
19.         19      E(au/eV):     3.82638     104.1212 Occ: 0.000000 Typ: A

复制代码

观察图形尝试对应原子轨道的结果：1-1s, 2-2s, 3-2px, 4-2py, 5-2pz, 6-3s(1), 7-3px(1), 8-3py(1), 9-3pz(1), 10-3s(2), 11-3px(2), 12-3py(2), 13-3pz(2), 14-dxx, 15-dyy, 16-dzz, 17-dxy, 18-dxz, 19-dyz，很合理。结合Levine的Quantum Chemistry (7 ed.)第15章15.4节的讲解，也能把基组定义文件里的收缩系数和指数对应到具体原子轨道上。
与此同时，用#p B3LYP/6-31G(d) 5D关键词算的结果在iloadasCart= 0的设置下载入multiwfn并做相同操作，得到的最后几个轨道变为

1.         14      E(au/eV):     0.68165      18.5486 Occ: 0.000000 Typ: A
   
2.         15      E(au/eV):     0.68165      18.5486 Occ: 0.000000 Typ: A
   
3.         16      E(au/eV):     0.72218      19.6516 Occ: 0.000000 Typ: A
   
4.         17      E(au/eV):     0.72218      19.6516 Occ: 0.000000 Typ: A
   
5.         18      E(au/eV):     0.73713      20.0583 Occ: 0.000000 Typ: A

复制代码

对应关系14-dz^2, 15-dxz, 16-dyz, 17-d(x^2-y^2), 18-dxy也能看出来。

编辑：然而对于例2的def2-QZVP（Gaussian默认5D 7F，对硫是96421即9s6p4d2f1g），这样会一下产生70个轨道，再这样尝试对应就复杂了