1.自然轨道的定义
对于一个单或多行列式波函数方法（例如RHF, MP2, CCSD, CASCI, CASSCF等等），可将电荷密度（charge density）
\rho(\boldsymbol{r})展开到一组正交归一的轨道
\psi_i(\boldsymbol{r})上
\rho(\boldsymbol{r})=\sum\limits_{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})=\sum\limits_{\mu}^{\rm nbf}C_{\mu i}\chi_{\mu}(\boldsymbol{r})S_{\mu\nu}=\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})=\sum\limits_{i}^{\rm nmo}U_{ik}\psi_{i}(\boldsymbol{r})\psi_{i}(\boldsymbol{r})=\sum\limits_{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})=\sum\limits_{ij}^{\rm nmo}D_{ij}\psi_i^*(\boldsymbol{r})\psi_j(\boldsymbol{r})=\sum\limits_{ij}^{\rm nmo}D_{ij}\sum\limits_{kl}^{\rm nmo}U_{ik}^{*}U_{jl}\tilde{\psi}_k^*(\boldsymbol{r})\tilde{\psi}_l(\boldsymbol{r})=\sum\limits_{kl}^{\rm nmo}(\sum\limits_{ij}^{\rm nmo}U_{ik}^{*}D_{ij}U_{jl})\tilde{\psi}_k^*(\boldsymbol{r})\tilde{\psi}_{l}(\boldsymbol{r})=\sum\limits_{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_i\delta_{ij}\delta_{ij}=
\left\{
\begin{array}{lr}
0, i\neq j \\
1, i=j
\end{array}
\right.
此即为自然轨道表象，对应的分子轨道称为自然轨道（natural orbital, NO），不妨专门记为
\psi_{i}^{\rm NO}(\boldsymbol{r})。此时
\rho(\boldsymbol{r})可简写为
\rho(\boldsymbol{r})=\sum\limits_{i}^{\rm nmo}n_i\psi_{i}^{\rm NO*}(\boldsymbol{r})\psi_{i}^{\rm NO}(\boldsymbol{r})=\sum\limits_{i}^{\rm nmo}n_i|\psi_{i}^{\rm NO}(\boldsymbol{r})|^2n_i即为
\tilde{D}_{ii}，称为第
i个自然轨道的占据数（natural orbital occupation number, NOON）。这式子很好理解：用电子数作为权重乘以轨道模方。所有轨道占据数加起来即为体系总电子数
\sum\limits_{i}^{\rm nmo}n_i=N
举个简单的例子，在RHF方法里
n_i取值只能是整数2/0，对应双占据/空轨道
n_i=\left\{
\begin{array}{lr}
2, i\in {\rm occupied} \\
0, i\in {\rm unoccupied}
\end{array}
\right.
电荷密度简化为
\rho(\boldsymbol{r})=2\sum\limits_{i}^{\rm occ}\psi_{i}^{*}(\boldsymbol{r})\psi_i(\boldsymbol{r})=2\sum\limits_{i}^{\rm occ}\sum\limits_{\mu \nu}^{\rm nbf}C_{\mu i}^{*}C_{\nu i}\chi_{\mu}^{*}(\boldsymbol{r})\chi_{\nu}(\boldsymbol{r})=\sum\limits_{\mu \nu}^{\rm nbf}P_{\mu \nu}\chi_{\mu}^{*}(\boldsymbol{r})\chi_{\nu}(\boldsymbol{r})P_{\mu \nu}=2\sum\limits_{i}^{\rm occ}C_{\mu i}^{*}C_{\nu i}
即求和指标从所有轨道（nmo）减小为双占据轨道（occ），接着将分子轨道展开至原子基函数（AO basis）上，便可出现大家在量子化学课上熟知的RHF（原子基）密度矩阵元
P_{\mu \nu}，写成矩阵形式
\rm P=2\boldsymbol{\rm CC}^{\dagger}
如果有一组轨道是当前轨道的酉变换，密度矩阵不会变，占据轨道的占据数也仍是2，即
\rm \boldsymbol{\rm \tilde{C}}=\boldsymbol{\rm CU}2\boldsymbol{\rm \tilde{C}\tilde{C}^{\dagger}}=2\boldsymbol{\rm CUU}^{\dagger}\boldsymbol{\rm C}^{\dagger}=2\boldsymbol{\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})=\sum\limits_{i}^{\rm occ}n_i\psi_i^*(\boldsymbol{r})\psi_i(\boldsymbol{r})=\sum\limits_{i}^{\rm occ}n_i\sum\limits_{\mu \nu}^{\rm nbf}C_{\mu i}^{*}C_{\nu i}\chi_{\mu}^{*}(\boldsymbol{r})\chi_{\nu}(\boldsymbol{r})=\sum\limits_{\mu \nu}^{\rm nbf}P_{\mu \nu}\chi_{\mu}^{*}(\boldsymbol{r})\chi_{\nu}(\boldsymbol{r})P_{\mu \nu}=\sum\limits_{i}^{\rm nmo}n_iC_{\mu i}^{*}C_{\nu i}
写成矩阵形式即为
\boldsymbol{\rm P}=\boldsymbol{\rm C\eta C}^{\dagger}
其中方阵
\boldsymbol{\rm \eta}的非对角元全是0，对角元
\eta_{ii}=n_i。在此一般形式的波函数下，若自然轨道经过酉变换
\boldsymbol{\tilde{\rm C}}=\boldsymbol{\rm CU}\boldsymbol{\rm P}=\boldsymbol{\rm C\eta 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})=\sum\limits_{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})=\sum\limits_{i}^{\rm nmo}n_i=N\qquad(0\leqslant x\leqslant1)
但是，仅
\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})=\sum\limits_{ij}^{\rm nmo}D_{ij}\tilde{\psi}_i^*(\boldsymbol{r})\tilde{\psi}_j(\boldsymbol{r})=\sum\limits_{ij}^{\rm nmo}D_{ij}\sum\limits_{\mu \nu}^{\rm nbf}\tilde{C}_{\mu i}^*\tilde{C}_{\nu j}\tilde{\chi}_{\mu}^*(\boldsymbol{r})\tilde{\chi}_{\nu}(\boldsymbol{r})=\sum\limits_{\mu \nu}^{\rm nbf}(\sum\limits_{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}=\sum\limits_{ij}^{\rm nmo}\tilde{C}_{\mu i}^*D_{ij}\tilde{C}_{\nu j} \boldsymbol{\rm P}=\boldsymbol{\rm \tilde{C}D\tilde{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语法）写成。本篇为理论篇，在不确定的将来某天还有实战篇，先看看反响如何?
