Gaussian、ORCA使用笔记

发表于 更新于 分类于
关于使用Gaussian、ORCA时的一些学习经验和技巧

本文主要记录一些本人在学习、使用Gaussian、ORCA等量子化学软件以及各类辅助软件和脚本的过程中积累的一些经验，会不断更新。

用Gaussian寻找过渡态

过渡态，也即transition-state，可简写为TS。根据我个人的体验，过渡态的寻找算得上是很难的，所以也很值得一写。

关键词的选择

关于各种理论和方法的实现细节，我一点都算不上了解。有兴趣的可以参考Sobereva的《过渡态、反应路径的计算方法及相关问题》。总得来说，寻找过渡态时，应优先考虑这样的关键词，而不要轻易尝试QST2、QST3：

1
# opt=(TS,calcfc,noeigen) freq

同时，如果你的体系里溶剂化效应非常明显，也请记得加上溶剂化，因为溶剂化模型会改变势能面的形状，程度随体系而变。一般来说我们最常用的是SMD模型。但是有时由于SMD模型的背景噪音，会导致在过渡态搜索的时候一直不收敛，或者在同一振动路径上能收敛出多个结构。拿这样算出来的结构去算IRC一定是跑不出来的，应该尝试将溶剂化模型换为默认的IEFPCM重新算过渡态和IRC。

如果你在之后遇到了困难（IRC跑不出来等）需要重新算TS，由于之前已经得到了比较接近的结构，可以使用较贵的关键词来得到更精确的结构：

1
# opt=(TS,calcfc,tight) freq

如果还不能解决问题，则需要考虑一些更昂贵的关键词：

1
2
# opt=(TS,calcall,tight) freq
# opt=(TS,recalc=3,tight) freq

初猜的获取

关键词知道了，怎么样获得合理的初猜呢？

  1. 鉴于哈蒙德假说，反应的过渡态与反应物和产物中能量较高的构型相似。故而我们可以试图在能量较高的构型中，将你关心的那一根键的键长增加30%到50%。
  2. 如果你的反应的过渡态不是非常直观，不好通过你的化学直觉找出来的话，也可以尝试首先使用柔性扫描得到过渡态时键长的定性范围。但这个方法不普适，同时多变量扫描时计算量非常大，需谨慎考虑。

关于做柔性扫描的方法，首先采用# opt=modredundant关键词，随后在坐标末尾空一行，再输入：

1
[原子序号1]  [原子序号2]  {[原子序号3]}  {[原子序号4]} S  [步数]  [步长]

输入两个原子对应键长，三个对于键角，四个对应二面角。详细的内容可参照Sobereva的《详谈使用Gaussian做势能面扫描》。

如果你的计划中采用的是M06-2X这样的泛函，有时候会出现找不到想要的过渡态的情况，比如得到的结构的虚频对应的总是分子骨架的振动。这个时候你可以考虑将方法换成B3LYP，并先适当地降低溶剂化级别的等级，如用IEFPCM代替SMD或者干脆暂时不用溶剂化模型，得到合理的结果后再将得到的结构作为M06-2X的初猜。如果这样还找不到，我还用过一种不知道原理上可不可靠的方法：先将你关心的几个原子距离调到柔性扫描时得到的最佳距离，再冻结住这几个原子，用B3LYP进行限制性优化；再将得到的结构作为初猜，使用低级别的溶剂化模型和B3LYP找TS；最后将上一步得到的结构作为初猜，用你本来的溶剂化级别和M06-2X找“准确”的过渡态。

过渡态的验证：IRC

计算正常完成之后，我们再继续验证计算得到的结构是否是我们想要的结构。如果计算报错了，可以参考《Gaussian FAQ；新手求助报错时的注意事项》以及《量子化学计算中帮助几何优化收敛的常用方法》等文章中的方法尝试解决。

首先使用GaussView（或者免费的Avogadro等，都可以）打开输出文件，并进行以下检查：

  1. 分子构型是否是合理的，有没有优化出明显和这个反应无关的结构
  2. 是否只有一个虚频？这个虚频的振动方向是否是反应进行的方向？

满足以上两点之后，可以使用IRC来判断计算的结果到底是不是你想要的过渡态。

所需的关键词为：

1
# IRC=calcfc

按我个人的经验，跑IRC是非常容易报错的，大家可以参阅Sobereva的《在Gaussian中计算IRC的方法和常见问题》来解决问题。

用ORCA计算单点能

由于ORCA有强大的RI，使得我们可以在机器不变、耗时不增加的情况下，使用高得多的方法和基组来得到更加精确的单点能。关于ORCA的安装和基础使用，可以参阅《如何安装量子化学软件ORCA》和《量子化学软件ORCA快速入门指南》。

自然，这一步是在优化完成之后才进行的。但是如何快捷地将大量的文件从Gaussian输出文件转化为ORCA输入文件呢？

为了解决这个问题，我写了一个脚本，可以单独或批量地读取高斯输出的.log文件，转换为同名的.inp文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
@echo off
mode con cols=55 lines=10
Title ORCA Input Generator

:start
echo 请输入文件名(不含后缀名, 遍历请输入".")
set /p tsknum=文件名：
If "%tsknum%"=="" Goto Start
If /I "%tsknum%"!="." Goto A
If /I "%tsknum%"=="." Goto B

:A
cd /d D:\multiwfn
echo 100 >> oitemp.txt
echo 2 >> oitemp.txt
echo 12 >> oitemp.txt
echo E:\zeropoint\%tsknum%.inp >> oitemp.txt
echo 8 >> oitemp.txt
Multiwfn.exe "E:\opt\%tsknum%.log" < oitemp.txt
del /Q oitemp.txt
goto end

:B
for /R E:\opt %%f in (*.log) do ( 
cd /d D:\multiwfn
echo 100 >> oitemp.txt
echo 2 >> oitemp.txt
echo 12 >> oitemp.txt
echo E:\zeropoint\%%f.inp >> oitemp.txt
echo 8 >> oitemp.txt
Multiwfn.exe "E:\opt\%%f.log" < oitemp.txt
del /Q oitemp.txt
)
goto end

:end
exit

这是一个.bat文件，只能在Windows上使用，linux用户很方便可以找到功能更强的shell脚本请不要担心。请大家在使用前首先确保自己的电脑安装了最新的Multiwfn，同时在settings.ini中设定为isilent=1，iloadGaugeom=1；再根据自己电脑上各类软件的安装目录对脚本进行修改后再使用。

计算热力学量

怎么计算标况下气相下自由能，初学者总是没完没了犯错，故有必要再提一下。计算方式分为以下三种情况，结果准确度依次提升。注意这里B3LYP/6-31G*仅用来泛指中低计算级别，最适合用什么中低计算级别视具体体系而定，参看《谈谈量子化学中基组的选择》（http://sobereva.com/336）和《简谈量子化学计算中DFT泛函的选择》（http://sobereva.com/272）： (1)巨懒无敌懒人、本科生：# B3LYP/6-31G* opt freq，读Sum of electronic and thermal Free Energies=后面的值。由于此时电子能量部分是在中低便宜级别下算的，所以自由能精度很垃圾，数据无法用于发表。 (2)一般人：# B3LYP/6-31G* opt freq，读取吉布斯热校正量（Thermal correction to Gibbs Free Energy=后面的值）。然后在此结构下计算高精度单点能（比如用B2PLYP-D3(BJ)/def2-TZVP级别，关键词是B2PLYPD3/def2TZVP，或更好的比如CCSD(T)/cc-pVTZ）。之后将二者相加。 (3)讲究的人：先用B3LYP/6-31G*优化体系，然后在考虑B3LYP/6-31G*级别下的ZPE的频率校正因子、从0K升温对焓的改变的频率校正因子、熵的频率校正因子的情况下分别计算出吉布斯热校正量的相应部分，然后都加到高精度单点上。若不懂相关知识的话，仔仔仔细细细阅读《谈谈谐振频率校正因子》（http://sobereva.com/221）。

另外，一种粗略但比较省事的计算溶剂下焓或自由能的方法如下： (1)用opt freq优化结构同时获得焓或自由能的热校正量（获得焓的热校正量时建议用ZPE校正因子）。溶剂模型是否需要用看情况 (2)在(1)的结构下结合溶剂模型用高级别计算单点 然后将1和2的结果相加，再加上1.89 kcal/mol即可。此种做法步骤少，但相当于是在高级别下计算了溶解自由能，会比在溶剂模型专门参数化的级别下算的溶解自由能精度低。

由于我很懒，所以没用过热力学组合方法和频率校正因子。虽然只是将吉布斯自由能校正量加上单点能已经很简单了，不过在算的结构很多的时候，一个个提取数值并相加这种体力活也非常繁琐；所以我写了几个Python脚本（在3.7.4下运行正常，不知道对python2兼容性如何，请注意），目前已经放到GitHub上了，大家可以自行下载使用。

使用这几个脚本可以自动提取高斯输出文件中的吉布斯自由能校正量以及ORCA的单点能，同时自动进行相加，并将结果汇总到Excel中。该脚本还附带了核时统计功能。具体使用方法见GitHub页面上的README。

任之 wechat
订阅公众号以获取最新内容！