谈谈Gaussian软件中的guess=mix
笔者经常碰到小伙伴在用Gaussian软件计算涉及自由基的反应时，不清楚何时该加关键词guess=mix，何时不该加；也可能会有师兄/老师这样告诉新手：碰到自由基一律用guess(mix,always)。前者可能量化基础不扎实，碰到这类问题不懂；后者则可能缺乏实际计算经验。趁假期有空，正好写上一篇，详细解释一下。当然，笔者写的绝对不是标准答案，只能力求合理性和正确性，仅供对这个问题不清楚的小伙伴们参考。
为方便起见，本文仅讨论RHF和UHF，所有内容对DFT同样适用。对UHF计算不熟悉的新手可以阅读本公众号发过的软件教程《用Gaussian做UHF计算》。本文中所涉及计算皆使用G16 A.03。
首先给出简单的结论：（1）guess=mix只在自旋极化单重态（即使用UHF方法在单重态下做计算，发现有严重自旋污染）时需要考虑加，其他情况（如二重态、三重态）等无需考虑这个问题；（2）always表示在结构优化的每一步中都执行guess=mix。顾名思义，这只在结构优化中可能有用，而在单点计算中无需加、加了也没用。
接下来我们一一解释。关于自旋极化单重态和guess=mix的含义，Sob老师的博文《谈谈片段组合波函数与自旋极化单重态》（http://sobereva.com/82）讲得十分详细，强烈推荐经常做此类计算、但又还没看过博文的小伙伴仔细阅读。这里再举一个简单的单点计算实例：在UHF/STO-3G水平下计算键长为2.0 Å的H2分子，对于如下三种关键词写法：
(1) #p RHF/STO-3G nosymm(2) #p UHF/STO-3G nosymm(3) #p UHF/STO-3G nosymm guess=mix这样的UHF计算和结果显然没有太大价值，多花了计算时间却得到与RHF一样的解。而guess=mix的作用就是混合HOMO与LUMO轨道，使UHF的初始轨道与RHF不同，这样迭代就有可能收敛到更低的解上。
写法（3）得到的结果为−0.937213 a.u.，<S**2> = 0.946，能量明显低许多，代价是自旋不严格等于零，即不再是纯态。打开其log文件搜索Initial guess，会发现波函数初猜就已经不是纯态<S**2> = 1.0。
mix具体的做法并没有统一规定，可以是只混合alpha列的HOMO与LUMO轨道，也可以是alpha、beta两列的HOMO与LUMO轨道各自混合、但混合方式不一样。从Sob博文中展示的例子可以看出高斯是采用后者。当然，我们也可以自己手动操作一番，加深理解。例如，将下列计算产生的fchk文件打开
h2_uhf_only.gjf文件内容：
%chk=h2_uhf_only.chk#p UHF/STO-3G nosymm guess=(only,save)
Title Card
0 1H 0.0 0.0 0.0H 0.0 0.0 2.0
h2_uhf_only.fchk文件部分内容（修改前）：
Alpha MO coefficients R N= 4 6.68386924E-01 6.68386924E-01 7.53443037E-01 -7.53443037E-01Beta MO coefficients R N= 4 6.68386924E-01 6.68386924E-01 7.53443037E-01 -7.53443037E-01h2_uhf_only.fchk文件部分内容（修改后）
Alpha MO coefficients R N= 4 1.00538561E+00 -6.01437543E-02 -6.01437543E-02 1.00538561E+00Beta MO coefficients R N= 4 6.68386924E-01 6.68386924E-01 7.53443037E-01 -7.53443037E-01容易看出，上述改动仅是将alpha一列轨道的HOMO与LUMO轨道混合了一下。更具体地说，就是按照如下公式进行混合（此即等比例混合）：
再使用命令
unfchk h2_uhf_only.fchk h2_uhf_only.chk将修改后的fchk转回chk文件。将gjf文件内关键词guess=(only,save)改为guess=read，重新提交计算，即可得出较低的UHF能量−0.937213 a.u.。当然，等比例混合只是一种做法，也可以用 (根号3)/2、1/2这种搭配，保持轨道的正交归一性即可。高斯内部源代码实现的就是这种功能，只不过我们举了个最简单的例子，用计算器就能算。
但若混合比例很离谱，比如在这个例子里直接把HOMO、LUMO轨道互换（这相当于把波函数直接置于鞍点上）
h2_uhf_only.fchk文件部分内容（再一次修改）：
Alpha MO coefficients R N= 4 7.53443037E-01 -7.53443037E-01 6.68386924E-01 6.68386924E-01Beta MO coefficients R N= 4 6.68386924E-01 6.68386924E-01 7.53443037E-01 -7.53443037E-01再转化为chk文件，使用guess=read提交计算，也能立即收敛出一个能量−0.665399 a.u.，但明显高很多，说明这种混合不合理。
我们还可以举出很多guess=mix有效的例子，如单重态O2分子、单重态卡宾、并苯体系（并的环越多，自旋污染越大）等等。
高斯的guess=mix虽然很傻瓜、黑箱化，很多时候很好用，但它不能保证没有其他更低的UHF解，也不能保证收敛到一个稳定的UHF波函数，少数时候甚至连对称性破缺波函数都得不到。这些情况笔者在平时计算中也见过几例：
（1）O2在键长1.25 Å时，仅使用
#p UHF/cc-pVDZ nosymm guess=mix得到的电子能量为−149.589479 a.u.，对应<S**2>=0.627，确实是对称破缺的，有自旋污染。然而如果加上关键词stable=opt算，会发现输出文件内提示内部不稳定性
The wave function has an internal instability
并自动使用二阶优化方法继续优化至稳定，能量为−149.596725 a.u.，对应<S**2>=1.015，能量更低，自旋偏离纯态更多。实际上在这个例子里，O-O键长大于1.16 Å的几乎所有点都会有这个问题。
（2）自由基引发剂tBuOOH（叔丁基过氧化氢）分子，在M062X/6-31G(d)水平下优化得平衡几何结构，然后将O-O键长拉至1.90 Å，仅使用关键词
#p UM062X/6-31G(d) nosymm guess=mix算单点得到的电子能量为−308.590951 a.u.，对应<S**2> = 0，显然是收敛到了RM062X的解上。但是一旦加上stable=opt，会发现能量降为−308.594232 a.u.，<S**2> = 0.439。
所以笔者的建议是，用guess=mix时总是加上stable=opt，即检验波函数稳定性，若发现不稳定则优化至稳定。
如果仅写stable，表示只检验波函数稳定性，若稳定则结束；若不稳定只会输出不稳定信息然后结束，不会继续优化至稳定。这时候还要手动补交guess=read stable=opt任务，相当于波函数稳定性检验了两次，多费了些时间，因此不如直接用stable=opt简单。
有三点要注意：（1）高斯中stable=opt不能与结构优化，限制性优化或IRC等一起使用（除非自己写脚本实现，但可以预见计算量会很大）；（2）stable=opt只能保证得到稳定的波函数，对一些极特殊情况仍不能保证没有其他更低的UHF解，这一点以后再展开；（3）没有opt=stable这种东西。
诸如二重态、三重态这些情况，alpha、beta两列初始轨道可以一模一样，但由于两列轨道中的电子数不一样，对应的密度矩阵也不一样，因此不存在上述问题，无需guess=mix。
Reference
1. 《谈谈片段组合波函数与自旋极化单重态》http://sobereva.com/82
2. https://chemistry.stackexchange.com/questions/42005/initial-guess-for-unrestricted-hartree-fock-calculation