片段组合波函数实例1.双原子分子

使用片段组合波函数作为SCF初猜是一种略为高级的计算技巧。笔者刚入门量子化学一两年时知道这个技巧时是有点震惊的：常规HF/DFT计算并不总是一个黑箱（black box），不是交个任务看到normal termination就完事了。特殊情况下需要计算者拥有一定的量子化学基本知识，并将其用于指定合适的SCF初猜，仔细分析计算结果。后来键解离和过渡金属算多了之后，笔者对这些情况就习以为常了。

片段组合波函数起初没有准确的一个英文单词来表达，Gaussian官网[1]将其描述为“Generate a guess built from fragment guesses or SCF solutions”，官网其他处[2]有时也简称为fragment guess，可以简单翻译为片段初猜。笔者第一次见到较为贴切的中文表达是在Sobereva的博文《谈谈片段组合波函数与自旋极化单重态》[3]中，他将其称之为“片段组合波函数”，较为形象。如果有读者见过这个词更早的中文表达，欢迎留言告诉大家。下文中我们将不加区分地使用片段初猜和片段组合波函数两个词。

一个十分经典的例子就是氮气N2分子单重态下N≡N三键解离的计算，这里我们先展示一个单点计算，文末再讨论解离曲线扫描。本文算例均使用UHF方法，所呈现的问题在UDFT里也存在，而所展示的计算技巧也同样适用于UDFT。以原子间距离2.0 Å的N2分子为例，一个表面上看似合理的Gaussian输入文件如下

%chk=N2_cc-pVDZ_2.0.chk
%mem=4GB
%nprocshared=4
#p UHF/cc-pVDZ nosymm guess=mix stable=opt

title

0 1
N   0.0   0.0   0.0
N   0.0   0.0   2.0

正常结束后可以在输出文件中找到两处SCF迭代收敛信息：

 SCF Done:  E(UHF) =  -108.616833307     A.U. after   16 cycles
 SCF Done:  E(UHF) =  -108.659688488     a.u. after    6 cycles

第一处SCF Done表示在guess=mix产生的SCF初猜下成功收敛，然而接着自动检验波函数稳定性表明该波函数其实是不稳定的。由于我们写的是stable=opt，碰到不稳定波函数会继续优化至稳定，所以就有了第二个SCF Done，明显能量更低。那么这就是单重态下键长2.0 Å的N2分子的最低能量了么？我们看一下自旋布居再说。

直接在输出文件中查找Mulliken charges and spin densities:即可。注意两个SCF Done对应两处自旋布居的信息，显然应该找第二处。可以发现

Mulliken charges and spin densities:
               1          2
     1  N    0.000000  -0.861075
     2  N   -0.000000   0.861075

第一列是Mulliken电荷，不是我们关心的。第二列是各个原子上的自旋布居，正负号分别代表alpha/beta自旋。注意不少人（尤其是实验为主、计算为辅的课题组）将此处照搬直译为自旋密度，这是错误的！

可以看到每个N原子只有接近1个电子的净自旋，这不太符合基本知识和化学直觉。因为N2平衡结构是N≡N三键，而解离时左右各自都应该像孤立N原子，所以无论是看平衡结构还是看解离极限，每个原子应各提供3个单电子。N原子不像过渡金属，旁边没有配体，不可能转移部分单电子到配体上。另外，根据洪特规则N原子中3个单电子应以自旋平行的方式分别占据不同的轨道，原子的能量最低。综上每个N原子应该是接近3个电子的净自旋更合理，如下图所示

既然我们有此猜想，能否将其直接写进输入文件中？Gaussian提供了这种功能，输入文件示例如下

%chk=N2_cc-pVDZ_2.0_frag.chk
%mem=4GB
%nprocshared=4
#p UHF/cc-pVDZ nosymm guess(fragment=2)

title

0 1 0 4 0 -4
N(fragment=1)   0.0   0.0   0.0
N(fragment=2)   0.0   0.0   2.0

--Link1--
%chk=N2_cc-pVDZ_2.0_frag.chk
%mem=4GB
%nprocshared=4
#p UHF chkbasis nosymm guess=read geom=allcheck stable=opt

在一开始的任务中，我们将体系划分为2个片段，通过guess(fragment=2)体现。在书写电荷和自旋多重度处，要依次写明整体的电荷，整体的自旋多重度，片段1的电荷，片段1的自旋多重度，片段2的电荷。。。依次类推。这个问题里整体、片段电荷都是零，没什么好说的。我们认为两个原子各有3个单电子，同时总自旋为单重态，所以一个原子的自旋多重度是4，另一个是-4，正负号用来表示alpha/beta自旋，保证整体是单重态，把两个原子的自旋多重度互换一下不影响计算结果。

第一个任务仅是计算每个片段的SCF并将其组装为整体的波函数，保存进chk文件作为整体SCF的初猜；我们还需通过Link1写上第二步任务，即读取构造好的初猜进行SCF计算，并检验波函数稳定性。如果读者细心翻找，会发现输出文件中有3处SCF收敛信息

 SCF Done:  E(UHF) =  -54.3911145622     A.U. after    7 cycles
 SCF Done:  E(UHF) =  -54.3911145622     A.U. after    7 cycles
 SCF Done:  E(UHF) =  -108.769405741     A.U. after   11 cycles

前两处即是孤立N原子的四重态UHF计算，两个原子能量相等，符合预期。最后一处是整体的UHF计算，该能量明显比前述-108.659688 a.u.要低非常多。那这个能量是合理的么？我们仍旧先看自旋布居，去输出文件里找到最后一处

Mulliken charges and spin densities:
               1          2
     1  N   -0.000000   2.901523
     2  N   -0.000000  -2.901523

每个原子的自旋布居均接近-3/3，说明一个原子有3个alpha单电子，另一个有3个beta单电子，符合我们的预期。这个例子目前的计算信息告诉我们：

（1）对于特殊体系（例如键解离），一定要检验波函数稳定性；

（2）一个体系可能拥有不止一个稳定的波函数，可以通过构建合适的SCF初猜（例如片段组合波函数）找到其他稳定解。

有的读者可能会问：那我不研究化学键解离，就研究一些平衡结构附近的分子，还要考虑这种问题么？答案是需要的。著名强关联体系Cr2分子，其中存在Cr-Cr六重键，在任何键长下都要用片段组合波函数构建初始猜测（否则容易算出偏高的能量）：

%chk=Cr2_cc-pVTZ_frag.chk
%mem=48GB
%nprocshared=24
#p UHF/cc-pVTZ nosymm guess(fragment=2) scf(maxcycle=128,xqc)

title

0 1 0 7 0 -7
Cr(fragment=1)  0.0   0.0   0.0
Cr(fragment=2)  0.0   0.0   2.0

--Link1--
%chk=Cr2_cc-pVTZ_frag.chk
%mem=48GB
%nprocshared=24
#p UHF chkbasis nosymm guess=read geom=allcheck scf(maxcycle=128,xqc) stable=opt

那读者可能又会问：N2和Cr2的例子有规律可循，分成的片段都是2个原子，各个键长下都采用类似的片段组合方式（N有3个单电子，Cr有6个单电子），这是普遍规律么？答案是：对于双原子分子确实经常是适用的。

不幸的是，存在极其特殊的情况：两个碳原子构成的C2分子，在不同键长下可能需要不同的片段组合方式。这个分子中存在的是C=C双键还是四重键，曾有广泛的争议和讨论，这里笔者无意去涉及这个争议，只想向大家介绍与本文有关的介绍技巧。对于单重态C2分子，显然有两种划分片段的方式：

(1)两个C原子均为三重态

%chk=C2_cc-pVTZ_2.0_frag1.chk
%mem=4GB
%nprocshared=4
#p UHF/cc-pVTZ nosymm guess(fragment=2)

title

0 1 0 3 0 -3
C(fragment=1)   0.0   0.0   0.0
C(fragment=2)   0.0   0.0   2.0

--Link1--
%chk=C2_cc-pVTZ_2.0_frag1.chk
%mem=4GB
%nprocshared=4
#p UHF chkbasis nosymm guess=read geom=allcheck stable=opt

(2)两个C原子均为五重态

%chk=C2_cc-pVTZ_2.0_frag2.chk
%mem=4GB
%nprocshared=4
#p UHF/cc-pVTZ nosymm guess(fragment=2)

title

0 1 0 5 0 -5
C(fragment=1)   0.0   0.0   0.0
C(fragment=2)   0.0   0.0   2.0

--Link1--
%chk=C2_cc-pVTZ_2.0_frag2.chk
%mem=4GB
%nprocshared=4
#p UHF chkbasis nosymm guess=read geom=allcheck stable=opt

在2.0 Å键长下，第（1）种划分方式得到的UHF能量更低，这符合我们的预期，解离之后更像两个C原子。然而在键长较短时（例如1.2 Å），以第（2）种划分方式得到的UHF能量更低。注意某一键长下的C2分子基态未必是单重态，“讨论C2基态自旋多重度是什么”偏离了本文核心内容，笔者用单重态举例仅是示例。

以上举的例子都是单点计算。对于键长扫描，笔者展示四种方案：

（1）最稳妥但操作略繁琐：不用量化软件自带的扫描功能，每个结构单独算单点，使用片段组合波函数做初猜，检验波函数稳定性，分析结果。

（2）稳妥且偷懒：使用量化软件自带的扫描功能从距离长扫至短，在UHF/UDFT下对第一个结构采用片段组合波函数作为SCF初猜、检验波函数稳定性后，读取该波函数进行后续的键长扫描。对常见的单键、多重键解离一般都能得到正确的解离曲线；而对极个别体系（例如C2分子），不同键长下可能需要不同的片段组合方式作为初猜，此时只能采用方式上述（1）。

（3）不靠谱：用UHF/UDFT配合guess(mix,always)关键词做解离曲线扫描，一般算单键解离的正确性还行，算多重键解离大概率会得到能量跳跃的曲线，不是最佳方案。

（4）错误：用RHF/RDFT直接从平衡位置扫描至解离，解离端的曲线是定性错误的。

当然，本文还有不少问题没有回答：

（1）双原子分子的电子态一般在实验中和文献上都有明确的对称性指认和归属。而本文在UHF/UDFT下结合片段组合波函数初猜，还加了nosymm，这能看出对称性么？

（2）本文介绍的体系涉及强关联问题，要定量算准解离曲线、基态与低激发态能级差这些问题，依靠UHF/UDFT及TDDFT一般是不行的。那该怎么算呢？

（3）含过渡金属的分子怎么划分片段，怎么指认片段自旋？反铁磁耦合的计算如何实现？

这些内容留待以后再介绍。

Reference

1. http://gaussian.com/guess

2. http://gaussian.com/afc

3. http://sobereva.com/82