1.自然轨道的定义
对于一个单或多行列式波函数方法（例如RHF, MP2, CCSD, CASCI, CASSCF等等），可将电荷密度（charge density）
rho(boldsymbol{r})展开到一组正交归一的轨道
psi_i(boldsymbol{r})上
rho(boldsymbol{r})=sumlimits_{ij}^{rm nmo}D_{ij}psi_i^*(boldsymbol{r})psi_j(boldsymbol{r})
注意电荷密度是个全空间的函数。其中nmo表示分子轨道数，
boldsymbol{rm D}是个厄米矩阵（实数下即实对称阵），此处表示分子轨道基下的密度矩阵，矩阵的迹必须等于体系电子数
{rm trace}(boldsymbol{rm D})=N
这里所说的
psi_i(boldsymbol{r})满足正交归一关系是指
psi_i(boldsymbol{r})=sumlimits_{mu}^{rm nbf}C_{mu i}chi_{mu}(boldsymbol{r})S_{munu}=int chi_{mu}^{*}(boldsymbol{r})chi_{nu}(boldsymbol{r}){rm d}boldsymbol{r}boldsymbol{rm C}^{dagger}boldsymbol{rm SC}=boldsymbol{rm I}
其中
chi_{mu}(boldsymbol{r})为原子基函数，nbf为基函数数目，
boldsymbol{rm S}为原子基函数之间的重叠积分。
对于一个给定的体系和确定的波函数方法，
rho(boldsymbol{r})是固定的，因此若换成另一组正交归一的轨道 ，便会对应一个新的矩阵
boldsymbol{rm tilde{D}}，写成公式就是
tilde{psi}_{k}(boldsymbol{r})=sumlimits_{i}^{rm nmo}U_{ik}psi_{i}(boldsymbol{r})psi_{i}(boldsymbol{r})=sumlimits_{k}^{rm nmo}U_{ik}tilde{psi}_{k}(boldsymbol{r})boldsymbol{rm UU}^{dagger}=boldsymbol{rm U}^{dagger}boldsymbol{rm U}=boldsymbol{rm I}rho(boldsymbol{r})=sumlimits_{ij}^{rm nmo}D_{ij}psi_i^*(boldsymbol{r})psi_j(boldsymbol{r})=sumlimits_{ij}^{rm nmo}D_{ij}sumlimits_{kl}^{rm nmo}U_{ik}^{*}U_{jl}tilde{psi}_k^*(boldsymbol{r})tilde{psi}_l(boldsymbol{r})=sumlimits_{kl}^{rm nmo}(sumlimits_{ij}^{rm nmo}U_{ik}^{*}D_{ij}U_{jl})tilde{psi}_k^*(boldsymbol{r})tilde{psi}_{l}(boldsymbol{r})=sumlimits_{kl}^{rm nmo}tilde{D}_{kl}tilde{psi}_k^*(boldsymbol{r})tilde{psi}_{l}(boldsymbol{r})
其中
tilde{boldsymbol{rm D}}=boldsymbol{rm U}^{dagger}boldsymbol{rm DU}
这其中有个特殊的酉变换尤其重要：存在一个特殊的
boldsymbol{rm U}使得
tilde{boldsymbol{rm D}}是对角矩阵，即
tilde{D}_{ij}=n_idelta_{ij}delta_{ij}=
left{
begin{array}{lr}
0, ineq j \
1, i=j
end{array}
right.
此即为自然轨道表象，对应的分子轨道称为自然轨道（natural orbital, NO），不妨专门记为
psi_{i}^{rm NO}(boldsymbol{r})。此时
rho(boldsymbol{r})可简写为
rho(boldsymbol{r})=sumlimits_{i}^{rm nmo}n_ipsi_{i}^{rm NO*}(boldsymbol{r})psi_{i}^{rm NO}(boldsymbol{r})=sumlimits_{i}^{rm nmo}n_i|psi_{i}^{rm NO}(boldsymbol{r})|^2n_i即为
tilde{D}_{ii}，称为第
i个自然轨道的占据数（natural orbital occupation number, NOON）。这式子很好理解：用电子数作为权重乘以轨道模方。所有轨道占据数加起来即为体系总电子数
sumlimits_{i}^{rm nmo}n_i=N
举个简单的例子，在RHF方法里
n_i取值只能是整数2/0，对应双占据/空轨道
n_i=left{
begin{array}{lr}
2, iin {rm occupied} \
0, iin {rm unoccupied}
end{array}
right.
电荷密度简化为
rho(boldsymbol{r})=2sumlimits_{i}^{rm occ}psi_{i}^{*}(boldsymbol{r})psi_i(boldsymbol{r})=2sumlimits_{i}^{rm occ}sumlimits_{mu nu}^{rm nbf}C_{mu i}^{*}C_{nu i}chi_{mu}^{*}(boldsymbol{r})chi_{nu}(boldsymbol{r})=sumlimits_{mu nu}^{rm nbf}P_{mu nu}chi_{mu}^{*}(boldsymbol{r})chi_{nu}(boldsymbol{r})P_{mu nu}=2sumlimits_{i}^{rm occ}C_{mu i}^{*}C_{nu i}
即求和指标从所有轨道（nmo）减小为双占据轨道（occ），接着将分子轨道展开至原子基函数（AO basis）上，便可出现大家在量子化学课上熟知的RHF（原子基）密度矩阵元
P_{mu nu}，写成矩阵形式
rm P=2boldsymbol{rm CC}^{dagger}
如果有一组轨道是当前轨道的酉变换，密度矩阵不会变，占据轨道的占据数也仍是2，即
rm boldsymbol{rm tilde{C}}=boldsymbol{rm CU}2boldsymbol{rm tilde{C}tilde{C}^{dagger}}=2boldsymbol{rm CUU}^{dagger}boldsymbol{rm C}^{dagger}=2boldsymbol{rm CC}^{dagger}=boldsymbol{rm P}
因此RHF正则占据轨道（及其任意酉变换）本身也是自然轨道。
对于UHF则有四种常见的密度矩阵：alpha自旋，beta自旋，自旋密度矩阵（即alpha-beta密度矩阵差），总密度（即alpha beta密度矩阵和），对应四种自然轨道：alpha自然轨道，beta自然轨道，自旋自然轨道(SNO)，UHF自然轨道（UNO）。前两种的轨道占据数严格为1，因此alpha/beta正则占据轨道（及其任意酉变换）本身亦是自然轨道；后两种则需要对角化相应的密度矩阵得到（见下文），轨道占据数是小数，SNO占据数范围[-1,1]，UNO占据数范围[0,2]。
在一般的方法里，例如MP2, CCSD, CAS等波函数含有多个行列式，其自然轨道没有严格的占据与空的概念，占据数一般也不是整数，也是[0,2]范围的小数。注意像MP2和CCSD的参考态RHF是单行列式的，但（请按/符号断句）/MP2里的一阶波函数/和/CCSD波函数/是多行列式的；而CASCI和CASSCF波函数本身就是多行列式的。
对一般的波函数而言，将电荷密度展开至原子基函数上（简便起见，省略上标NO）
rho(boldsymbol{r})=sumlimits_{i}^{rm occ}n_ipsi_i^*(boldsymbol{r})psi_i(boldsymbol{r})=sumlimits_{i}^{rm occ}n_isumlimits_{mu nu}^{rm nbf}C_{mu i}^{*}C_{nu i}chi_{mu}^{*}(boldsymbol{r})chi_{nu}(boldsymbol{r})=sumlimits_{mu nu}^{rm nbf}P_{mu nu}chi_{mu}^{*}(boldsymbol{r})chi_{nu}(boldsymbol{r})P_{mu nu}=sumlimits_{i}^{rm nmo}n_iC_{mu i}^{*}C_{nu i}
写成矩阵形式即为
boldsymbol{rm P}=boldsymbol{rm Ceta C}^{dagger}
其中方阵
boldsymbol{rm eta}的非对角元全是0，对角元
eta_{ii}=n_i。在此一般形式的波函数下，若自然轨道经过酉变换
boldsymbol{tilde{rm C}}=boldsymbol{rm CU}boldsymbol{rm P}=boldsymbol{rm Ceta C}^{dagger}=boldsymbol{rm C}(boldsymbol{rm UU}^{dagger})boldsymbol{eta}(boldsymbol{rm UU}^{dagger})boldsymbol{rm C}^{dagger}=(boldsymbol{rm CU})(boldsymbol{rm U}^{dagger}boldsymbol{rm eta}boldsymbol{rm U})(boldsymbol{rm CU})^{dagger}=boldsymbol{tilde{rm C}tilde{rm eta}tilde{rm C}}^{dagger}boldsymbol{rm tilde{eta}}=boldsymbol{rm U}^{dagger}boldsymbol{rm eta}boldsymbol{rm U}
中间的矩阵
boldsymbol{rm tilde{eta}}不再是对角矩阵，只是对称矩阵。因此对于一般形式的波函数，自然轨道是唯一的，正则轨道或其他类型的轨道不能充当自然轨道的角色。假设我们知道原子基密度矩阵
boldsymbol{rm P}（例如Gaussian的fchk文件里就有），如何求自然轨道占据数和自然轨道呢？
2.从密度矩阵求自然轨道
直接对角化矩阵
boldsymbol{rm P}是不行的，因为（1）自然轨道
boldsymbol{rm C}不是酉矩阵；（2）没法保证矩阵
boldsymbol{rm P}本征值的和等于总电子数
N。注意到我们有正交归一关系
boldsymbol{rm C}^{dagger}boldsymbol{rm SC}=boldsymbol{rm I}，我们可以给矩阵
boldsymbol{rm P}左右各乘一个
boldsymbol{rm S}^{1/2}boldsymbol{rm S}^{1/2}boldsymbol{rm PS}^{1/2}=boldsymbol{rm S}^{1/2}boldsymbol{rm C eta C}^{dagger}boldsymbol{rm S}^{1/2}=(boldsymbol{rm S}^{1/2}boldsymbol{rm C})boldsymbol{rm eta}(boldsymbol{rm C}^{dagger}boldsymbol{rm S}^{1/2})
关于
boldsymbol{rm S}^{1/2}可阅读公众号本期另一文《
boldsymbol{rm S}^{1/2}的一些性质》。此矩阵的迹
{rm trace}(boldsymbol{rm S}^{1/2}boldsymbol{rm PS}^{1/2})={rm trace}(boldsymbol{rm S}^{1/2}boldsymbol{rm C}(boldsymbol{rm eta C}^{dagger}boldsymbol{rm S}^{1/2}))={rm trace}((boldsymbol{rm eta C}^{dagger}boldsymbol{rm S}^{1/2})boldsymbol{rm S}^{1/2}boldsymbol{rm C})={rm trace}(boldsymbol{rm eta C}^{dagger}boldsymbol{rm S}^{1/2}boldsymbol{rm S}^{1/2}boldsymbol{rm C})={rm trace}(boldsymbol{rm eta C}^{dagger}boldsymbol{rm SC})={rm trace}(boldsymbol{rm eta})=sumlimits_{i}^{rm nmo}n_i=N
便是总电子数，符合要求。可能有读者会有疑问，非得乘
boldsymbol{rm S}^{1/2}？事实上，仅考虑电子数的话，不止一种方式，如下通式均满足要求
{rm trace}(boldsymbol{rm S}^{1-x}boldsymbol{rm P}boldsymbol{rm S}^{x})=sumlimits_{i}^{rm nmo}n_i=Nqquad(0leqslant xleqslant1)
但是，仅
boldsymbol{rm S}^{1/2}boldsymbol{rm PS}^{1/2}是对称矩阵。接着事情就很简单了，我们可以将这个对称矩阵对角化，
boldsymbol{rm S}^{1/2}boldsymbol{rm PS}^{1/2}=boldsymbol{rm U eta U}^{dagger}
对比上述刚乘
boldsymbol{rm S}^{1/2}时的形式可以发现
boldsymbol{rm U}=boldsymbol{rm S}^{1/2}boldsymbol{rm C}
则自然轨道系数矩阵为
boldsymbol{rm C}=boldsymbol{rm S}^{-1/2}boldsymbol{rm U}
在实际编程中求
boldsymbol{rm S}^{-1/2}时需要舍弃接近零的值，即处理线性依赖。相应地，本征值得自己从大到小排序（MKL库函数输出是从小到大），取到自然分子轨道数目即止。若有本征值被舍弃，则
boldsymbol{rm U}的对应本征矢也应该舍弃，保证最后自然轨道系数矩阵的维度是基函数*自然轨道数。
注意，自然轨道数小于等于总轨道数。例如在CASCI和CASSCF方法中，若提供的密度矩阵是活性空间密度矩阵，则求出来的自然轨道数只能等于活性轨道数。若提供的密度矩阵是总密度矩阵，则自然轨道数等于总轨道数。
我们已经知道如何求自然轨道及其占据数，接着回到原有的情形：假设有一套普通的正交归一轨道
boldsymbol{rm tilde{C}}（例如，Boys局域轨道，UNO轨道等等），它是自然轨道
boldsymbol{rm C}的酉变换
boldsymbol{rm tilde{C}}=boldsymbol{rm CU}，则对应的占据数矩阵为
boldsymbol{rm tilde{eta}}=boldsymbol{rm U}^{dagger}boldsymbol{rm eta U}
可见，只需先求出变换关系
boldsymbol{rm U}，（可以调MKL库中解线性方程组的函数），然后做两次矩阵乘法即可得到
boldsymbol{rm tilde{eta}}。由于该轨道不是自然轨道，没有严格的占据数概念，其占据数矩阵不是对角的，在非对角元上也有值。我们可以将
boldsymbol{rm tilde{eta}}的对角元
tilde{eta}_{ii}“近似地看成”第
i个轨道的“占据数”。假设
boldsymbol{rm tilde{eta}}的非对角元较小、对角元接近本征值，便可认为这套轨道与自然轨道较为接近，可以作为一种衡量接近自然轨道程度的指标。
最后，回到本文一开始的公式，假设我们现在将电荷密度展开在这组普通的正交归一轨道
boldsymbol{rm tilde{C}}上
rho(boldsymbol{r})=sumlimits_{ij}^{rm nmo}D_{ij}tilde{psi}_i^*(boldsymbol{r})tilde{psi}_j(boldsymbol{r})=sumlimits_{ij}^{rm nmo}D_{ij}sumlimits_{mu nu}^{rm nbf}tilde{C}_{mu i}^*tilde{C}_{nu j}tilde{chi}_{mu}^*(boldsymbol{r})tilde{chi}_{nu}(boldsymbol{r})=sumlimits_{mu nu}^{rm nbf}(sumlimits_{ij}^{rm nmo}tilde{C}_{mu i}^*D_{ij}tilde{C}_{nu j})tilde{chi}_{mu}^*(boldsymbol{r})tilde{chi}_{nu}(boldsymbol{r})P_{mu nu}=sumlimits_{ij}^{rm nmo}tilde{C}_{mu i}^*D_{ij}tilde{C}_{nu j} boldsymbol{rm P}=boldsymbol{rm tilde{C}Dtilde{C}}^{dagger}
对比上文的
boldsymbol{rm P}=boldsymbol{tilde{rm C}tilde{rm eta}tilde{rm C}}^{dagger}，可以发现
boldsymbol{tilde{rm eta}}就是
boldsymbol{rm D}。
Acknowledgement
感谢清癯、zhigang、Acid、yuqiwang、暖云大师和食肉动物等人的审阅和建议。
References
- Modern Quantum Chemistry, A. Szabo and N. S. Ostlund, p139.
- 《在Multiwfn中基于fch产生自然轨道的方法与激发态波函数、自旋自然轨道分析实例》http://sobereva.com/403
- http://gaussian.com/population
后记
本文使用mdnice的Chrome浏览器插件（支持Markdown和LaTex语法）写成。本篇为理论篇，在不确定的将来某天还有实战篇，先看看反响如何?
