谈谈轨道成份的计算方法
注：本文属于深入、详细的讨论。如果只是着急想马上得到轨道成分计算结果，直接看本文第6节或Multiwfn手册的4.8节的例子即可。
谈谈轨道成份的计算方法
On the calculation methods of orbital composition
文/Sobereva @北京科音
First release: 2012-Feb-24 Last update: 2020-Jun-3
1 前言
经常有人发帖问怎么计算轨道成份，文献中也经常给出轨道成份，用于讨论电子结构、分析电荷转移、根据前线轨道理论分析反应位点等等。然而轨道成份的计算方法却极少被讨论，时常看到错误的说法在网上以讹传讹。在《分子轨道成份的计算》（化学学报, 69, 2393）一文中作者深入详细地讨论了轨道成份计算方法的原理并通过实际体系进行了比较，此文可以在此下载：/usr/uploads/file/20150610/20150610022849_33146.pdf。但是对于初学者来说可能一些讨论并不很容易理解，而且程序的具体使用步骤未在其中涉及。本文算是那篇文章的补充，本文对各种计算方法的原理介绍比较粗略，忽略了不少细节讨论，但是对初学者来说更为易懂，并会详细说明实际计算时的过程，重点说明了Multiwfn中分析轨道成份功能的基本用法。希望对于怀有“轨道成份怎么计算”问题的人在看完此文以及《分子轨道成份的计算》后都能找到答案。
笔者的另一篇文章《通过轨道离域指数(ODI)衡量轨道的空间离域程度》（http://sobereva.com/525）和此文内容关系较密切，其中笔者提出了一种名为ODI的指数，基于轨道中的原子所占成份来计算，可以很可靠地衡量轨道离域程度，感兴趣的读者建议一看。
提醒：按照本文使用Multiwfn计算轨道成份的数据发表文章时必须按照Multiwfn启动时显示的信息引用Multiwfn原文，也请同时这样引用前述的《分子轨道成份的计算》: Tian Lu, Feiwu Chen, Calculation of Molecular Orbital Composition, Acta Chim. Sinica, 69, 2393-2406 (2011) (in Chinese)
2 轨道成份的概念
计算轨道成份也就是计算构成轨道的各个组成部分的贡献。层次从小到大可以分为：(1)计算基函数的贡献 (2)计算原子轨道的贡献 (3)计算原子的贡献 (4)计算某个片段的贡献。
正如同蛋糕有不同切法，计算轨道成份的方法也有很多种，没有一种方法是绝对正确的，毕竟轨道成份不是可观测量，定义完全是人为的。很多文献中给出轨道成份时不注明是怎么计算的，这是不应该的，这导致数据难以被后人重现，也没法得知作者所用的轨道成份计算方法及结果是否合理。
不管什么成份计算方法，最基本的要求是能够归一，即所有部分的贡献和为100%。另一个基本要求是成份必须是正值，负值是没有意义的。在满足这两个要求下，还应当尽量满足其它一些要求：(1)满足旋转不变性。因为分子的朝向是任意的，不能因为分子在空间中旋转了，成份的计算结果产生过大变化。(2)结果应当能随基组增大而收敛，否则相当于没法得到一个确定的结果 (3)计算速度应当快，算法应当简单，易于实现 (4)物理意义应当清晰 (5)结果应当和分子轨道图形定性上相符 (6)适用体系应当广泛 (7)数值稳定性好，避免出现异常值。
接下来谈谈计算基函数、原子轨道、原子以及片段的成份的基本原理和具体计算步骤。本文所指的“轨道”可以是任何类型轨道，比如分子轨道、定域化轨道、自然轨道等等，但举例都用分子轨道举例，对于其它类型轨道的处理方法是一样的。
3 计算基函数所占的成份
基函数本身没有直接的化学意义，它只是求解量子化学问题的数学工具。之所以要提到它是因为它可以认为是构成轨道的最小单元，在计算轨道中其它组成部分的成份时会涉及到它。另外，对于诸如STO-3G这种极小基来说，一个基函数就对应于一个原子轨道，计算基函数的成份此时就是计算原子轨道的成份，是有化学意义的。
这里先来考虑一个最简单的异核双原子体系HeH+，用Gaussian在HF/STO-3G计算后的pop=full输出是
Molecular Orbital Coefficients:
1 2
O V
Eigenvalues -- -1.52378 -0.26764
1 1 He 1S 0.90013 -0.64981
2 2 H 1S 0.19438 1.09302
可知占据的那条分子轨道波函数可写为ψ1=0.9*χ_1 + 0.1944*χ_2。其中χ是基函数符号。根据轨道的归一化条件可写为
1=<ψ1|ψ1>=0.9^2*<χ_1|χ_1>+2*0.9*0.1944*<χ_1|χ_2>+0.1944^2*<χ_2|χ_2>=0.9^2*S11+2*0.9*0.1944*S12+0.1944^2*S22。
上式中，0.9^2和0.1944^2分别是χ_1和χ_2的“定域项”，表明了它们独自对ψ1的贡献。2*0.9*0.1944*S12是交叉项，代表了χ_1和χ_2对ψ1的联合贡献。
由于基函数是归一化了的，所以S11和S22都为1，因此<ψ1|ψ1>=0.8478+2*0.9*0.1944*S12。很显然，很多人用轨道系数的平方来计算轨道成份是大错特错，0.9^2+0.1944^2等于0.8478而明显不得1！即这么算出来的轨道成份和不是100%。这是因为基函数间不是正交的，重叠积分S12不为零，所以计算轨道成份必须考虑怎么把交叉项划分给各个基函数。
计算基函数成份的方法有三种，分别是Mulliken、Stout-Politzer和SCPA方法。本质上，这三种方式的差异仅在于怎么划分交叉项。
Mulliken方法在计算χ_n对轨道的贡献时，将χ_n的定域项全部算作χ_n的贡献，而将χ_n与其它基函数间的交叉项的一半划给χ_n。对于上例，He对轨道的贡献是(0.9^2+0.9*0.1944*S12)*100%。利用iop(3/33=1)关键词可输出重叠矩阵可知S12=0.4343，因此He的贡献为88.6%。实际上，计算轨道成份的Mulliken方法就对应于Mulliken布居用到的方法。Mulliken布居方法中，某条轨道上某基函数的布居数就是这个基函数的贡献值乘上轨道占据数，因此χ_1在ψ1这条双占据轨道上的Mulliken布居数就是88.6%*2=1.772。由于以He为中心的只有χ_1这么一个基函数，而且此体系只有ψ1是占据轨道，因此He原子的Mulliken布居数就是1.772。由于He的核电荷为2，因此它的Mulliken电荷是2-1.772=0.228，和Gaussian输出的一致。
Mulliken方法对交叉项的平分做法被一些人批评，因为这显得没物理意义。于是有些人提出了改进方法。Stout-Politzer方法考虑了基函数的不均等性，它将交叉项按照相应两个基函数的系数平方比来划分。对于上例，χ_1所划得的交叉项就是2*0.9*0.1944*S12* 0.9^2/(0.9^2+0.1944^2)。虽然Stout-Politzer方法看起来比Mulliken方法考虑得更周到，但是实际结果反倒往往有所恶化，数值不稳定性增加，出现负值的倾向加剧，还破坏了Mulliken方法的旋转不变性。此方法不建议使用。
很多人在用系数平方来计算轨道成份时也注意到了必须满足归一条件，于是就乘了个归一化因子，因此χ_1的贡献就是0.9^2/(0.9^2+0.1944^2)*100%。这个做法实际上对应于SCPA布居方法所用到的轨道划分方法。这个方法的计算公式经过简单推导（见《分子轨道成份的计算》），可知本质上和Mulliken、Stout-Politzer方法的差异也仅是划分交叉项的做法不同而已。Mulliken和Stout-Politzer时常出现没意义的负的成份值，对于虚轨道尤为严重，而SCPA方法的好处是结果一定为正。另外，SCPA方法计算简单，不需要考虑重叠积分，通过比较发现数值合理性也比Mulliken和Stout-Politzer方法都略强，是推荐使用的计算基函数贡献的方法。
4 计算原子轨道所占的成份
在分子体系中实际上没有清晰的原子轨道的概念，但是以某些方法将原子轨道概念还原出来，并计算其在各轨道中的成份，有助于分析电子结构特征、成键本质。
4.1 基函数加和方法
第一种计算方法是利用基组与原子轨道的对应关系。对于极小基，每个基函数描述一个原子轨道，因此就按上一节讨论的方法就得到了原子轨道的贡献。而对于扩展基，则是以一个或多个基函数来描述一个原子轨道体现的效应，因此可以把原子轨道对应的基函数的贡献相加作为原子轨道的贡献。
例如对于HF/6-31G*下计算水分子，Gaussian输出的占据轨道系数是这样的：
1 2 3 4 5
(A1)--O (A1)--O (B2)--O (A1)--O (B1)--O
Eigenvalues -- -20.55790 -1.34610 -0.71418 -0.57083 -0.49821
1 1 O 1S 0.99462 -0.20953 0.00000 -0.07310 0.00000
2 2S 0.02117 0.47576 0.00000 0.16367 0.00000
3 2PX 0.00000 0.00000 0.00000 0.00000 0.63927
4 2PY 0.00000 0.00000 0.50891 0.00000 0.00000
5 2PZ -0.00134 -0.09475 0.00000 0.55774 0.00000
6 3S 0.00415 0.43535 0.00000 0.32546 0.00000
7 3PX 0.00000 0.00000 0.00000 0.00000 0.51183
8 3PY 0.00000 0.00000 0.30390 0.00000 0.00000
9 3PZ 0.00046 -0.04982 0.00000 0.40482 0.00000
10 4XX -0.00394 -0.00103 0.00000 0.01187 0.00000
11 4YY -0.00421 0.02692 0.00000 0.00079 0.00000
12 4ZZ -0.00409 0.02132 0.00000 -0.04630 0.00000
13 4XY 0.00000 0.00000 0.00000 0.00000 0.00000
14 4XZ 0.00000 0.00000 0.00000 0.00000 -0.03417
15 4YZ 0.00000 0.00000 -0.05088 0.00000 0.00000
16 2 H 1S 0.00032 0.13302 0.23243 -0.14007 0.00000
17 2S -0.00021 0.00173 0.10729 -0.08280 0.00000
18 3 H 1S 0.00032 0.13302 -0.23243 -0.14007 0.00000
19 2S -0.00021 0.00173 -0.10729 -0.08280 0.00000
注意这样的表里每行对应一个基函数，而基函数前的数字编号没有物理意义，它只是为了区分以相同原子为中心的类型相同的基函数，编号越靠前的通常分布区域越接近原子核。这些行并非直接对应原子轨道，数字也不是原子轨道的主量子数。很多初学者容易在这点上混淆。
对于6-31G*基组，每个内层原子轨道用一个收缩度为6的基函数来体现，故这个基函数的贡献就是相应内层原子轨道的贡献。也就是说，氧的1s原子轨道成份就是这里的1S基函数的成份（本文中原子轨道符号用小写字母，基函数符号用大写字母）。而6-31G*中每个价层原子轨道用一个收缩度为3和一个收缩度为1的相应类型的基函数一起来表现，这两个基函数的贡献之和可认为是对应的原子轨道的贡献。也就是说，氧的2px原子轨道的贡献可认为是这里2PX和3PX基函数的贡献之和，而氢的1s原子轨道的贡献就是其1S和2S基函数贡献之和。表中的4XX、4YY等等这些D型基函数没有物理意义，显然不可能对应于本来就不存在的2d原子轨道，这些D型基函数存在的意义只是为了在自洽场计算中更充分地描述出实际的电子密度的极化效果
按照基函数与原子轨道对应关系来计算原子轨道贡献的方法有很大局限性。有些基组基函数很多，不容易搞清楚对应关系（但可以巧妙利用Multiwfn的布居分析获得对应关系，见《利用布居分析判断基函数与原子轨道的对应关系》http://sobereva.com/418）。对于弥散、极化函数也没法处理，考虑它们的话没法确定应该将它们归在哪个原子轨道里，而不考虑的话则原子轨道贡献和不为100%（但可以对已考虑了的原子轨道贡献值做归一化处理）。而且由于基函数之间的混合、缺乏对基函数在原子轨道中权重的考虑，导致这种方法在原理上不甚严格。虽然缺点能列举出很多，但实际应用中这种计算方法还是比较便利的，对于大多数情况来说没什么问题。但并不建议在含有弥散基函数的情况下，或者基组不平衡较强（即某些原子的基组质量和相邻原子基组质量相差悬殊。如某原子用6-311G(2d)，而旁边的原子只用3-21G）时使用本方法。由于用这种方法计算原子轨道贡献时先要计算基函数的贡献，这就牵扯到用什么方法来计算基函数的贡献，我仍建议用SCPA方法。
用本节的方法对于计算d及以上角动量原子轨道的贡献时有一点需要注意，也就是应当尽量用球谐型高斯函数，而不建议用笛卡尔型高斯函数。这点不打算在这里多谈，请参阅《谈谈5d、6d型d壳层基函数与它们在Gaussian中的标识》一文（http://sobereva.com/51）以及《分子轨道成份的计算》中的讨论。
4.2 NAO方法
另一种计算原子轨道成份的方法是借助于NBO分析。NBO分析可以从体系的密度矩阵信息中得到自然原子轨道(NAO)，NAO的数目与原先基函数相同。NAO当中被称为“极小集”的那部分就对应于一般意义上的原子轨道。并且NBO分析还可以将分子轨道以NAO来展开。由于NAO彼此间都是正交的，不用考虑怎么划分交叉项，因此NAO在相应轨道中的系数平方就是它在轨道中的成份了。
在Gaussian中，用了pop=nboread，输入文件末尾空一行写上$NBO NAOMO $END，在NBO模块中就会输出以NAO为基的MO的展开系数，还是以HF/6-31G*的水为例：
MOs in the NAO basis:
NAO 1 2 3 4 5 6 7 8
---------- ------- ------- ------- ------- ------- ------- ------- -------
1. O 1 (S) 0.9958 -0.0896 0.0000 -0.0177 0.0000 0.0046 0.0000 0.0000
2. O 1 (S) 0.0847 0.8834 0.0000 0.2928 0.0000 0.2288 0.0000 0.0000
3. O 1 (S) 0.0003 0.0081 0.0000 -0.0246 0.0000 -0.2297 0.0000 0.0000
4. O 1 (S) 0.0000 0.0000 0.0000 0.0000 0.0000 0.0554 0.0000 0.0000
5. O 1 (px) 0.0000 0.0000 0.0000 0.0000 0.9993 0.0000 0.0000 0.0000
6. O 1 (px) 0.0000 0.0000 0.0000 0.0000 -0.0173 0.0000 0.0000 0.0000
7. O 1 (py) 0.0000 0.0000 0.8554 0.0000 0.0000 0.0000 -0.3288 -0.3600
8. O 1 (py) 0.0000 0.0000 0.0334 0.0000 0.0000 0.0000 0.0250 0.4768
9. O 1 (pz) 0.0023 -0.1585 0.0000 0.9170 0.0000 -0.2598 0.0000 0.0000
10. O 1 (pz) -0.0011 -0.0006 0.0000 -0.0124 0.0000 0.0254 0.0000 0.0000
11. O 1 (d1) 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
12. O 1 (d2) 0.0000 0.0000 0.0000 0.0000 -0.0342 0.0000 0.0000 0.0000
13. O 1 (d3) 0.0000 0.0000 -0.0388 0.0000 0.0000 0.0000 0.0177 -0.1517
14. O 1 (d4) 0.0001 -0.0176 0.0000 0.0075 0.0000 0.0212 0.0000 0.0000
15. O 1 (d5) 0.0000 0.0057 0.0000 -0.0352 0.0000 0.0106 0.0000 0.0000
16. H 2 (S) 0.0244 0.3046 0.3644 -0.1882 0.0000 -0.4212 0.3712 0.4027
17. H 2 (S) 0.0028 -0.0165 0.0095 -0.0079 0.0000 -0.4839 0.5547 -0.3845
18. H 3 (S) 0.0244 0.3046 -0.3644 -0.1882 0.0000 -0.4212 -0.3712 -0.4027
19. H 3 (S) 0.0028 -0.0165 -0.0095 -0.0079 0.0000 -0.4839 -0.5547 0.3845
下面是在一开始输出的NAO布居分析
NAO Atom No lang Type(AO) Occupancy Energy
----------------------------------------------------------
1 O 1 S Cor( 1S) 1.99992 -20.39645
2 O 1 S Val( 2S) 1.74644 -1.14691
3 O 1 S Ryd( 3S) 0.00135 1.42530
4 O 1 S Ryd( 4S) 0.00000 3.93946
5 O 1 px Val( 2p) 1.99707 -0.49473
6 O 1 px Ryd( 3p) 0.00060 1.16772
7 O 1 py Val( 2p) 1.46327 -0.31260
8 O 1 py Ryd( 3p) 0.00223 1.30664
9 O 1 pz Val( 2p) 1.73207 -0.41806
10 O 1 pz Ryd( 3p) 0.00031 1.19893
11 O 1 dxy Ryd( 3d) 0.00000 2.03061
12 O 1 dxz Ryd( 3d) 0.00234 2.06405
13 O 1 dyz Ryd( 3d) 0.00301 2.91411
14 O 1 dx2y2 Ryd( 3d) 0.00073 2.47969
15 O 1 dz2 Ryd( 3d) 0.00254 2.02361
16 H 2 S Val( 1S) 0.52321 0.33308
17 H 2 S Ryd( 2S) 0.00086 0.70497
18 H 3 S Val( 1S) 0.52321 0.33308
19 H 3 S Ryd( 2S) 0.00086 0.70497
表中Cor对应于内核NAO，Val对应于价层NAO，标着这两类标签的NAO就是极小集部分，可以简单认为它们就是一般意义的原子轨道。而Ryd型NAO是用来描述那些没有被Cor和Val型NAO所描述的信息，它们对占据的分子轨道贡献极小而一般可以忽略，但是对于高能态的虚轨道则可能有很大贡献。通常我们分析的都是占据轨道和能量最低的几个虚轨道，所以一般不必考虑Ryd型NAO。
将上面两个表对照，就能知道怎么计算各个原子轨道成份了。例如，要计算第三个MO中氧的2py轨道的贡献，从上面表中看到py Val( 2p)对应的NAO标号是7，因此就将系数表中第三列第七行的数值0.8554取平方乘100%即可，即73.2%。
4.3 AOIM方法
AOIM方法有两个版本，第一个版本是刘文剑和黎乐民提出的，而这里涉及的是后来由Haigang Lu、黎乐民等人提出的版本。相应的分析程序AOIM可以在http://faculty.sxu.cn/luhg/aoim.html下载（已失效）。
AOIM是通过正交投影的方法在分子体系中还原出原子轨道，每个原子轨道以一个STO函数表达，这个过程中令误差函数最小化来尽可能地避免波函数信息的损失。同时也得到了每个分子轨道向各个原子轨道的展开系数。AOIM相当于将原始的在扩展基下表达的波函数转化为了在极小基下表达，因此按照本文第三节的方法，就可以计算各轨道中每个原子轨道的成份了，在AOIM中这一步用的是Mulliken方法。
AOIM程序需要读入Gaussian的fch文件来进行分析。在DOS下执行比如aoim < ltwd.fch就可以对ltwd.fch进行AOIM分析。还是那个水分子的例子，AOIM输出的轨道成份信息部分是
MO components in correctted STOs(>0.05):
MO: 1 E= -20.55786 A.U.
1.00 1 O1S
MO: 2 E= -1.34606 A.U.
0.81 1 O2S
0.08 2 H1S
0.08 3 H1S
MO: 3 E= -0.71414 A.U.
0.60 1 O2Py
0.20 2 H1S
0.20 3 H1S
MO: 4 E= -0.57083 A.U.
0.09 1 O2S
0.81 1 O2Pz
...（略）
可见第三条MO的氧的2py轨道的贡献是60%。和之前用NAO方法的结果73.2%有一定差异，这是正常的，不同计算方法不可能总在定量上很一致，尤其是轨道涉及多个原子间交叠的情况差异会更明显，但多数情况下本节介绍的三种方法结果还是定性一致的。
AOIM程序只会输出贡献值大于5%的原子轨道，输出的分子轨道只包含占据轨道和LUMO及LUMO+1轨道，如果想考察贡献值更低的基函数以及更高阶的虚轨道，可以自行修改代码，这很容易做到。
NAO和AOIM方法的基组稳定性都明显好于加和基函数贡献的方法，能随着基组质量增加而很快收敛，而且原理都比较严格。可惜AOIM程序年久失修，对Gaussian09的fch文件不兼容，对大体系分析速度慢，不支持f及以上的角动量基函数，也不支持球谐型基函数，且在个别体系比如氧化铍上有严重问题。对于计算原子轨道成份的方法，整体来说我最推荐NAO方法。
5 计算原子以及片段所占的成份
如果只是想定性了解各个原子对轨道贡献的相对大小，最简单的办法就是显示轨道图形，将isovalue调到一个能够区分不同位置轨道波函数数值大小的值，然后图形观察各个原子附近的等值面涵盖范围的大小即可，范围越大表明贡献越大。
而定量计算原子所占成份最直接的办法就是将原子所拥有的全部原子轨道所占成份加和，这里我也是建议用NAO方法来计算原子轨道所占成份。
还有一类基于实空间划分的方法可用于计算原子所占成份。这类方法都是将整个分子空间划分为一个个原子区域，然后积分某个轨道在相应原子空间内的模的平方来作为此原子在此轨道中的成份。一种广为熟知的划分方法就是Bader提出的QTAIM划分方法，是根据电子密度零通量面来划分，但这种划分一方面是积分困难，另一方面也没有很强的化学意义。比较适合用于计算原子贡献的划分方式是Hirshfeld划分，这种划分是模糊的，原子的空间是相互交叠的，没有明确的边界，具体来说就是空间中每个位置上各个原子都占有一定权重。设ρ_pro(r)为分子中每个原子都处在自由状态时电子密度的叠加（也叫Promolecule密度），r代表空间坐标，而ρ_A(r)是A原子在自由状态下的密度，那么A原子的权重函数就是ρ_A(r)/ρ_pro(r)。在离A原子核越近的区域中A原子的权重会越大。这种Hirshfeld划分的化学意义比较直观明确，积分也比较容易。在计算结果上与NAO方法定性一致，但计算量还是要大于NAO方法。对于一些虚轨道，尤其是高能态时，Ryd型NAO所占成份不可忽视，也就是说加和极小集NAO的成份会明显偏离100%，此时就不建议用NAO方法了，而建议用Hirshfeld方法。
利用Becke函数也可以进行模糊式划分，形式上和Hirshfeld划分类似，结果也比较相似。相对于Hirshfeld划分的好处是原子权重不需要通过原子的密度来构建，而只需要提供原子半径即可。Hirshfeld还有种重要的变体叫Hirshfeld-I，它通过迭代方式逐渐调整原子权重函数，使之可以响应原子实际所处的化学环境，基于Hirshfeld-I划分也可以计算原子对轨道的贡献。
计算分子中某个片段的方法没什么可说的，就是将其中原子所占成份加和即可。
6 在Multiwfn中计算轨道成份
在Multiwfn（http://sobereva.com/multiwfn）中的主功能8里提供了分析轨道成份的功能，支持Mulliken、Stout-Politzer、SCPA方法，还独家支持Hirshfeld、Hirshfeld-I、Becke、NAO、AIM方法。功能强大、操作简单，输出信息清晰易读，比其它程序好用得多。
使用Mulliken、Stout-Politzer、SCPA方法进行分析的话必须用包含基函数信息的文件作为输入文件，如.mwfn、.fch、.molden、GAMESS-US输出文件等，不清楚的话详见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。以Becke、Hirshfeld、Hirshfeld-I、AIM方法分析原子贡献的话任何带有波函数信息的文件都可以用（.fch/.molden/.wfn/.wfx等）。以NAO方法分析的话需要提供NBO程序的输出文件。这里还是以分析HF/6-31G*分子轨道成份为例进行说明。选项的详细介绍请阅读手册的3.10节。
6.1 在Multiwfn中使用Mulliken、Stout-Politzer、SCPA方法进行分析
启动Multiwfn，输入.fch或.molden等含有基函数信息的文件路径后，选8就进入了成份分析界面。如果你的经验不丰富的话，可以先选-1进入分子片段设定界面，然后输入all来查看一下当前体系中的所有基函数信息（之后输入q退回）：
Basis: 1 Shell: 1 Center: 1(O ) Type:S
Basis: 2 Shell: 2 Center: 1(O ) Type:S
Basis: 3 Shell: 3 Center: 1(O ) Type:X
Basis: 4 Shell: 3 Center: 1(O ) Type:Y
Basis: 5 Shell: 3 Center: 1(O ) Type:Z
Basis: 6 Shell: 4 Center: 1(O ) Type:S
Basis: 7 Shell: 5 Center: 1(O ) Type:X
Basis: 8 Shell: 5 Center: 1(O ) Type:Y
Basis: 9 Shell: 5 Center: 1(O ) Type:Z
Basis: 10 Shell: 6 Center: 1(O ) Type:XX
Basis: 11 Shell: 6 Center: 1(O ) Type:YY
Basis: 12 Shell: 6 Center: 1(O ) Type:ZZ
Basis: 13 Shell: 6 Center: 1(O ) Type:XY
Basis: 14 Shell: 6 Center: 1(O ) Type:XZ
Basis: 15 Shell: 6 Center: 1(O ) Type:YZ
Basis: 16 Shell: 7 Center: 2(H ) Type:S
Basis: 17 Shell: 8 Center: 2(H ) Type:S
Basis: 18 Shell: 9 Center: 3(H ) Type:S
Basis: 19 Shell: 10 Center: 3(H ) Type:S
每个基函数对应的壳层，所属中心，基函数类型都一目了然。每个壳层包含一组基函数，即S壳层包含一个S基函数，P壳层包含X,Y,Z三个基函数，6D型壳层包含XX,YY,ZZ,XY,XZ,YZ共六个基函数...
比如想用Mulliken方法分析第3条分子轨道的成份，就选择1，然后输入轨道编号3，就得到了以下结果：
Threshold of absolute value: > 0.500000%
Orbital: 3 Energy(a.u.): -0.71414263 Occ: 2.00000000 Type: Alpha&Beta
Basis Type Atom Shell Local Cross term Total
4 Y 1(O ) 3 25.89789% 14.57170% 40.46960%
8 Y 1(O ) 5 9.23602% 16.82836% 26.06437%
15 YZ 1(O ) 6 0.25896% 0.83797% 1.09693%
16 S 2(H ) 7 5.40298% 7.23541% 12.63839%
17 S 2(H ) 8 1.15080% 2.39536% 3.54616%
18 S 3(H ) 9 5.40298% 7.23541% 12.63839%
19 S 3(H ) 10 1.15080% 2.39536% 3.54616%
Sum up listed above: 48.50044% 51.49956% 100.00000%
Sum up all basis functions: 48.50044% 51.49956% 100.00000%
Composition of each shell, threshold of absolute value: > 0.500000%
Shell 3 Type: P in atom 1(O ) : 40.46960%
Shell 5 Type: P in atom 1(O ) : 26.06437%
Shell 6 Type: D in atom 1(O ) : 1.09693%
Shell 7 Type: S in atom 2(H ) : 12.63839%
Shell 8 Type: S in atom 2(H ) : 3.54616%
Shell 9 Type: S in atom 3(H ) : 12.63839%
Shell 10 Type: S in atom 3(H ) : 3.54616%
Composition of each atoms:
Atom 1(O ) : 67.630897%
Atom 2(H ) : 16.184551%
Atom 3(H ) : 16.184551%
默认情况下只有贡献值大于0.5%的项被输出，这个阈值由Multiwfn目录下的settings.ini里的compthres参数控制。由于这条轨道中各个基函数的贡献要么为0要么就大于0.5%，因此Sum up listed above（显示出的基函数的贡献和）的值与Sum up all basis functions（全部基函数的贡献和）相同。
Total项是Local和Cross term的总和。Local就是之前的说的定域项的贡献，Cross term是指这个基函数与其它所有基函数间交叉项中归属于这个基函数的贡献。
Multiwfn还输出了每个壳层对轨道的贡献，即壳层内基函数的贡献之和。注意除非分子的朝向已经清楚指定了，才可以在文中讨论某个基函数的贡献，否则讨论基函数的贡献其实是没有意义的，不少文献都犯了这个错误。因为分子在旋转后，各个基函数的贡献也会发生变化（除了S型以外，因为它是球对称的）。在不指明朝向时，只能对某个壳层的贡献进行讨论，因为Mulliken方法得出的壳层的贡献值是不会因为分子旋转而改变的。对于原子轨道的贡献来说也是一样，必须指明了朝向才能说某个原子轨道的贡献，否则只能说原子轨道壳层的贡献。注意在Stout-Politzer方法下，即便是壳层的贡献值也存在旋转依赖性。而SCPA方法下，必须使用的是球谐型高斯函数才能严格满足旋转不变性，不过使用笛卡尔型高斯函数时旋转依赖性并不大。这个问题在这里不做继续探讨，可参见《分子轨道成份的计算》的讨论。
由于4、8号基函数对应于氧的2py轨道，因此40.46960%+26.06437%=66.53%就是2py对第三条分子轨道的贡献，和前文用AOIM和NAO方法算出来的60%和73.2%都在比较相近。由于第三个和第五个壳层的贡献分别和第4号、第8号基函数贡献相同，也就是说在此MO当中氧的2px和2pz的成份都为0。
输出信息最后给出了通过基函数贡献加和得到的每个原子的贡献。
输入0退回上一级菜单。下面我们要进行片段分析。Multiwfn可以十分灵活、方便地自定义片段，片段不仅可以包含几个原子，也可以将特定的基函数、壳层、原子混合在一个片段里，片段的贡献值就是通过其中包含的基函数的贡献加和得到。这里我们先把氧的2p原子轨道壳层定义为第一个片段。选择-1进入第一个片段的定义界面，从all命令给出的列表以及6-31G*的基组定义中可以得知第三个和第五个壳层对应于氧的2p原子轨道壳层，于是就输入s 3,5。此时如果输入list命令来查看当前片段内已有的基函数列表，就会看到这两个壳层内包含的3,4,5,7,8,9号基函数都已经被加入了。输入q保存并退出片段定义界面。
然后选-2进入第二个片段定义界面。这里将两个氢原子（原子编号为2和3）定义为第二个片段，于是就输入a 2,3。输入q保存并退回。
选择选项4，Multiwfn就会使用Mulliken方法将第一个片段对所有轨道的贡献一次性都输出出来：
Orb# Type Energy Occ c^2 Int.cross Ext.cross Total
1 AB -20.558 2.000 0.000% 0.000% 0.000% 0.000%
2 AB -1.346 2.000 1.146% 0.474% 0.926% 2.546%
3 AB -0.714 2.000 35.134% 15.513% 15.887% 66.534%
4 AB -0.571 2.000 47.503% 22.651% 9.568% 79.722%
5 AB -0.498 2.000 67.063% 32.820% 0.000% 99.883%
6 AB 0.213 0.000 29.745% 10.680% -36.963% 3.462%
7 AB 0.307 0.000 81.284% 27.464% -105.805% 2.942%
8 AB 1.032 0.000 50.954% 6.230% -35.251% 21.933%
9 AB 1.133 0.000 79.565% -35.898% -3.568% 40.099%
10 AB 1.168 0.000 200.037% -100.056% 0.000% 99.980%
11 AB 1.178 0.000 63.463% -23.135% -4.472% 35.855%
12 AB 1.385 0.000 345.645% -160.782% -76.700% 108.163%
13 AB 1.431 0.000 159.364% -58.247% -66.716% 34.400%
14 AB 2.021 0.000 1.705% -0.037% -1.019% 0.650%
15 AB 2.031 0.000 0.000% 0.000% 0.000% 0.000%
16 AB 2.067 0.000 0.110% 0.027% 0.000% 0.136%
17 AB 2.636 0.000 55.219% -2.688% -50.347% 2.184%
18 AB 2.966 0.000 80.748% -0.574% -79.747% 0.427%
19 AB 3.978 0.000 11.297% -3.714% -6.502% 1.082%
片段1的贡献值Total是由c^2、Int.cross、Ext.cross三项加和得到的。c^2指的是片段内基函数的系数平方和。Int.cross（内交叉项）指的是片段1包含的基函数之间的交叉项之和。Ext.cross（外交叉项）指的是片段1的基函数与体系中其余基函数间的交叉项中被划归到片段1的部分。可以看到虚轨道当中出现了大量负值，这是Mulliken方法的弊病之一，因此如果需要分析虚轨道，建议用保证结果为正的SCPA，但此时将不会输出交叉项信息。轨道的Type显示的AB代表这个轨道可以同时占alpha电子和Beta电子，也即闭壳层轨道。
由于已经定义了片段2，所以下面还会输出片段1、2间的交叉项（如果只定义了片段1，则只输出上面那些内容）
Cross term between fragment 1 and 2 and their individual parts:
Orb# Type Energy Occ Frag1 part Frag2 part Total
1 AB -20.558 2.000 0.0000% 0.0000% 0.0000%
2 AB -1.346 2.000 0.9263% 0.9263% 1.8525%
3 AB -0.714 2.000 15.8871% 15.8871% 31.7743%
4 AB -0.571 2.000 9.5681% 9.5681% 19.1363%
5 AB -0.498 2.000 0.0000% 0.0000% 0.0000%
6 AB 0.213 0.000 -36.9633% -36.9633% -73.9266%
7 AB 0.307 0.000 -105.8051% -105.8051% -211.6102%
8 AB 1.032 0.000 -35.2511% -35.2511% -70.5022%
9 AB 1.133 0.000 -3.5678% -3.5678% -7.1357%
10 AB 1.168 0.000 0.0000% 0.0000% 0.0000%
11 AB 1.178 0.000 -4.4724% -4.4724% -8.9447%
12 AB 1.385 0.000 -76.6996% -76.6996% -153.3991%
13 AB 1.431 0.000 -66.7162% -66.7162% -133.4324%
14 AB 2.021 0.000 -1.0185% -1.0185% -2.0370%
15 AB 2.031 0.000 0.0000% 0.0000% 0.0000%
16 AB 2.067 0.000 0.0000% 0.0000% 0.0000%
17 AB 2.636 0.000 -50.3467% -50.3467% -100.6933%
18 AB 2.966 0.000 -79.7466% -79.7466% -159.4933%
19 AB 3.978 0.000 -6.5016% -6.5016% -13.0032%
由于Mulliken方法是将交叉项平分，所以片段1和片段2所占的交叉项是相同的。两个片段间总交叉项描述了两个片段间的相互作用程度，这与Mulliken键级直接相关。比如第三条MO，0.3177乘上轨道的占据数2得到的0.6354实际上就是此MO对片段1、2间的Mulliken键级的贡献。
6.2 在Multiwfn中使用Hirshfeld方法分析原子的贡献
下面介绍一下用Hirshfeld方法计算原子的贡献。生成各个原子的Hirshfeld权重时要用到原子在自由状态下的密度。Multiwfn提供了两种方法计算：(1)直接用程序内置的原子密度，详见手册附录3，这种做法对用户来说最方便 (2)基于原子波函数文件计算原子密度，详见手册3.7.3节或《使用Multiwfn作电子密度差图（http://sobereva.com/113）的第6~8节。两种做法结果虽然有异，但差别很小。为了简便，本文例子使用方法(1)。
在轨道成份分析功能里选8进入Hirshfeld分析，然后选1用程序内置的原子密度，之后Multiwfn会先进行初始化。然后比如要分析第3条分子轨道就输入3，将看到以下结果：
After normalization:
Atom 1(O ) : 64.926881%
Atom 2(H ) : 17.536559%
Atom 3(H ) : 17.536559%
由于数值积分精度原因，直接算出来的结果的加和值可能偏离100%一些，所以程序会自动做一下归一化。
如果你想获得某个原子在指定分子轨道范围内的贡献的话，就输入-2，之后输入原子编号和轨道范围即可，比如输入1和1,5就得到了氧原子在1至5号MO中的贡献：
Orb# Type Energy Occ Composition Population
1 Alpha&Beta -20.558 2.000 99.676867% 1.993537
2 Alpha&Beta -1.346 2.000 80.133155% 1.602663
3 Alpha&Beta -0.714 2.000 64.926646% 1.298533
4 Alpha&Beta -0.571 2.000 82.437155% 1.648743
5 Alpha&Beta -0.498 2.000 88.761277% 1.775226
Population of this atom in these orbitals: 8.318702
这里还给出了Population值，也就是原子对轨道的贡献值乘以轨道的占据数。最后给出了它在这些轨道的布居数加和值。由于1至5号就是全部的占据轨道，因此氧原子在水分子中的总电子数就是8.318702，也因此它的Hirshfeld电荷就是8-8.318702=-0.3187。
还可以选-9来定义片段，然后比如输入2,3就把两个氢都加进片段里了。之后输入轨道序号后，就会看到片段的贡献（Fragment contribution）。也可以选-3，然后输入轨道范围，来得到此片段对这些轨道各自的贡献。
6.3 在Multiwfn中使用Becke、Hirshfeld-I、AIM方法分析原子的贡献
在轨道成份分析功能里选9就能进入Becke分析。由于使用Becke方法和使用Hirshfeld方法分析原子的贡献的操作完全相同，因此这里不再举例。一般而言，这两种方法所得结果是比较接近的，用哪个都可以。对于大体系的话Hirshfeld方法更便宜，而且由于Hirshfeld方法的物理意义更强，我更倾向于用Hirshfeld方法。
使用Multiwfn基于Hirshfeld-I划分计算轨道成份在本文也不举例了，因为这种做法虽然比Hirshfeld方法更严格，但需要额外多花不少时间来通过迭代改进原子空间，而从结果来看也没显著优势。如果你对结果的物理意义要求特别高的话倒是可以考虑用这个。一般的用法是把Multiwfn自带的examples目录下的atmrad目录拷到当前目录下，然后进入轨道成分分析模块后，选择Hirshfeld-I方法，然后选相应选项开始计算即可。
Multiwfn还支持基于AIM划分计算轨道中原子的成份，由于不常用，这里也不给出例子了。这种做法计算的耗时很高，因为产生AIM盆的步骤非常昂贵，而结果并不比基于Hirshfeld或Becke方法算的好。如果确实想用这种做法，看Multiwfn手册4.8.6节。
6.4 在Multiwfn中使用NAO方法分析原子轨道、壳层、原子的贡献
如4.2节的例子所示，计算原子轨道的贡献就是将以NAO为基的分子轨道系数矩阵的相应矩阵元取平方乘100%。Multiwfn自身并不会生成这个矩阵，用户需要将含有这个矩阵的NBO程序的输出文件作为Multiwfn的输入文件。对于Gaussian，将route section写上pop=nboread，将分子坐标末尾空一行写上$NBO NAOMO $END，则这个Gaussian输出文件就能作为Multiwfn进行NAO分析的输入文件了（不要误用.fch文件作为输入文件！）。如果你用的是独立的NBO程序，即GENNBO，则应在.47文件的$NBO和$END中间添上NAOMO关键词，产生的输出文件也将可以作为Multiwfn的输入文件。
进入Multiwfn后选8，再选7，就进入了NAO分析模块。选0再输入分子轨道序号就能输出这个分子轨道的详细成份信息。输出哪些项可以通过2 Select output mode选项来控制，默认是只输出贡献大于0.5%的Cor和Val型的NAO和壳层。比如分析第3条轨道会输出：
NAO# Center Label Type Composition
7 1(O ) py Val( 2p) 73.154%
16 2(H ) S Val( 1S) 13.279%
18 3(H ) S Val( 1S) 13.279%
Condensed NAO terms to shells:
Atom: 1(O ) Shell: 5( 2p Val) 73.154%
Atom: 2(H ) Shell: 8( 1S Val) 13.279%
Atom: 3(H ) Shell: 10( 1S Val) 13.279%
Condensed NAO terms to atoms:
Center Composition
1(O ) 73.416%
2(H ) 13.288%
3(H ) 13.288%
Core composition: 0.000%
Valence composition: 99.711%
Rydberg composition: 0.280%
输出信息简单易懂。比如可见主要贡献来自于O的2p壳层，更具体来说是2py。当前轨道的Rydberg composition只有0.28%，比较小，说明对这个轨道做NAO轨道成份分析是有意义的。
在NAO分析模块中也同样可以定义片段来计算它对一批轨道的贡献，我们这里计算两个氢上的NAO对1至8号分子轨道上的贡献。因此输入
-1 //进入片段设定界面
a 2,3 //将2、3原子上的NAO都加入片段
q //保存片段
1 //计算片段对一批轨道的贡献
1-8 //要考察的轨道序号范围
结果如下所示
Orb.# Core Valence Rydberg Total
1 0.000% 0.119% 0.002% 0.121%
2 0.000% 18.556% 0.054% 18.611%
3 0.000% 26.557% 0.018% 26.576%
4 0.000% 7.076% 0.012% 7.089%
5 0.000% 0.000% 0.000% 0.000%
6 0.000% 35.482% 46.832% 82.314%
7 0.000% 27.543% 61.561% 89.104%
8 0.000% 32.433% 29.553% 61.986%
由于氢没有内核轨道，所以Core项全部为零。Valence、Rydberg项分别是指这个片段内所有Val型、Ryd型NAO的贡献。Total就是Core、Valence、Rydberg项的总和。前五个分子轨道是占据轨道，与之前讨论的一致，Ryd的成份非常低，这时可以说这两个氢对MO3的贡献是26.557%。而对于后面的虚轨道，则Ryd成份很大，不能忽略。由于Ryd型轨道物理意义不清，与原子轨道没有直接对应关系，难以讨论，所以此时就不建议再用NAO方法来分析虚轨道成份了，而更适合Hirshfeld方法分析原子的贡献。
7 总结
本文将轨道成份的计算原理、要注意的问题和具体操作过程都做了讨论。计算轨道成份不能稀里糊涂，一定要思路清楚。很多人急于算出轨道成份，在网上随便找一个来路不明的号称能计算轨道成份的软件就用于自己的计算，或者网上道听途说一个计算方法，比如什么“轨道成份就是系数的平方”，就按照这个方法算，这都是很危险的。
对于计算基函数的贡献，最建议用SCPA方法；对于计算原子轨道的贡献，最建议用NAO方法；对于计算原子的贡献，Hirshfeld、Hirshfeld-I、Becke方法以及将NAO方法得到的原子轨道的贡献进行加和的方法都推荐使用。注意NAO方法不适合分析虚轨道。Multiwfn也支持AIM方法计算原子的贡献，但由于此方法耗时高，所以一般不用。
加和SCPA方法得到的基函数的贡献来得到原子轨道的贡献，以及进一步加和得到原子轨道的贡献，虽然原理上不如上面建议的方法，但是只要没有弥散函数，且基组不平衡性较小时，且分析的不是较高能态虚轨道时，结果一般没什么问题。在使用SCPA方法分析原子轨道成份时最好（但不是必须）使用球谐型高斯函数以避免旋转依赖性，也避免笛卡尔型高斯函数引入的额外的基函数影响分析结果。
最后再次提醒一下，讨论某个基函数、某个原子轨道贡献时必须指明分子朝向，否则没有意义。不指定朝向时只能讨论壳层的贡献。