计算化学公社
标题: 将Gaussian与ORCA联用搜索过渡态、产生IRC、做振动分析 [打印本页]
作者Author: sobereva 时间: 2018-5-31 21:36
标题: 将Gaussian与ORCA联用搜索过渡态、产生IRC、做振动分析
将Gaussian与ORCA联用搜索过渡态、产生IRC、做振动分析
Use Gaussian with ORCA to search for transition states, generate IRC, and perform vibration analysis
文/Sobereva @北京科音 First release: 2018-May-31 Last update: 2022-Jul-13
1 前言
在《将Gaussian与Grimme的xtb程序联用搜索过渡态、产生IRC、做振动分析》（http://sobereva.com/421）中笔者介绍了Gaussian的external功能的使用，以及如何利用这个功能，借用xtb程序产生的能量、受力、Hessian来进行过渡态搜索、IRC、振动分析任务。类似地，只要自己写一个接口，也可以令Gaussian借用知名的、用户数目仅次于Gaussian的ORCA程序产生的这些信息来做这些任务，本文就介绍具体怎么实现。看本文前一定先把上文看了。ORCA的一些粗浅介绍可以参看《大体系弱相互作用计算的解决之道》（http://sobereva.com/214）。
这种Gaussian与ORCA联用的做法可以带来以下好处：
(1)可以令Gaussian做上述任务的时候支持ORCA才支持的理论方法，比如PBEh-3c、B97-3c、PWPB95-D3、NEVPT2、MRCI、DLPNO-CCSD(T)等。虽然其中有的理论方法并没有解析梯度，但原理上，我们可以通过自写脚本以有限差分方式得到Gaussian所需的梯度信息
(2)虽然Gaussian支持的DFT泛函很广，速度也快，但是由于ORCA的RI做得十分出色，借用ORCA来产生能量和导数可以使得上述任务耗时更低，尤其是对于大体系、大基组而言
(3)虽然ORCA也能直接做几何优化和找过渡态，但Gaussian在这方面算法上明显更成熟、更稳健，选项也更丰富、更易用
(4)ORCA目前没有IRC功能，和Gaussian联用使得ORCA中的理论方法也可以用来产生IRC
(5)使得ORCA产生的Hessian所对应的振动模式能够通过gview来观看
实际上，ORCA有个关键词ExtOpt，即在ORCA优化过程中通过外部文件读入能量、导数信息，这和Gaussian的external关键词颇为相似，但ORCA的optimizer的算法跟Gaussian比还是有差距，而且也只能用来优化，所以没太大意义。
本文使用Gaussian 09 E.01版，ORCA是4.0.1.2版，系统是CentOS 7.2。本文涉及的所有文件都可以在此处下载：http://sobereva.com/attach/422/gau_orca.zip
如果大家的实际研究中使用了本文的接口，请这样引用：Tian Lu, gau_orca: A Gaussian interface for ORCA program, http://sobereva.com/422 (accessed month day, year)
2 Gaussian与ORCA联用的external脚本的写法
文件包里的orca.sh就是笔者编写的接口脚本，Gaussian输入文件里写上external='./orca.sh'关键词就代表需要能量和导数信息的时候会调用当前目录下的orca.sh。下面解释一下此脚本的内容。
脚本一开始的部分，需要由用户自己填写ORCA运行时候是几核并行、内存用多大，用什么计算级别，以及数值精度方面的设定，比如SCF收敛限、积分格点精度。为了导数计算准确，脚本默认带了tightscf。还要设ORCA可执行文件路径。这些信息都需要用户使用前根据实际情况修改
#Set the number of parallel cores and maximum memory utilized by each core (MB)
nprocs=1
maxcore=1000
#Set calculation level
level="BLYP def2-svp def2/J"
#Set parameters for numerical aspects
numset="tightscf"
#ORCA executable path
orcapath="/sob/orca/orca"
之后，用read命令从Gaussian自动产生的InputFile中读取原子数、需要的导数阶数、净电荷、自旋多重度到相应变量
read atoms derivs charge spin < $2
下面的语句判断当前任务要用的关键词。ORCA里面engrad（即energy+gradient）关键词和Gaussian里的force关键词等价，用来计算当前结构下能量和受力，并输出到当前目录下与输入文件同名的.engrad文件中。如果写上freq，则会计算Hessian并且输出到当前目录下与输入文件同名的.hess文件中。这俩文件都是文本文件。
if [ $derivs == "2" ] ; then
task="engrad freq"
elif [ $derivs == "1" ] ; then
task="engrad"
fi
之后产生ORCA输入文件mol.inp。由于从InputFile读过来的坐标单位是Bohr，所以用了BOHRS关键词。
#Create ORCA input file
echo "Generating mol.inp"
cat >> mol.inp <<EOF
! $level $numset $task BOHRS nopop
%pal nprocs $nprocs end
%maxcore $maxcore
* xyz $charge $spin
$(sed -n 2,$(($atoms+1))p < $2 | cut -c 1-72)
*
EOF
值得一提的是，ORCA在计算上面语句产生的mol.inp的时候，会在当前目录下产生mol.gbw。ORCA默认是开启了autostart设定的，即如果计算时候发现当前目录下已经有了后缀为.gbw的与当前输入文件同名的文件，就会自动从中读取波函数作为初猜波函数。整个Gaussian任务在运行期间是把mol.gbw一直保留在当前目录的，因此如果Gaussian和ORCA联用做几何优化，ORCA在计算当前结构的时候会自动读取上一次运算时候产生的.gbw里的波函数当做初猜，这样一方面节约了计算时间（毕竟比重新产生的初猜好），另一方面有助于保持所处的电子态的连续性。（稍有Gaussian知识的人都知道，Gaussian做几何优化的时候，每一步的初猜用的就是上一步收敛的波函数，通过orca.sh将Gaussian与ORCA联用时也借助gbw文件同样发挥了这种效果）
再往后是调用ORCA运行mol.inp
#Run ORCA
echo "Running ORCA..."
$orcapath mol.inp > mol.out
echo "ORCA running finished!"
下面通过笔者用Fortran自写的extorca程序从ORCA产生的.engrad和.hess文件中提取能量以及导数信息，写入到与$3参数对应的OutputFile文件中
echo "Extracting data from ORCA outputs via extderi"
./extorca $3 $atoms $derivs
extorca的源程序是extorca.f90，内容并不复杂，就不多说了。
最后，orca.sh用以下命令清空ORCA运行中途产生的各种mol和mol_开头的文件以保持当前目录的清洁。但为了保留下来mol.gbw，用rm之前临时改了个名字。
mv mol.gbw tmp.gbw -f
rm mol.* mol_* -f
mv tmp.gbw mol.gbw -f
通过Gaussian与ORCA联用可以留下chk文件，但里面显然是没有波函数信息的，因此之后没法用Multiwfn等程序来看轨道、做波函数分析。不过，由于mol.gbw文件被保留了下来，这里面记录了ORCA计算最后一个结构时产生的波函数信息，因此可以转换为.molden格式，然后载入到Multiwfn里看轨道和做波函数分析。怎么把.gbw转换成.molden看《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。
3 实例：HCN->CNH异构化的过渡态搜索、振动分析以及IRC生成
在本文的压缩包的example目录里有HCN->CNH异构化的所有输入输出文件，用的是RI-BLYP/def2-SVP级别。由于这个体系非常小，用的计算级别也较低，因此实际上完全发挥不出借用ORCA计算能量和梯度的优势，反倒比纯粹用Gaussian还慢。而对于较大体系，ORCA仗着开挂般的RI，就完全不是同样的光景了。
首先运行找过渡态的任务，要把TS.gjf、orca.sh、extorca三个文件都放到当前目录下，并且根据实际情况恰当修改orca.sh开头的变量，特别是ORCA的可执行文件路径。TS.gjf中的坐标是初猜的过渡态结构，文件前三行为
%chk=TS.chk
%nproc=1
#P opt(nomicro,calcfc,ts,noeigen) external='./orca.sh'
用诸如g09 < TS.gjf |tee TS.out运行即可。注意Gaussian用external功能时应当确保以串行方式计算，因此刻意写了%nproc=1，否则在通过orca.sh调用ORCA期间，Link401会有很高无意义的CPU占用。
过渡态任务完成后，可以做振动分析检验有无虚频。对应的输入文件freq.gjf总共就三行，其它信息通过geom=allcheck从找过渡态后产生的TS.chk里读
%chk=TS.chk
%nproc=1
# freq geom=allcheck external='./orca.sh'
用gview观看输出文件freq.out，可确认有且只有一个虚频，振动模式正是期望的，虚频数值为1089.0cm-1。同样在BLYP/def2-SVP下，完全由Gaussian计算所给出的虚频为1088.4 cm-1，可见相符极好。
最后来跑一下IRC，对应的输入文件IRC.gjf内容也只有以下三行。此任务从TS.chk里读取过渡态结构，由于用了%oldchk又避免了此任务改写TS.chk。
%oldchk=TS.chk
%nproc=1
# IRC=calcfc geom=allcheck external='./orca.sh'
读者可用gview打开IRC.out绘制IRC曲线图，可以看到曲线很光滑，而且IRC轨迹正常，证明Gaussian和ORCA联用非常成功。
作者Author: beefly 时间: 2018-6-12 17:00
gaussian16发布了一个接口程序包，不知是否能实现类似的效果
http://gaussian.com/interfacing/
作者Author: 胡说 时间: 2018-12-13 09:43
请问sob老师， 如果优化想要考虑溶剂，是不是只在orca任务里加入溶剂模型就可以了？
作者Author: sobereva 时间: 2018-12-13 19:56
是
作者Author: 胡说 时间: 2019-1-7 14:42
谢谢老师。
另外还有个问题，ORCA在梯度计算时候也会有RI加速效果吗？
还有请问利用该方法是否也可以用来优化激发态呢？ORCA里加入tddft计算，指定优化的iroot，通过关键词engrad计算激发态梯度。
作者Author: sobereva 时间: 2019-1-7 23:30
有
TDDFT支持解析梯度
作者Author: 胡说 时间: 2019-1-8 21:23
谢谢老师
试了下，发现ORCA还不支持CPCM下的TDDFT解析梯度
作者Author: 胡说 时间: 2019-2-26 15:39
老师好！最近用这种方法来优化大体系结构。
比较高斯和orca优化任务发现，高斯scf-502算完单点，701/702/703模块算导数非常快，即使体系很大，是502的零头。虽然orca的scf借助RI可以比高斯快，但之后算梯度模块基本和scf时间差不多。比如一个约400原子体系，带溶剂，b3lyp/def2sv(p), 高斯scf花了90分钟，701/702/703加起来也就三四五分钟；ORCA开RIJCOSX，scf50几分钟，梯度约40几分钟。这样一来用ORCA杂化泛函来优化倒没什么优势（ORCA纯泛函blyp开RI还是很快的，大概20+17分钟）
不明白为什么高斯的能量一阶导数算得如此之快，时间和其scf相比都可以忽略了，几乎scf的时间就是一轮opt的时间。但orca一轮opt的时间大概是1.8~2倍的scf时间。理论上orca算导数支持RI不应该快些吗？
作者Author: sobereva 时间: 2019-2-26 16:31
(, 下载次数 Times of downloads: 158)
基组不够大时RIJCOSX的优势不显著
要图便宜用RI-纯泛函
作者Author: 胡说 时间: 2019-2-27 09:13
嗯 确实是。
老师 请问图中engrad的时间是包含了SP的时间的总计算时间吗？
作者Author: sobereva 时间: 2019-2-27 14:52
是
作者Author: gauss98 时间: 2019-6-13 15:24
老师好！
我按照文中安装后，还是没有运行成功
log文件显示：
External calculation of energy, first and second derivatives.
Running external command "./orca.sh R"
input file "/home/gauss/g16/tmp/Gau-2383.EIn"
output file "/home/gauss/g16/tmp/Gau-2383.EOu"
message file "/home/gauss/g16/tmp/Gau-2383.EMs"
fchk file "/home/gauss/g16/tmp/Gau-2383.EFC"
mat. el file "/home/gauss/g16/tmp/Gau-2383.EUF"
Failed to open output file from external program.
Error termination via Lnk1e in /home/gauss/g16/l402.exe at Thu Jun 13 23:16:40 2019.
Job cpu time: 0 days 0 hours 0 minutes 0.6 seconds.
Elapsed time: 0 days 0 hours 0 minutes 0.7 seconds.
File lengths (MBytes): RWF= 22 Int= 0 D2E= 0 Chk= 2 Scr= 2
好像是没有能打开output文件
在gauss的tmp文件夹里边 也只有 EIn 文件 没有 EOu文件
是哪里设得不对？
我的orca.sh 文件只改了下面部分：
#Set the number of parallel cores and maximum memory utilized by each core (MB)
nprocs=32
maxcore=3000
#Set calculation level
level="b97-3c"
#Set parameters for numerical aspects
numset="tightscf grid4"
#ORCA executable path
orcapath="/home/gauss/orca411/orca"
作者Author: sobereva 时间: 2019-6-13 22:38
确保orca独立运行正常，并且先测试我的脚本里原本的关键词
还搞不清楚就在.sh里加入一些调试语句帮助弄清楚
作者Author: gauss98 时间: 2019-6-17 15:37
老师好！折腾了好几天了，有一些进展，还是有些问题没有解决。
目前，测试文件都通过了。顺利的连用高斯和orca进行几何优化和过渡态计算。
但是，很奇怪的是，只能在集群主节点运行。 分析是 extorca文件只能在主节点运行，在计算节点算，都能顺利启动orca，但是到了提取数据这项 就报错退出。我采取了很多措施，包括 在。bashrc里边用
alias extorca='/home/gauss/bin/extorca'
后来又将程序拷到 orca目录里边(orca能在计算节点顺利运行
#orca
export PATH=$PATH:/home/gauss/orca411
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/gauss/orca411
alias orc='/home/gauss/orca411/orca'
alias extorca='/home/gauss/orca411/extorca'
都不能成功调用extorca文件。很是奇怪，不知道什么原因。
作者Author: moonlight2690 时间: 2019-11-11 19:07
您好，请问您上述问题如何解决的，我也遇到和你一模一样的问题，望指点一二，多谢
作者Author: gauss98 时间: 2019-11-15 11:12
不用了，继续用高斯
作者Author: moonlight2690 时间: 2019-11-15 17:12
我问题解决了，如果你后续还需要的话咱们可以讨论哈~
作者Author: 吴国卿 时间: 2019-11-26 10:31
你好，我也是卡在了L402报错了，一直解决不了，请问你是怎么解决的呢？谢谢
作者Author: skdmax 时间: 2019-12-18 17:07
你好,请问最后是怎么解决的?我还是报这样的错误.请您指导一下,谢谢
作者Author: 函数与激情 时间: 2020-2-16 14:16
本帖最后由 函数与激情 于 2020-2-16 14:17 编辑
sob老师您好，我将您文中的脚本稍作修改，试图使用openmolcas的CASPT2梯度用高斯的optimizer优化，修改后似乎是可以优化，但是并未在gaussina.log中更新优化圈数，直到最后收敛到四个YES后gaussian.log中也就只有这最后一圈的信息，log中只优化初猜结构和优化收敛的结构。不知道上述是在哪里出了问题，希望sob老师指点一下，不胜感激！
作者Author: 天南海北 时间: 2020-3-17 18:51
我也是这个问题，如何解决的？
作者Author: chenw 时间: 2020-9-8 21:31
请问，如果将opt和freq同时设定在orca.sh里，是不是不合理的？会不会每一步结构优化都算一遍freq?
作者Author: snljty 时间: 2020-9-8 22:37
1.opt和freq应该设在Gaussian的输入文件而不是orca.sh里，orca只负责提供能量和能量一阶二阶导数，优化算法等是Gaussian做的。
2. 应该不会，因为Gaussian每步external传递给外部的输入文件第一行都有个参数，指出这步任务要求做几阶导数。而orca.sh脚本里面也写了，只有什么情况才调用orca做一阶/二阶导数计算并回传信息。不需要的时候自然不会去算。
作者Author: chenw 时间: 2020-9-8 23:25
非常感谢，是我看错了；
此外我还有问题想请教一下，社长的orca.sh脚本 计算级别中没写RIJCOSX，我自己加上了，合理吗？
此外我还使用了溶剂模型，写的是CPCM(Toluene)
总的内容是：level="B3LYP D3 def2-svp def2/J RIJCOSX CPCM(Toluene)" 。怎么检验溶剂模型已被正确使用？
作者Author: sobereva 时间: 2020-9-9 06:34
想加就加
这些本来就是直接传递给ORCA的关键词，必然会被使用。不放心的话检查一下中间产生的ORCA输出文件便知
作者Author: 喵星大佬 时间: 2020-9-9 09:05
所以你后来问题解决了嘛
作者Author: 函数与激情 时间: 2020-9-9 14:51
后来我写信咨询了高斯官方，问题得到了解决。
这里的问题是因为External开发时其主要目的是ONIOM计算的low layer，其假设是所调用的计算非常迅速以至于其耗时几乎可以忽略。在这种情况下，如果像通常优化一样每步都要经过几个链接（Link），打印梯度、位移等收敛限信息，调用Link本身的时间虽然很低，但如果External调用的计算几乎不花时间，相对来说Link之间的传递会明显减慢计算速度。因此，使用External优化时默认采用类似分子力场计算的MicroIteration，全部在Link 402中完成，所以不会打印中间信息。External调用CASPT2的梯度计算，花费时间远远超出在几个链接之间跳转的时间，因此可以在优化中关闭 MicroIteration，关键词是“Opt=NoMicro”，这样每步优化都会经过几个Link并打印您要的信息。
作者Author: 喵星大佬 时间: 2020-9-9 16:37
我记得sob老师这个帖子里和联用xtb的帖子里都是说优化要加opt=nomicro的
也就是你那个脚本可以正常用高斯调用OpenMolca的高级别方法咯
作者Author: 函数与激情 时间: 2020-9-9 21:21
对的
作者Author: Scarlett-ww 时间: 2021-1-16 22:01
应该是你的orca目录写错了，再仔细检查一下~
作者Author: 任豹 时间: 2021-11-11 22:15
社长您好，单独的opt或者freq都可以进行计算，opt和freq一起时不能同时进行啊，这是因为软件受限制嘛？
作者Author: RAL 时间: 2021-11-11 22:21
因为Gaussian在做opt freq时，freq默认会从chk文件里读取波函数初猜，然后发现chk里什么都没有，报错...
作者Author: 任豹 时间: 2021-11-12 10:39
明白了，谢谢
作者Author: 任豹 时间: 2021-11-14 18:13
使用orca的版本5.0.1,进行过渡态搜索考虑溶剂模型时
想使用SMD模型!
在orca.sh添加关键词
%cpcm
smd true
SMDsolvent "water"
end
计算失败!
直接使用orca计算是成功的!
报错如下图1
关键词如下图2
但是使用CPCM模型
添加关键词
! CPCM(Water)
计算成功
这是怎么回事?
社长在http://sobereva.com/327中提到隐式溶剂模型下进行自由能计算中说绝对不要在Gaussian里用CPCM
我后面是要计算自由能的!
作者Author: sobereva 时间: 2021-11-15 02:50
之后你再单独用ORCA结合SMD算能量就完了
本来就不建议opt freq过程中用SMD。在Gaussian里用CPCM和在ORCA里用CPCM完全是两回事，而且对于opt freq的目的在ORCA里本来就建议用CPCM
作者Author: 任豹 时间: 2021-11-15 09:14
好的,谢谢sob社长!
作者Author: beefly 时间: 2021-11-19 12:19
从Gaussian 16b开始，opt和freq已经可以一起调用external算了
作者Author: beefly 时间: 2021-11-19 12:30
Gaussian 16.c有个重大修改（也可能是bug）。在oniom计算类型中，旧版EOu/EIn文件只包含external所对应层的原子。而G16.c要求必须包含全部原子（其他层的原子序数为0作为区分），所以目前的各种接口（orca、xtb、cfour、mm）在oniom计算中都有问题。
作者Author: mfdsrax2 时间: 2022-1-15 10:09
能在Windows下实现联用吗？
作者Author: sobereva 时间: 2022-1-15 20:38
不能（WSL、虚拟机另谈）
作者Author: StormSpirts 时间: 2022-8-26 19:17
虽然SMD不适合opt和freq，不过我改了一个能用SMD的脚本，大家有需要的自取吧。
作者Author: nirvana 时间: 2022-12-9 11:25
请问这种方法，版权上允许吗？（有Gaussian版权）
作者Author: 胡说 时间: 2022-12-13 13:21
老师，您好，请教一个问题。我发现通过Gau-ORCA接口来做freq和ORCA直接进行freq任务，两者输出的频率以及ZPE都是不一样的，低频的会差异大些，基于同样的Hessian，传给Gaussian后，为什么会得到和ORCA不一样的频率？
作者Author: sobereva 时间: 2022-12-14 17:33
谁管这个
作者Author: sobereva 时间: 2022-12-14 17:36
得看具体差多少
从Hessian转化到频率中牵扯到一些零碎细节，如坐标变换、单位转换、原子质量等，这些的差异都可能对结果产生一定影响
作者Author: wzkchem5 时间: 2022-12-14 19:09
当然允许，高斯的external optimization功能是高斯自己预留的，就是为了和别的程序接起来做结构优化等计算。如果这样做不允许的话，高斯一开始就不会做external optimization的功能。
此外，这个帖子是2018年的，随着ORCA这两年的发展，高斯ORCA联用的必要性越来越小了，例如ORCA已经实现了IRC功能，同时也支持了高斯不支持的NEB等方法可以高效地寻找过渡态。一般来说单用ORCA就足够了，高斯ORCA联用或许有时会比单用ORCA收敛快，但是因为ORCA的进步速度是大于高斯的，高斯ORCA联用比单用ORCA更快的情况应该是越来越少了。
作者Author: 胡说 时间: 2022-12-15 09:34
明白了，谢谢老师。差异是很小的。
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |