理论计算干货：利用常见的程序做轨道局域化

本文主要介绍一下如何利用常见的程序做轨道局域化（也称定域化），为后续多参考态计算系列篇做点铺垫。至于如何自己写轨道局域化的代码，这又是另一个话题，以后再介绍。

1. Gaussian

高斯做 Boys 局域化十分简单，例如对苯分子的轨道局域化：

%chk=ben_6-31Gd.chk%mem=2GB%nprocshared=4#p RHF/6-31G(d)
Title Card Required
0 1C     1.05889731   -0.35087719    0.00000000C     2.45405731   -0.35087719    0.00000000C     3.15159531    0.85687381    0.00000000C     2.45394131    2.06538281   -0.00119900C     1.05911631    2.06530481   -0.00167800C     0.36151531    0.85709881   -0.00068200H     0.50913831   -1.30319419    0.00045000H     3.00356531   -1.30339019    0.00131500H     4.25127531    0.85695381    0.00063400H     3.00414131    3.01752581   -0.00125800H     0.50899431    3.01758581   -0.00263100H    -0.73808869    0.85728181   -0.00086200
--Link1--%chk=ben_6-31Gd.chk%mem=2GB%nprocshared=4#p RHF/6-31G(d) guess=(read,local,only,save) geom=allcheck

即只需在分步任务 Link1 中写 guess=(read,local,only,save)。read表示读取上一步收敛的HF波函数；local 表示进行 Boys 局域化；only表示不做 SCF 能量计算，此处即做完局域化后就退出程序；save 表示将局域化后的轨道存入 .chk 文件中。四个词缺一不可。注意如果 HF 计算中写了 nosymm 关键词，Link1 中最好也写上 nosymm。geom=allcheck 则是读取上一步的 Title Card 行、电荷、自旋多重度和坐标。高斯默认分别在占据轨道空间与空轨道空间内部做局域化。

局域化完成之后用 formchk 把 .chk 转化为 .fchk 文件，用GaussView 打开，绘制轨道（读者也可以用 Multiwfn 等软件看、绘制轨道）。在占据轨道 1-21 号中可以很容易找到 6 根局域的 C−H σ键、3 根 C−C σ 键，且在空轨道 22-102 号中可以很容易找到对应的反键轨道。然而另外 3 根 C−C σ 键和 3 个 C−C π 键却找不到，因为同一处空间位置上的 C−C σ 键和 π 键混合在一起了。无法区分哪个是 σ 键，哪个是 π 键，这便是 Boys 局域化方法的一个著名现象：σ-π 混合（不同角度讲既可以看作缺点也可以看作优点）。

▲ Boys局域轨道举例

如果想要不混合的局域轨道，得用 Pipek-Mezey  (简称 PM )局域化，只需在上面输入文件再续一个任务

--Link1--%chk=ben_6-31Gd.chk%mem=2GB%nprocshared=4#p RHF/6-31G(d,p) guess=(read,only,save) geom=allcheck iop(4/9=20212)

就可以看到清楚的 π 键了。需要注意的是，轨道局域化后已经没有轨道能级的概念了，当然也没有 HOMO、LUMO 这些说法了。不过仍然可以依据分子轨道 Fock 矩阵对角元值的大小来将局域轨道排序，这一点留待以后讲局域化代码时再介绍。

▲ PM 局域轨道举例（其余轨道如C-H键与Boys局域轨道无异）

读者可能会有疑问，为何这个 PM 局域化要接在 Boys 局域化后做。这与高斯这个功能写得比较差劲有关，比如苯的 Boys 局域化，基组从 6-31G(d) 换成 6-31G(d,p) 或 6-311G(d) 就会收敛失败；若用PM局域化，6-31G 或 6-31G(d)都会收敛失败。而且这种收敛失败并不会在 log 文件末尾报错体现，而是末尾正常结束 Normal termination，仔细往上翻找可以找到 Localization failed after …，打开 .fchk 文件看轨道会发现一团乱麻，所以现今几乎没人会选择用高斯做局域化。此处笔者仅为展示 PM 局域化关键词如何写，及轨道形状的 σ-π 分离现象，只好想方设法绕着弯在高斯里做出 PM 局域轨道。如若使用别的程序大可不必如此，可直接在 HF 计算后直接做 PM 局域化。

2. PySCF

PySCF 的输入文件示范如下

from pyscf import gto, scf, lofrom pyscf.tools import molden





    
mol = gto.M()mol.atom = '''[原子坐标]'''
mol.basis = '6-31G(d)'mol.max_memory = 2000 # MBmol.verbose = 4mol.build()
#RHF calculationmf = scf.RHF(mol)mf.kernel()
#localize the occ/unocc separatelyocc_idx = range(0,21)vir_idx = range(21,96)occ_loc_orb = lo.Boys(mol, mf.mo_coeff[:,occ_idx]).kernel()mf.mo_coeff[:,occ_idx] = occ_loc_orbvir_loc_orb = lo.Boys(mol, mf.mo_coeff[:,vir_idx]).kernel()mf.mo_coeff[:,vir_idx] = vir_loc_orb#localization done
#output Boys LMOsmolden.from_mo(mol, 'ben_6-31Gd_boys.molden', mf.mo_coeff)

注意 # 符号后都是 python 的注释内容，可以省去。 PySCF 支持对指定的某些轨道做局域化，十分灵活，更多实例可去 pyscf/examples/local_orb 目录下找。 此处笔者写的是 (0,21) 和 (21,96)，分别表示整个占据轨道空间和整个空轨道空间（注意 Python 默认是从 0 开始计数的，且区间左闭右开）。 最后输出 .molden 文件，用于观看轨道。

可能有读者对 96 这个数字有疑问，对于 6-31G(d) 基组，高斯默认采用Cartesian 型 6D，而 PySCF 默认采用球谐型 5D，因此这里 PySCF 算的轨道和基函数均只有 96 个，而在高斯里有 102 个（每个C原子多出一个基函数，6 个 C 则多出 6 个基函数，96+6=102。 此处无基组线性依赖，因此分子轨道数目也是 102）。

若想在 PySCF 中也使用 Cartesian 型基函数，则需在基组后加上一行 mol.cart = True。 注意： 对于 PySCF-1.6.4 之前的版本，在使用 Cartesian 型时，输出的 molden 文件内容有误，需用 1.6.4（含）之后的版本。

3. Multiwfn

Multiwfn 做轨道局域化的详细内容可阅读其作者Sob的博文 http://sobereva.com/380

这里仅简要介绍一下。 做完 HF 计算之后得到 .fchk 文件，用 Multiwfn 打开（笔者写稿时用的是 2019-Jul-31 版），输入 19 选择轨道局域化，默认是 PM局域化且使用 Mulliken 电荷，可以通过输入 -6 改为其他局域化方法，如 Boys 局域化、或改为 PM 局域化且使用 Löwdin 电荷，界面上还有其他诸多选项，如是否局域化芯轨道等等，一看便知。

接着输入 1 表示只局域化占据轨道，输入 2 则表示分别局域化占据轨道和空轨道。 局域化完成之后轨道自动保存进 new.fch 文件（与之前载入的 .fchk 文件同一目录），且程序会自动载入局域化后的轨道，可输入 0 观看局域轨道。

4. GAMESS

GAMESS 做轨道局域化输入文件如下

 $CONTRL SCFTYP=RHF RUNTYP=ENERGY ICHARG=0 MULT=1 NOSYM=1LOCAL=BOYS $END $SYSTEM MWORDS=1000 $END $SCF DIRSCF=.TRUE. DIIS=.TRUE. $END $BASIS GBASIS=N31 NGAUSS=6 NDFUNC=1 $END $LOCAL FCORE=.F. $END $DATA C12H12C1      1[原子坐标]$END

所有的 $ 符号不能顶格，需空一格。 在 $CONTRL 里指定 LOCAL=BOYS, POP, RUEDNBRG 分别可做 Boys 局域化、PM 局域化和 ER 局域化。 局域化细节可以在 $LOCAL 中给出，例如是否冻芯 FCORE=.T. 或 .F.，只对部分轨道局域化 NINA 和 MOINA 关键词等。 GAMESS 默认只能对（部分或全部）占据轨道做局域化，笔者没找出如何用关键词做空轨道的局域化，若有小伙伴知道可以留言告知。

5. ORCA

ORCA 做轨道局域化输入文件示例（如 ben.inp）如下

%locLocMet FBT_core -1000Tol 1e-8OCC trueVIRT falseend! RHF 6-31G** xyz 0 1[原子坐标]*

在 %loc 中指定局域化细节，LocMet 表示选取局域化方法，常用的有 PM 和 FB两种，FB 也即 Boys 局域化。 除此之外还有 NEWBOYS, AHFB, IAOBOYS, IAOIBO 等等，不如前两种常用，但各有特点。 T_core 表示芯轨道能量截断阈值，低于此阈值的会被冻结、不参与局域化，而将其设为 -1000 则可认为是所有的芯轨道参与局域化。 Tol 是局域化函数值收敛限。 OCC/VIRT 表示占据轨道空间/空轨道空间，后面 true/false 表示是否局域化该空间。 其他详细信息请参看ORCA 手册。 完成后会生成两份轨道文件： ben.gbw（存放正则轨道）和ben.loc（存放局域轨道），可以利用 orca_2mkl ben -molden 命令将ben.gbw 转化为 molden 文件观看正则轨道。 而观看局域轨道则需将 ben.loc 后缀手动改为 gbw，例如 ben_boys.gbw，再用 orca_2mkl ben_boys -molden 命令生成相应的 molden 文件。 Multiwfn 软件就支持用 molden 文件观看轨道。

6. Molcas/OpenMolcas

Molcas/OpenMolcas 做轨道局域化输入文件示例（如 ben.input）如下

&GATEWAYCoord=ben.xyzBasis=6-31G*Group=C1
&SEWARD&SCF&LOCALISATIONBoys

其中原子坐标存放在另一个文件 ben.xyz（第一行写原子数，第二行写Angstrom，第三行开始写坐标）里。 局域化的所有参数和选项可以看 Molcas 手册 LOCALISATION 程序部分。 局域化方法的关键词有 Pipek-Mezey, Boys, Edmiston-Ruedenberg, Cholesky, PAO 五种可选，默认是 Pipek-Mezey。 可以用 Occupied, Virtual, All 指定局域化哪个子空间或全部轨道，默认是Occupied（PAO 默认是 Virtual）。 默认不冻结轨道，可以用 Nfrozen 指定冻结的轨道，如 Nfrozen=6 即表示冻结前 6 个芯轨道。

PS1： 本文均以 RHF 为例，DFT 的轨道也可以做局域化，过程一样；