摘抄：Gaussian、ORCA使用笔记

发表于 更新于 分类于 阅读次数：
还没用过

关键词的选择

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

# opt=(TS,calcfc,noeigen) freq

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

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

# opt=(TS,calcfc,tight) freq

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

# opt=(TS,calcall,tight) freq

# opt=(TS,recalc=3,tight) freq

初猜的获取

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

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

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

[原子序号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 来判断计算的结果到底是不是你想要的过渡态。

所需的关键词为：# IRC=calcfc

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

【后面两个part因为无关所以不做调整】

用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；再根据自己电脑上各类软件的安装目录对脚本进行修改后再使用。

计算热力学量

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

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

原文地址：https://alexander-qi.github.io/2020/gaussnote/