从密度矩阵产生自然轨道_实战篇（上）
本公众号之前发过自然轨道的原理介绍，详见《从密度矩阵产生自然轨道-理论篇》和《S^(1/2)的一些性质》。对其中原理和公式熟悉的读者，可以自己编写代码从密度矩阵产生自然轨道，无需阅读本文。笔者也写过这部分代码，可以提供一个参考
https://gitlab.com/jxzou/mokit/-/blob/master/src/lo.f90见其中第一个subroutine no。
其实基本上每个量化程序自己都能产生自然轨道，不过一般在各软件手册里提及较少，导致有些人不知道、或者以为需要借助外部程序才能实现。本文就介绍一下如何在常见量化程序中获得自然轨道，以Gaussian软件为主，其他软件顺带提一下。RDFT/UDFT的操作与RHF/UHF的操作一样，不做单独说明。本文为上篇，分为4小节，分别是：
一、RHF自然轨道
二、UHF自然轨道
三、MP2和CCSD自然轨道
四、UMP2和UCCSD自然轨道
只对某一部分感兴趣的读者可以直接下划。更多内容请见下篇。
1
RHF自然轨道
在之前的原理介绍篇里讲过，RHF占据轨道的任意酉变化均是其自然轨道，即自然轨道有无数组，这样的轨道没有什么实际意义。
2
UHF自然轨道
UHF自然轨道简称为UNO。以单重态O2为例，Gaussian输入文件示例如下
%chk=O2_uhf.chk
#p UHF/def2TZVP nosymm guess=mix stable=opt
sp of singlet O2 [标题行，随便写]
0 1
O 0.0 0.0 0.0
O 0.0 0.0 1.1616
--Link1--
%chk=O2_uhf.chk
#p chkbasis nosymm guess(read,only,save,NaturalOrbitals) geom=allcheck可以看出，其基本逻辑是先做一个UHF计算，然后读取波函数（实际上是读取密度）产生自然轨道。chkbasis表示读取基组（若有赝势也可自动读），geom=allcheck表示读取电荷、自旋多重度和坐标，only表示不做能量计算，其他关键词顾名思义即可。获得自然轨道这一步的耗时一般只需几秒。另外，不习惯写--Link1--的读者可以拆成两个gjf文件来写。
注意，复杂体系的UHF/UDFT计算可能会得到内部不稳定（internal instability）的波函数，这些波函数做布局分析和产生自然轨道没有意义。因此这里先确保波函数是稳定的（即用stable=opt），然后才用来产生自然轨道。另外，对于过渡金属体系可能存在多个稳定的UHF解，能量各不一致，一般取能量最低的解来产生自然轨道。
对于单重态UHF计算，还应使用guess=mix来产生对称性破缺的初始猜测，对此不熟悉的小伙伴请阅读公众号发过的介绍《谈谈Gaussian软件中的guess=mix》。注意算完UHF后应先查看输出文件中自旋平方的期望值<S**2>，若其十分接近或等于零，此时UHF解实际上是RHF波函数，UHF能量等于RHF能量，即体系是闭壳层的，此时看自然轨道没有什么意义。
在计算完成后，用高斯自带的formchk小程序将chk文件转化为fchk文件，然后用GaussView打开fchk文件查看自然轨道（见下图）
（GV 6.0主面板上Tools->MOs->Visualize->选中或输入轨道序号->Update。注意不要点击Calculation。若Visualize是灰色的，说明电脑上没装高斯）
注意电子箭头旁的数据不是轨道能量，而是轨道占据数，为0~2之间的小数。若算出来有很明显的负数，说明自己很可能操作有误。自然轨道没有HOMO、LUMO这种说法，若非要一个称呼，可以称为HONO（最高占据自然轨道）、LUNO（最低未占据自然轨道）。从图中可以看出HONO、LUNO占据数都是1.0，是简并的，而其他自然轨道的占据数极其接近2/0，相当于RHF轨道。占据数情况表明，这是一个几乎完美的双自由基，双自由基特征y=1.0。关于双自由基我们以后再出一期详细的介绍。
当然，可视化轨道的软件远不止GaussView一个，读者可以用自己偏好的软件看轨道。高斯在UNO这个功能上有个缺憾（某种程度上可以说是个bug）：原理上UNO轨道只有1列，不区分alpha与beta，但由于前一步是UHF任务（有两列轨道），高斯没能及时把chk文件里两列的信息删减为1列，导致产生的fchk文件里有两列一模一样的轨道，读者随便看1列即可。
其他软件就不会有这个问题，例如在OpenMolcas软件中，UHF计算完后默认就会产生UNO，存于后缀.UnaOrb的文件中，且只有1列轨道，占据数可以看#OCHR底下的数据。在GAMESS里是在$SCF中加上UHFNOS=.true.，产生的自然轨道及其占据数在.dat文件中的UHF轨道后。而在ORCA中则是在%scf里写上UHFNO true，结果在O2_uhf.uno文件中，这实际上也是一个gbw文件，你可以重命名为.gbw后缀后再通过ORCA自带的orca_2mkl小程序生成mkl和molden文件。
各个软件的关键词简要总结
有小伙伴可能听说过、或阅读本文后察觉到，UNO轨道可以作为CASSCF计算的初始轨道。由于有占据数大小排序（明显偏离2/0的可以认为是活性轨道），还可以根据轨道形状来进一步挑选，可以认为是一种不错的初始。更为详细的介绍留待以后的CASSCF计算技巧系列篇。
3
MP2和CCSD自然轨道
此处特指基于RHF的MP2和CCSD方法，基于UHF的请阅读下一部分。以基态平衡位置附近的N2分子为例，Gaussian输入文件示例如下
%chk=N2_cc-pVTZ.chk
#p MP2/cc-pVTZ density
title
0 1
N 0.0 0.0 0.0
N 0.0 0.0 1.08606
--Link1--
%chk=N2_cc-pVTZ.chk
#p chkbasis guess(read,only,save,NaturalOrbitals) geom=allcheck唯一要注意的地方就是post-HF方法要加上density关键词，才会计算当前方法的密度，才能做布局分析和产生自然轨道，否则仍是HF下的结果。算CCSD的话就把MP2换成CCSD即可。高斯不支持CCSD(T)方法的密度，因此高斯无法产生CCSD(T)自然轨道。
对于OpenMolcas软件，RMP2计算默认不产生自然轨道，需要在&MBPT2中写上Prpt关键词（类似于高斯的density），生成的自然轨道存于后缀.MP2Orb的文件中。在GAMESS里则要在
各个软件的关键词简要总结
若无特别强调，表中均为弛豫（relaxed）密度。对于不优化轨道的方法（如MP2, CCSD, TDDFT等等），密度有弛豫与非弛豫密度（unrelaxed）之分，对应的自然轨道也略有不同。计算性质（如偶极矩、极化率等）推荐用弛豫密度，但其对应的自然轨道占据数理论上允许超出[0,2]的限制，可以说是个缺点。非弛豫密度则满足占据数[0,2]的限制。大部分程序只支持其中一种密度。这部分更详细的讨论请阅读
http://classic.chem.msu.su/cgi-bin/ceilidh.exe/gran/gamess/forum/?C34df668afbHW-7216-1405+00.htm和ORCA手册（以4.2.1版为例）中8.1.3.2节Coupled-Cluster Densities。
另外，开源免费的PySCF程序也支持MP2和CCSD非弛豫密度和对应的自然轨道，不过它的Python语句其实是自己写公式算，而非简单的一两个关键词。这里我们举一例输入文件
from pyscf import gto, scf, mp
import numpy as np
mol = gto.Mole()
mol.build(
atom = 'C 0 0 0; O 0 0 1.1281', # in Angstrom
basis = 'ccpvtz',
verbose = 4)
mf = scf.RHF(mol).run() # 算RHF能量
mypt2 = mp.RMP2(mf)
mypt2.frozen = 2
mypt2.kernel() # 算MP2能量
rdm1 = mypt2.make_rdm1() # 一阶约化密度矩阵
natocc, natorb = np.linalg.eigh(rdm1)
natorb = np.dot(mf.mo_coeff, natorb)
print(natocc) # 打印MP2自然轨道占据数4
UMP2和UCCSD自然轨道
仍以单重态O2为例，Gaussian输入文件示例如下
%chk=O2_uccsd.chk
#p UHF/def2TZVP nosymm guess=mix stable=opt
sp of singlet O2 [标题行，随便写]
0 1
O 0.0 0.0 0.0
O 0.0 0.0 1.1616
--Link1--
%chk=O2_uccsd.chk
#p UCCSD/chkbasis nosymm guess=read geom=allcheck density
--Link1--
%chk=O2_uccsd.chk
#p chkbasis nosymm guess(read,only,save,NaturalOrbitals) geom=allcheck笔者特意从UHF开始写，意在提醒读者，应先检验UHF波函数稳定性，确保稳定了才能算UMP2或UCCSD。看过《谈谈Gaussian软件中的guess=mix》一文的读者肯定知道，对于单重态直接仅写个UMP2是没有意义的，算出来必然与RMP2一模一样，白费时间。
自然键轨道（NBO）不是自然轨道，不在本文介绍范围内。下期预告：
五、CASSCF自然轨道
六、CASCI自然轨道
七、自然跃迁轨道
八、激发态自然轨道