专栏首页量子化学谈谈Gaussian软件中的guess=mix

谈谈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

其中,写法(1)和(2)得到的电子能量是一样的,均为−0.783792 a.u.,<S**2> = 0,SCF迭代圈数均为1圈。打开(2)的log文件搜索Initial guess,会发现其波函数初猜是自旋纯态<S**2> = 0。这说明UHF默认的初始轨道与RHF的初始轨道是一样的,alpha、beta两列轨道一模一样。另外,由于是单重态,alpha、beta两列轨道里的电子数也是一样的。这些因素导致了迭代过程和结果也一样。

这样的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 1
H   0.0   0.0   0.0
H   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-01
Beta MO coefficients                       R   N=           4
  6.68386924E-01  6.68386924E-01  7.53443037E-01 -7.53443037E-01

h2_uhf_only.fchk文件部分内容(修改后)

Alpha MO coefficients                      R   N=           4
  1.00538561E+00 -6.01437543E-02 -6.01437543E-02  1.00538561E+00
Beta 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-01
Beta 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。

References

1. 《谈谈片段组合波函数与自旋极化单重态》http://sobereva.com/82

2. https://chemistry.stackexchange.com/questions/42005/initial-guess-for-unrestricted-hartree-fock-calculation

本文分享自微信公众号 - 量子化学(quantumchemistry),作者:jxzou

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-02-07

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • S^(1/2)的一些性质

      本文所述为量子化学电子结构理论中的基础知识,为本公众号同期另一文《从密度矩阵产生自然轨道_理论篇》一文的补充,对此基础内容熟悉的读者可以直接略过。

    用户7592569
  • 激发态计算中的溶剂效应

    已经有不少小伙伴在催更,非常感谢大家的支持,也有了更新的动力。关于隐式溶剂模型的介绍,可参见《理论计算中的溶剂效应模型》一文。本文着重介绍在激发态计算中使用隐式...

    用户7592569
  • 在Windows CMD里“使用”常见Linux命令

    相信不少小伙伴都曾经用过/偶尔使用Windows下的命令行终端(可按键盘组合键win+R然后输入cmd启动)

    用户7592569
  • 教你使用Python简单暴力爬取大量妹子图片

    当我们在我们的浏览器上输入www.baidu.com这个url后按下回车后,就向百度的服务器端发起请求,请求百度搜索的主页面资源,此时百度的服务器端收到请求,处...

    xujjj
  • SVD奇异值分解的数学涵义及其应用实例

    SVD(Singular Value Decomposition, 奇异值分解)是线性代数中既优雅又强大的工具, 它揭示了矩阵最本质的变换. 使用SVD对矩阵进...

    SIGAI学习与实践平台
  • python:解析URL

    这个函数的性能实在太差了。10000次用了整整45s。 在不严格的情况下,自己用split进行判定会好很多。快了12倍。

    超级大猪
  • underscore.js,jquery.js源码阅读 杨龙飞

    windseek
  • LARAVEL的那个RCE最有趣的点在这里

    直接把上述请求body中的viewFile参数的值替换为一个恶意ftp地址就可以实现rce

    tnt阿信
  • Django搭建博客(九):为博客添加代码高亮显示和 md文档支持

    特别需要注意的是:代码块必须使用三个 '`' 符号包裹起来才能正确识别,语言标记可有可无,但是三个 '`' 必须单独成行。

    渔父歌
  • mc参数备忘&java-json备忘

    mc参数(摘自 http://www.blogjava.net/jzone/articles/302991.html) 查看方法 telnet进去 或 ech...

    财主刀刀

扫码关注云+社区

领取腾讯云代金券