自然跃迁轨道分析
在《分析激发态的跃迁类型》一文中我们介绍了如何分析电子激发的跃迁性质。在TD-DFT框架下，通过程序输出的轨道跃迁系数即可计算出相应的跃迁贡献，进而可以根据占主要贡献的轨道的特征来分析跃迁属性。在某些体系中，可能无法用一两对轨道的跃迁来简单描述跃迁的属性，此时可以借助自然跃迁轨道（natural transition orbital, NTO）来进行分析。关于NTO的原理，可以参考其原始文献J. Chem. Phys. 2003, 118, 4775。在进行TD-DFT计算后，可以得到基态到某激发态之间的跃迁密度矩阵T，其维度为nocc×nvir。将T进行SVD分解，即
UTV† = Λ
得到两个酉矩阵U和V。将正则的占据轨道通过U矩阵进行变换，即得到一套占据NTO；将正则的虚轨道通过V矩阵进行变换，得到一套虚NTO。NTO的本征值的范围为[0,1]。绝大多数情况下，体系的虚轨道的数目远多于占据轨道数目。在虚NTO中有nocc个的本征值与占据NTO的本征值分别相同，它们一一匹配构成了NTO跃迁，剩余的虚NTO的本征值为0。假设NTO的最大本征值为0.9，则意味着该组跃迁占了90%的贡献。简言之，NTO是由正则分子轨道通过某种变换得到，并试图将电子激发描述为一对或少数几对轨道的跃迁。
很多量子化学程序都自带NTO分析。本文介绍如何使用Gaussian及Multiwfn进行NTO分析，具体程序版本为Gaussian 16 C.01和Multiwfn 3.8(dev)。示例体系为苯丙氨酸：
首先进行结构优化：
%chk=PhA-opt.chk
#p opt freq b3lyp/6-311G(d,p)
Title Card Required
0 1
N -1.87585202 0.43021225 0.00000000
H -0.87629202 0.40054625 0.00000000
C -2.56280602 1.70602425 0.00000000
H -3.18763902 1.78282725 -0.88982300
C -3.44942902 1.85214625 1.23214100
C -1.57232702 2.86163525 0.00000000
H -4.19366402 1.05582425 1.24137300
H -2.83665602 1.78649825 2.13119600
C -4.20885802 3.15424325 1.32108500
O -0.33540102 2.63015025 0.00000000
C -5.05434102 3.39852125 2.40988500
C -4.06748202 4.11720725 0.31474800
H -5.16441302 2.64878525 3.19339000
C -5.75844802 4.60576325 2.49235000
C -4.77158902 5.32444925 0.39721300
H -3.40921302 3.92701925 -0.53296100
H -6.41671702 4.79595125 3.34005900
C -5.61707202 5.56872725 1.48601300
H -4.66151702 6.07418525 -0.38629200
H -6.16526902 6.50865225 1.55021700
H -2.40132388 -0.42059873 0.00000000
O -2.04731713 4.21044376 0.00000000
H -1.36550372 4.79029233 -0.34713989
然后做TD计算：
%oldchk=PhA-opt.chk
%chk=PhA-TD.chk
#p td(nstates=10) b3lyp/6-311G(d,p) guess=read geom=allcheck
此时，我们可以看到第一激发态的信息如下：
Excited State 1: Singlet-A 5.3156 eV 233.25 nm f=0.0109 <S**2>=0.000
42-> 45 -0.27579
42-> 46 -0.20117
42-> 47 -0.14103
43-> 46 0.34319
44-> 45 0.46482
44-> 46 -0.14454计算可知，44→45的贡献为43.2%，43→46的贡献为23.6%，42→45的贡献为15.2%，没有占绝对优势贡献的跃迁。
下面进行NTO分析，基于已经完成的TD计算，输入文件如下：
%oldchk=PhA-TD.chk
%chk=PhA-NTO.chk
#p b3lyp chkbasis guess=(only,read) geom=allcheck
pop=SaveNTO density=(check,transition=1)其中，density=transition=1表示使用第一个激发态的跃迁密度矩阵，pop=saveNTO表示将NTO轨道保存到chk文件中，便于使用GaussView或其他程序观看轨道。
在输出文件中，会按本征值从小到大和从大到小分别列出占据NTO和虚NTO的本征值。本例中部分轨道的本征值如下：
Alpha occ. eigenvalues -- 0.00005 0.00005 0.00005 0.00007 0.00010
Alpha occ. eigenvalues -- 0.00012 0.00014 0.00018 0.00021 0.00025
Alpha occ. eigenvalues -- 0.00042 0.04340 0.35180 0.60392
Alpha virt. eigenvalues -- 0.60392 0.35180 0.04340 0.00042 0.00025
Alpha virt. eigenvalues -- 0.00021 0.00018 0.00014 0.00012 0.00010
Alpha virt. eigenvalues -- 0.00007 0.00005 0.00005 0.00005 0.00004可以看到从最高占据跃迁轨道（highest occupied transition orbital, HOTO）到最低未占据跃迁轨道（lowest unoccupied transition orbital, LUTO）的贡献占了60.4%，而HOTO−1→LUTO+1的贡献占了35.2%，其他贡献可忽略不计。将chk文件转换为fchk文件后用GaussView打开，轨道能量排序也是按照上述规则：
绘制出4个轨道的等值面（isovalue=0.03）如下所示：
因此该电子激发可指认为局域π→π*跃迁。
本例中，我们先进行了普通的TD-DFT计算，再做NTO分析。也可以在TD计算时，直接进行NTO分析，并将NTO信息保存下来，其关键词如下：
#p td(nstates=10) b3lyp/6-311G(d,p) guess=read geom=allcheck
pop=(NTO,SaveNTO) density(transition=1) 使用Gaussian进行NTO分析的一个不便之处是每次计算只能分析一个激发态，也就是说每次需要通过指定的density=transition=n来进行计算。而使用Multiwfn则可以方便地分析多个激发态，具体可参考Multiwfn的手册以及卢天老师的教程：
http://sobereva.com/377
需要注意的是，使用Multiwfn分析需要更精细的跃迁系数，因此在TD计算时要加入iop(9/40=4)，以输出绝对值大于10−4的系数。此处也展示一下分析流程。
首先载入PhA-TD.fchk文件，选择18-电子激发分析，接着选择6-生成自然跃迁轨道，此时会提示需要提供输出文件的路径，如果输出文件与fchk文件同名且在同一目录下，直接按回车即可。接着输入1，表示分析第一个激发态，此时本征值最大的10个占据NTO和虚NTO的本征值会输出在屏幕上，排序方式与Gaussian相同：
The highest 10 eigenvalues of occupied NTOs:
0.000099 0.000122 0.000143 0.000176 0.000207
0.000249 0.000417 0.043400 0.351797 0.603922
The highest 10 eigenvalues of virtual NTOs:
0.603922 0.351797 0.043400 0.000417 0.000249
0.000207 0.000176 0.000143 0.000122 0.000099下面还提示可将NTO轨道信息写入molden、fchk或mwfn文件，此处我们选择2，写入fch文件，提出保存路径即可。假设我们保存为S1.fch，重新打开Multiwfn并载入之，则可绘制出轨道：
