利用常见的程序做轨道局域化

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

1. Gaussian

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

%chk=ben_6-31Gd.chk
%mem=2GB
%nprocshared=4
#p RHF/6-31G(d)

Title Card Required

0 1
C     1.05889731   -0.35087719    0.00000000
C     2.45405731   -0.35087719    0.00000000
C     3.15159531    0.85687381    0.00000000
C     2.45394131    2.06538281   -0.00119900
C     1.05911631    2.06530481   -0.00167800
C     0.36151531    0.85709881   -0.00068200
H     0.50913831   -1.30319419    0.00045000
H     3.00356531   -1.30339019    0.00131500
H     4.25127531    0.85695381    0.00063400
H     3.00414131    3.01752581   -0.00125800
H     0.50899431    3.01758581   -0.00263100
H    -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, lo
from pyscf.tools import molden

mol = gto.M()
mol.atom = '''
[原子坐标]
'''

mol.basis = '6-31G(d)'
mol.max_memory = 2000 # MB
mol.verbose = 4
mol.build()

#RHF calculation
mf = scf.RHF(mol)
mf.kernel()

#localize the occ/unocc separately
occ_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_orb
vir_loc_orb = lo.Boys(mol, mf.mo_coeff[:,vir_idx]).kernel()
mf.mo_coeff[:,vir_idx] = vir_loc_orb
#localization done

#output Boys LMOs
molden.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=1
LOCAL=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
 C12H12
C1      1
[原子坐标]
$END

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

5. ORCA

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

%loc
LocMet FB
T_core -1000
Tol 1e-8
OCC true
VIRT false
end
! 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）如下

&GATEWAY
Coord=ben.xyz
Basis=6-31G*
Group=C1

&SEWARD
&SCF
&LOCALISATION
Boys

其中原子坐标存放在另一个文件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的轨道也可以做局域化，过程一样；而UHF和UDFT有alpha, beta两组轨道，理论上有分别在alpha占据轨道内做局域化、beta占据轨道内做局域化等等类别，较少用到。

PS2：一般而言，轨道局域化默认是分别在占据轨道空间与空轨道空间内部做的，不会涉及一个占据轨道与一个空轨道（即两个空间）之间的轨道旋转，因此不会改变HF能量。初学者可以读入局域轨道再算一次HF，可以发现SCF会1圈（有的程序是第0圈）收敛，收敛之后的轨道又会被重新对角化为离域的正则轨道。这便是HF能量的酉不变性。像在PySCF等程序中支持对指定的轨道作局域化，如果用户指定的轨道中包含了占据轨道和空轨道，则局域化后再算HF就不会一圈收敛了（这并不是缺点，有时会用到这种灵活功能）。类似地，对CASSCF分别在双占据空间/活性空间/空轨道空间每个空间内部做局域化，也不会改变CASSCF能量。

PS3：对于芯轨道（即最内层轨道如1s, 2s等）是否参与局域化，应视具体情况而定，如果要求所有的轨道都是局域的，那么芯轨道最好也参与局域化，因为芯轨道并不总是长得像原子轨道，比如该文中的苯，1-6号正则轨道其实是分布在6个碳原子上的，只有经过局域化才能各自归属给1个碳。如果所研究的问题不涉及芯轨道，则无需参与局域化。

PS4：一般只局域化占据轨道空间，很少直接对整个空轨道空间进行局域化，虽然硬要做也可以，不过对大体系、大基组这么做非常耗时，且局域出的空轨道局域程度较差。实际上有一些其他办法解决/绕开这个问题，以后有机会再介绍。

PS5：局域化方法的入门文献，强烈推荐Pipek和Mezey于1989年发表在JCP上的文章（DOI: 10.1063/1.456588），文中对几种局域化方法都做了简洁明了的介绍，比更早的一些Boys局域化文章讲得更为清楚。

相关程序的安装教程：