计算化学公社

 找回密码 Forget password
 注册 Register
Views: 16782|回复 Reply: 39

[ORCA] 将Gaussian与ORCA联用搜索过渡态、产生IRC、做振动分析

  [复制链接 Copy URL]

3万

帖子

99

威望

4万

eV
积分
81661

管理员

公社社长+计算化学玩家

发表于 Post on 2018-5-31 21:36:40 | 显示全部楼层 Show all |阅读模式 Reading model
将Gaussian与ORCA联用搜索过渡态、产生IRC、做振动分析

文/Sobereva @北京科音
First release: 2018-May-31  Last update: 2020-Feb-12

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 grid4。还要设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 grid4"
#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联用非常成功。

评分 Rate

参与人数
Participants 8
eV +36 收起 理由
Reason
北大-陶豫 + 5 好物!
zsu007 + 5
Allan + 3 牛!
shuimiaomiao110 + 3 牛!
alonewolfyang + 5
captain + 5 好物!
asdf + 5 谢谢
我本是个娃娃 + 5 好物!

查看全部评分 View all ratings

北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

568

帖子

18

威望

3866

eV
积分
4794

Level 6 (一方通行)

发表于 Post on 2018-6-12 17:00:53 | 显示全部楼层 Show all
gaussian16发布了一个接口程序包,不知是否能实现类似的效果
http://gaussian.com/interfacing/

40

帖子

0

威望

2797

eV
积分
2837

Level 5 (御坂)

发表于 Post on 2018-12-13 09:43:45 | 显示全部楼层 Show all
请问sob老师, 如果优化想要考虑溶剂,是不是只在orca任务里加入溶剂模型就可以了?

3万

帖子

99

威望

4万

eV
积分
81661

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2018-12-13 19:56:00 | 显示全部楼层 Show all
胡说 发表于 2018-12-13 09:43
请问sob老师, 如果优化想要考虑溶剂,是不是只在orca任务里加入溶剂模型就可以了?

北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

40

帖子

0

威望

2797

eV
积分
2837

Level 5 (御坂)

发表于 Post on 2019-1-7 14:42:47 | 显示全部楼层 Show all

谢谢老师。

另外还有个问题,ORCA在梯度计算时候也会有RI加速效果吗?
还有请问利用该方法是否也可以用来优化激发态呢?ORCA里加入tddft计算,指定优化的iroot,通过关键词engrad计算激发态梯度。

3万

帖子

99

威望

4万

eV
积分
81661

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2019-1-7 23:30:48 | 显示全部楼层 Show all
胡说 发表于 2019-1-7 14:42
谢谢老师。

另外还有个问题,ORCA在梯度计算时候也会有RI加速效果吗?


TDDFT支持解析梯度
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

40

帖子

0

威望

2797

eV
积分
2837

Level 5 (御坂)

发表于 Post on 2019-1-8 21:23:26 | 显示全部楼层 Show all
sobereva 发表于 2019-1-7 23:30

TDDFT支持解析梯度

谢谢老师
试了下,发现ORCA还不支持CPCM下的TDDFT解析梯度

40

帖子

0

威望

2797

eV
积分
2837

Level 5 (御坂)

发表于 Post on 2019-2-26 15:39:00 | 显示全部楼层 Show all
老师好!最近用这种方法来优化大体系结构。
比较高斯和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不应该快些吗?

3万

帖子

99

威望

4万

eV
积分
81661

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2019-2-26 16:31:53 | 显示全部楼层 Show all
胡说 发表于 2019-2-26 15:39
老师好!最近用这种方法来优化大体系结构。
比较高斯和orca优化任务发现,高斯scf-502算完单点,701/702/7 ...

Clipboard01.jpg

基组不够大时RIJCOSX的优势不显著
要图便宜用RI-纯泛函

北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

40

帖子

0

威望

2797

eV
积分
2837

Level 5 (御坂)

发表于 Post on 2019-2-27 09:13:00 | 显示全部楼层 Show all
sobereva 发表于 2019-2-26 16:31
基组不够大时RIJCOSX的优势不显著
要图便宜用RI-纯泛函

嗯 确实是。
老师 请问图中engrad的时间是包含了SP的时间的总计算时间吗?

3万

帖子

99

威望

4万

eV
积分
81661

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2019-2-27 14:52:20 | 显示全部楼层 Show all
胡说 发表于 2019-2-27 09:13
嗯 确实是。
老师 请问图中engrad的时间是包含了SP的时间的总计算时间吗?

北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

190

帖子

0

威望

1091

eV
积分
1281

Level 4 (黑子)

发表于 Post on 2019-6-13 15:24:02 | 显示全部楼层 Show all
老师好!
我按照文中安装后,还是没有运行成功
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"

                                                                                                   

3万

帖子

99

威望

4万

eV
积分
81661

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2019-6-13 22:38:30 | 显示全部楼层 Show all
gauss98 发表于 2019-6-13 15:24
老师好!
我按照文中安装后,还是没有运行成功
log文件显示:

确保orca独立运行正常,并且先测试我的脚本里原本的关键词
还搞不清楚就在.sh里加入一些调试语句帮助弄清楚
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班。这些培训是计算化学快速入门以及全面系统性提升研究水平的最佳途径,培训各种相关信息见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

190

帖子

0

威望

1091

eV
积分
1281

Level 4 (黑子)

发表于 Post on 2019-6-17 15:37:46 | 显示全部楼层 Show all
sobereva 发表于 2019-6-13 22:38
确保orca独立运行正常,并且先测试我的脚本里原本的关键词
还搞不清楚就在.sh里加入一些调试语句帮助弄 ...

老师好!折腾了好几天了,有一些进展,还是有些问题没有解决。

目前,测试文件都通过了。顺利的连用高斯和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文件。很是奇怪,不知道什么原因。



26

帖子

0

威望

542

eV
积分
568

Level 4 (黑子)

发表于 Post on 2019-11-11 19:07:02 | 显示全部楼层 Show all
gauss98 发表于 2019-6-13 15:24
老师好!
我按照文中安装后,还是没有运行成功
log文件显示:

您好,请问您上述问题如何解决的,我也遇到和你一模一样的问题,望指点一二,多谢

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )

GMT+8, 2022-5-22 18:13 , Processed in 0.276503 second(s), 27 queries .

快速回复 返回顶部 返回列表 Return to list