专栏首页量子化学从密度矩阵产生自然轨道_实战篇(上)

从密度矩阵产生自然轨道_实战篇(上)

本公众号之前发过自然轨道的原理介绍,详见《从密度矩阵产生自然轨道-理论篇》和《S^(1/2)的一些性质》。对其中原理和公式熟悉的读者,可以自己编写代码从密度矩阵产生自然轨道,无需阅读本文。笔者也写过这部分代码,可以提供一个参考

https://gitlab.com/jxzou/mokit/-/blob/master/src/lo.f90

见其中第一个subroutine no。

其实基本上每个量化程序自己都能产生自然轨道,不过一般在各软件手册里提及较少,导致有些人不知道、或者以为需要借助外部程序才能实现。本文就介绍一下如何在常见量化程序中获得自然轨道,以Gaussian软件为主,其他软件顺带提一下。RDFT/UDFT的操作与RHF/UHF的操作一样,不做单独说明。本文为上篇,分为4小节,分别是:

一、RHF自然轨道

二、UHF自然轨道

三、MP2和CCSD自然轨道

四、UMP2和UCCSD自然轨道

只对某一部分感兴趣的读者可以直接下划。更多内容请见下篇

1

RHF自然轨道

在之前的原理介绍篇里讲过,RHF占据轨道的任意酉变化均是其自然轨道,即自然轨道有无数组,这样的轨道没有什么实际意义。

2

UHF自然轨道

UHF自然轨道简称为UNO。以单重态O2为例,Gaussian输入文件示例如下

%chk=O2_uhf.chk
#p UHF/def2TZVP nosymm guess=mix stable=opt

sp of singlet O2 [标题行,随便写]

0 1
O   0.0   0.0   0.0
O   0.0   0.0   1.1616

--Link1--
%chk=O2_uhf.chk
#p chkbasis nosymm guess(read,only,save,NaturalOrbitals) geom=allcheck

可以看出,其基本逻辑是先做一个UHF计算,然后读取波函数(实际上是读取密度)产生自然轨道。chkbasis表示读取基组(若有赝势也可自动读),geom=allcheck表示读取电荷、自旋多重度和坐标,only表示不做能量计算,其他关键词顾名思义即可。获得自然轨道这一步的耗时一般只需几秒。另外,不习惯写--Link1--的读者可以拆成两个gjf文件来写。

注意,复杂体系的UHF/UDFT计算可能会得到内部不稳定(internal instability)的波函数,这些波函数做布局分析和产生自然轨道没有意义。因此这里先确保波函数是稳定的(即用stable=opt),然后才用来产生自然轨道。另外,对于过渡金属体系可能存在多个稳定的UHF解,能量各不一致,一般取能量最低的解来产生自然轨道。

对于单重态UHF计算,还应使用guess=mix来产生对称性破缺的初始猜测,对此不熟悉的小伙伴请阅读公众号发过的介绍《谈谈Gaussian软件中的guess=mix》。注意算完UHF后应先查看输出文件中自旋平方的期望值<S**2>,若其十分接近或等于零,此时UHF解实际上是RHF波函数,UHF能量等于RHF能量,即体系是闭壳层的,此时看自然轨道没有什么意义。

在计算完成后,用高斯自带的formchk小程序将chk文件转化为fchk文件,然后用GaussView打开fchk文件查看自然轨道(见下图)

(GV 6.0主面板上Tools->MOs->Visualize->选中或输入轨道序号->Update。注意不要点击Calculation。若Visualize是灰色的,说明电脑上没装高斯)

注意电子箭头旁的数据不是轨道能量,而是轨道占据数,为0~2之间的小数。若算出来有很明显的负数,说明自己很可能操作有误。自然轨道没有HOMO、LUMO这种说法,若非要一个称呼,可以称为HONO(最高占据自然轨道)、LUNO(最低未占据自然轨道)。从图中可以看出HONO、LUNO占据数都是1.0,是简并的,而其他自然轨道的占据数极其接近2/0,相当于RHF轨道。占据数情况表明,这是一个几乎完美的双自由基,双自由基特征y=1.0。关于双自由基我们以后再出一期详细的介绍。

当然,可视化轨道的软件远不止GaussView一个,读者可以用自己偏好的软件看轨道。高斯在UNO这个功能上有个缺憾(某种程度上可以说是个bug):原理上UNO轨道只有1列,不区分alpha与beta,但由于前一步是UHF任务(有两列轨道),高斯没能及时把chk文件里两列的信息删减为1列,导致产生的fchk文件里有两列一模一样的轨道,读者随便看1列即可。

其他软件就不会有这个问题,例如在OpenMolcas软件中,UHF计算完后默认就会产生UNO,存于后缀.UnaOrb的文件中,且只有1列轨道,占据数可以看#OCHR底下的数据。在GAMESS里是在$SCF中加上UHFNOS=.true.,产生的自然轨道及其占据数在.dat文件中的UHF轨道后。而在ORCA中则是在%scf里写上UHFNO true,结果在O2_uhf.uno文件中,这实际上也是一个gbw文件,你可以重命名为.gbw后缀后再通过ORCA自带的orca_2mkl小程序生成mkl和molden文件。

各个软件的关键词简要总结

有小伙伴可能听说过、或阅读本文后察觉到,UNO轨道可以作为CASSCF计算的初始轨道。由于有占据数大小排序(明显偏离2/0的可以认为是活性轨道),还可以根据轨道形状来进一步挑选,可以认为是一种不错的初始。更为详细的介绍留待以后的CASSCF计算技巧系列篇。

3

MP2和CCSD自然轨道

此处特指基于RHF的MP2和CCSD方法,基于UHF的请阅读下一部分。以基态平衡位置附近的N2分子为例,Gaussian输入文件示例如下

%chk=N2_cc-pVTZ.chk
#p MP2/cc-pVTZ density

title

0 1
N   0.0   0.0   0.0
N   0.0   0.0   1.08606

--Link1--
%chk=N2_cc-pVTZ.chk
#p chkbasis guess(read,only,save,NaturalOrbitals) geom=allcheck

唯一要注意的地方就是post-HF方法要加上density关键词,才会计算当前方法的密度,才能做布局分析和产生自然轨道,否则仍是HF下的结果。算CCSD的话就把MP2换成CCSD即可。高斯不支持CCSD(T)方法的密度,因此高斯无法产生CCSD(T)自然轨道。

对于OpenMolcas软件,RMP2计算默认不产生自然轨道,需要在&MBPT2中写上Prpt关键词(类似于高斯的density),生成的自然轨道存于后缀.MP2Orb的文件中。在GAMESS里则要在

各个软件的关键词简要总结

若无特别强调,表中均为弛豫(relaxed)密度。对于不优化轨道的方法(如MP2, CCSD, TDDFT等等),密度有弛豫与非弛豫密度(unrelaxed)之分,对应的自然轨道也略有不同。计算性质(如偶极矩、极化率等)推荐用弛豫密度,但其对应的自然轨道占据数理论上允许超出[0,2]的限制,可以说是个缺点。非弛豫密度则满足占据数[0,2]的限制。大部分程序只支持其中一种密度。这部分更详细的讨论请阅读

http://classic.chem.msu.su/cgi-bin/ceilidh.exe/gran/gamess/forum/?C34df668afbHW-7216-1405+00.htm

和ORCA手册(以4.2.1版为例)中8.1.3.2节Coupled-Cluster Densities。

另外,开源免费的PySCF程序也支持MP2和CCSD非弛豫密度和对应的自然轨道,不过它的Python语句其实是自己写公式算,而非简单的一两个关键词。这里我们举一例输入文件

from pyscf import gto, scf, mp
import numpy as np

mol = gto.Mole()
mol.build(
    atom = 'C 0 0 0; O 0 0 1.1281', # in Angstrom
    basis = 'ccpvtz',
    verbose = 4)
mf = scf.RHF(mol).run() # 算RHF能量

mypt2 = mp.RMP2(mf)
mypt2.frozen = 2
mypt2.kernel()  # 算MP2能量

rdm1 = mypt2.make_rdm1()  # 一阶约化密度矩阵
natocc, natorb = np.linalg.eigh(rdm1)
natorb = np.dot(mf.mo_coeff, natorb)
print(natocc)  # 打印MP2自然轨道占据数

4

UMP2和UCCSD自然轨道

仍以单重态O2为例,Gaussian输入文件示例如下

%chk=O2_uccsd.chk
#p UHF/def2TZVP nosymm guess=mix stable=opt

sp of singlet O2 [标题行,随便写]

0 1
O   0.0   0.0   0.0
O   0.0   0.0   1.1616

--Link1--
%chk=O2_uccsd.chk
#p UCCSD/chkbasis nosymm guess=read geom=allcheck density

--Link1--
%chk=O2_uccsd.chk
#p chkbasis nosymm guess(read,only,save,NaturalOrbitals) geom=allcheck

笔者特意从UHF开始写,意在提醒读者,应先检验UHF波函数稳定性,确保稳定了才能算UMP2或UCCSD。看过《谈谈Gaussian软件中的guess=mix》一文的读者肯定知道,对于单重态直接仅写个UMP2是没有意义的,算出来必然与RMP2一模一样,白费时间。

自然键轨道(NBO)不是自然轨道,不在本文介绍范围内。下期预告:

五、CASSCF自然轨道

六、CASCI自然轨道

七、自然跃迁轨道

八、激发态自然轨道

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

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

原始发表时间:2020-10-17

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 从密度矩阵产生自然轨道-理论篇

      对于一个单或多行列式波函数方法(例如RHF, MP2, CCSD, CASCI, CASSCF等等),可将电荷密度(charge density)

    用户7592569
  • 用ORCA做DLPNO-CCSD(T)计算

    DLPNO-CCSD(T)是一种基于局域轨道的近似CCSD(T)方法,其所得相关能与标准CCSD(T)的相关能比值可达99.9%以上,误差可控制在1 kcal/...

    用户7592569
  • Kuhn-Munkres配对算法

    生活或工作中,我们常常碰到分配问题。比如公司有n个任务,由n个工人来做,每个工人不同程度地擅长一个或几个任务。如果你是管理层,如何布置任务最大程度地发挥大家所长...

    用户7592569
  • Oracle Package的使用

    我们在Oracle的数据库里面在逻辑处理的时候可能会写大量的存储过程,由于数据多了以后,找起来比较麻烦,用package不仅能把存储过程分门别类,而且在pack...

    Vaccae
  • Excel实战技巧61: 处理剪切、复制和粘贴操作,使它们不会破坏已设置的单元格格式

    这是《Professional Excel Development》中介绍的一个技巧,特整理分享于此。

    fanjy
  • Python for 循环语句

    ? 文 | 糖豆 来源 | 菜鸟教程 糖豆贴心提醒,本文阅读时间4分钟,文末有秘密! Python for循环可以遍历任何序列的项目,如一个列表或者一...

    小小科
  • 11g Active DataGuard初探(r5笔记第54天)

    原本dataguard中日志应用和数据库只读查询是一个互斥的关系,两者不能并存。如果需要应用日志,则数据库只能在Mount状态下 使用recover manag...

    jeanron100
  • Linux 日常常用指令

       最近搞了一个阿里ECS,CentOS7,涉及到一些基本的Linux指令,在这里总结一下,在搭环境中常用的一些指令,熟悉这些指令就基本能够使用CentOS进...

    Rekent
  • SAP 销售条件表增强栏位

           有时遇到一个比较特殊的业务,比如公司间免费订单,既要让价格为0,不读取VK11里创建的价格, 又要让公司间的价格读取VK11,这实际上是有矛盾的,...

    SAP梦心
  • 一名Clojurian的Emacs配置

    我是一名热衷于函数式编程的Clojurian(Clojure粉),网络ID是lambeta(λβ),读作/‘læmeitə/,个人的博客网站是https://l...

    lambeta

扫码关注云+社区

领取腾讯云代金券