计算化学公社

 找回密码 Forget password
 注册 Register
Views: 8229|回复 Reply: 2

[ORCA] 在ORCA中做counterpoise校正并计算分子间结合能的例子

[复制链接 Copy URL]

4万

帖子

99

威望

4万

eV
积分
89946

管理员

公社社长+计算化学玩家

发表于 Post on 2020-3-13 08:25:07 | 显示全部楼层 Show all |阅读模式 Reading model
在ORCA中做counterpoise校正并计算分子间结合能的例子

文/Sobereva@北京科音
First release: 2020-Mar-13  Last update: 2022-Jan-9


在《谈谈BSSE校正与Gaussian对它的处理》(http://sobereva.com/46)一文中介绍过BSSE问题与counterpoise校正。Counterpoise校正在高精度弱相互作用计算时通常需要考虑。Gaussian有现成的counterpoise关键词,因此做counterpoise校正很简单。但是ORCA程序(至少对于撰文时最新的4.2.1版来说)没有提供现成的counterpoise关键词,只能通过定义鬼原子的方式手工实现,因此实现起来稍微麻烦。不熟悉ORCA者看http://sobereva.com右侧ORCA分类里我写过的大量相关文章。

之前在《详谈Multiwfn产生ORCA量子化学程序的输入文件的功能》(http://sobereva.com/490)中笔者专门介绍了Multiwfn产生ORCA量子化学程序的输入文件的功能。为了使得用ORCA做counterpoise任务更方便,笔者于2020-Mar-13给Multiwfn的这个功能增加了相应选项,在这里通过一个实例简单演示一下。最新版Multiwfn可以在http://sobereva.com/multiwfn免费下载,相关知识见《Multiwfn FAQ》(http://sobereva.com/452)。

注意按照下文例子算出来的实际上是复合物结构中单体间相互作用能。而一般意义的结合能在计算时是要考虑单体形成复合物过程中形变造成的能量改变(变形能)。由于下面的例子里单体是刚性的,变形能可以忽略不计,因此下文仍称作结合能。

本文用的示例体系是下图这个体系,我们要在DLPNO-CCSD(T)/aug-cc-pVTZ下计算其结合能,并考虑counterpoise校正。


1.png

此体系的xyz格式的结构文件,后文中涉及的输入输出文件以及要用的shell脚本都可以在这里下载:http://sobereva.com/attach/542/file.rar。并不是必须得用xyz格式,只要是Multiwfn支持的含有结构信息的格式(如pdb、mol2...)都可以作为输入文件,见《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》(http://sobereva.com/379)。

启动Multiwfn,输入dimer.xyz的路径载入之。如果你还不知道两个单体的原子序号范围,可以进入Multiwfn的主功能0来查询。进入主功能0后,选择Tools - Select fragment,输入其中一个片段中的任意一个原子的序号比如16,整个片段就被高亮,而且序号也显示出来了,如下所示。

2.png

如此这般,把两个片段序号范围都拷出来备用,然后点Return退出主功能0,之后依次输入以下内容:
oi  //产生ORCA的输入文件
[按回车] //输出文件与输入文件同名但用inp作为后缀
0  //选择任务类型
7  //考虑counterpoise校正的结合能
1-15  //片段1的序号
16-32  //片段2的序号
-2  //加弥散函数
8  //用DLPNO-CCSD(T) with NormalPNO

在当前目录下马上就有了dimer.inp,内容如下

%pal nprocs   4 end
! DLPNO-CCSD(T) normalPNO RIJK aug-cc-pVTZ aug-cc-pVTZ/JK aug-cc-pVTZ/C tightSCF noautostart miniprint nopop
%maxcore 1000
* xyz   0   1
C      0.79991408   -1.02205164    0.68773696
H      0.85355588   -1.12205101   -0.39801435
H      1.49140210   -1.74416936    1.11972040
C      1.11688700    0.42495279    1.09966205
H      1.83814230    0.89014504    0.43045256
H      1.55556959    0.43982464    2.09708356
C     -0.24455916    1.16568959    1.10297714
H     -0.25807760    2.00086313    0.40532333
H     -0.44880450    1.57699582    2.09098447
C     -1.29871418    0.10381191    0.73930899
H     -1.47356078    0.10524338   -0.33800545
H     -2.25673428    0.27804118    1.22715843
C     -0.64687993   -1.22006836    1.13630660
H     -1.12443918   -2.08762702    0.68299327
H     -0.68601864   -1.34528332    2.22022006
C      0.04984615    0.09420760    5.61627735
C     -0.04649805   -0.05787837    7.13191782
H      0.94604832   -0.07334458    7.58427505
H     -0.60542282    0.77000613    7.57035274
H     -0.55366275   -0.98654445    7.39726741
C      0.76389939    1.40111272    5.28065247
H      0.84541894    1.53461185    4.20097059
H      0.22042700    2.25580115    5.68615385
H      1.77150393    1.41176313    5.69888547
C     -1.35516567    0.11403225    5.01895782
H     -1.31823408    0.23122219    3.93510886
H     -1.93746520    0.94145581    5.42730374
H     -1.88506873   -0.81375459    5.24028712
C      0.83774596   -1.07927730    5.03893917
H      0.34252564   -2.02626804    5.25918232
H      0.93258913   -0.99209454    3.95580439
H      1.84246405   -1.11668194    5.46268763
*

$new_job
! DLPNO-CCSD(T) normalPNO RIJK aug-cc-pVTZ aug-cc-pVTZ/JK aug-cc-pVTZ/C tightSCF noautostart miniprint nopop Pmodel
%maxcore 1000
* xyz   0   1
C      0.79991408   -1.02205164    0.68773696
H      0.85355588   -1.12205101   -0.39801435
H      1.49140210   -1.74416936    1.11972040
C      1.11688700    0.42495279    1.09966205
H      1.83814230    0.89014504    0.43045256
H      1.55556959    0.43982464    2.09708356
C     -0.24455916    1.16568959    1.10297714
H     -0.25807760    2.00086313    0.40532333
H     -0.44880450    1.57699582    2.09098447
C     -1.29871418    0.10381191    0.73930899
H     -1.47356078    0.10524338   -0.33800545
H     -2.25673428    0.27804118    1.22715843
C     -0.64687993   -1.22006836    1.13630660
H     -1.12443918   -2.08762702    0.68299327
H     -0.68601864   -1.34528332    2.22022006
C:     0.04984615    0.09420760    5.61627735
C:    -0.04649805   -0.05787837    7.13191782
H:     0.94604832   -0.07334458    7.58427505
H:    -0.60542282    0.77000613    7.57035274
H:    -0.55366275   -0.98654445    7.39726741
C:     0.76389939    1.40111272    5.28065247
H:     0.84541894    1.53461185    4.20097059
H:     0.22042700    2.25580115    5.68615385
H:     1.77150393    1.41176313    5.69888547
C:    -1.35516567    0.11403225    5.01895782
H:    -1.31823408    0.23122219    3.93510886
H:    -1.93746520    0.94145581    5.42730374
H:    -1.88506873   -0.81375459    5.24028712
C:     0.83774596   -1.07927730    5.03893917
H:     0.34252564   -2.02626804    5.25918232
H:     0.93258913   -0.99209454    3.95580439
H:     1.84246405   -1.11668194    5.46268763
*

$new_job
! DLPNO-CCSD(T) normalPNO RIJK aug-cc-pVTZ aug-cc-pVTZ/JK aug-cc-pVTZ/C tightSCF noautostart miniprint nopop Pmodel
%maxcore 1000
* xyz   0   1
C:     0.79991408   -1.02205164    0.68773696
H:     0.85355588   -1.12205101   -0.39801435
H:     1.49140210   -1.74416936    1.11972040
C:     1.11688700    0.42495279    1.09966205
H:     1.83814230    0.89014504    0.43045256
H:     1.55556959    0.43982464    2.09708356
C:    -0.24455916    1.16568959    1.10297714
H:    -0.25807760    2.00086313    0.40532333
H:    -0.44880450    1.57699582    2.09098447
C:    -1.29871418    0.10381191    0.73930899
H:    -1.47356078    0.10524338   -0.33800545
H:    -2.25673428    0.27804118    1.22715843
C:    -0.64687993   -1.22006836    1.13630660
H:    -1.12443918   -2.08762702    0.68299327
H:    -0.68601864   -1.34528332    2.22022006
C      0.04984615    0.09420760    5.61627735
C     -0.04649805   -0.05787837    7.13191782
H      0.94604832   -0.07334458    7.58427505
H     -0.60542282    0.77000613    7.57035274
H     -0.55366275   -0.98654445    7.39726741
C      0.76389939    1.40111272    5.28065247
H      0.84541894    1.53461185    4.20097059
H      0.22042700    2.25580115    5.68615385
H      1.77150393    1.41176313    5.69888547
C     -1.35516567    0.11403225    5.01895782
H     -1.31823408    0.23122219    3.93510886
H     -1.93746520    0.94145581    5.42730374
H     -1.88506873   -0.81375459    5.24028712
C      0.83774596   -1.07927730    5.03893917
H      0.34252564   -2.02626804    5.25918232
H      0.93258913   -0.99209454    3.95580439
H      1.84246405   -1.11668194    5.46268763
*


这个输入文件包含了DLPNO-CCSD(T)/aug-cc-pVTZ with normalPNO级别下的三个单点任务。第一个子任务是对二聚体做计算,输出的能量我们叫E_AB;第二个子任务是对第1个片段在二聚体基组下做计算,输出的能量我们叫E_A(AB);第三个子任务是对第2个片段在二聚体基组下做计算,输出的能量我们叫E_B(AB)。原子名后面带冒号的说明这个原子是鬼原子。计算前别忘了把maxcore和nproc都设为恰当情况。其中nproc只要在开头定义一次就够了。

在笔者的2*Intel E5-2696v3(36核)机子上,通过36核并行,每个核给6000MB内存,花了41分钟跑完。在输出文件中搜索FINAL SINGLE POINT ENERGY,总共会搜索到三个:
FINAL SINGLE POINT ENERGY      -393.613585294163
FINAL SINGLE POINT ENERGY      -196.197449643018

FINAL SINGLE POINT ENERGY      -197.412738820804
依次为E_AB、E_A(AB)、E_B(AB)。然后按照计算E_AB - E_A(AB) - E_B(AB)就得到考虑counterpoise校正的结合能了。如果不理解这点,看我的培训班(http://www.keinsci.com/workshop/KBQC_content.html)里的一页ppt
3.png

对于当前例子,counterpoise校正的结合能因此为627.51*(-393.613585294163+196.197449643018+197.412738820804)= -2.13 kcal/mol。

本文的这个体系其实是著名的S66弱相互作用测试集里的,用的结构也是其补充材料里的结构。在S66文中通过CCSD(T)/CBS计算出来的结果是-2.40 kcal/mol,可见本文的结果和金标准符合得非常好。

为了让结合能计算更省事,我还专门写了个脚本,是本文文件包里的getEbind.sh。此脚本在运行后会处理当前目录下的所有out文件,自动提取整体和两个片段的能量,并求差给出不同单位的结合能。例如当前目录下有三个.out文件,都是上文的方法用Multiwfn产生的任务的输出文件,getEbind.sh输出的信息如下,可见提取结合能非常方便。
H2O_ext.out
Etot = -761.79953946  E1 = -685.35948656  E2 = -76.43844639 Hartree
-0.0016065 Hartree
    -1.008 kcal/mol
    -4.218 kJ/mol


H2O_int.out
Etot = -761.80467188  E1 = -685.35943188  E2 = -76.43928426 Hartree
-0.0059557 Hartree
    -3.737 kcal/mol
   -15.637 kJ/mol


H2.out
Etot = -686.52275741  E1 = -685.35832899  E2 = -1.16228297 Hartree
-0.0021455 Hartree
    -1.346 kcal/mol
    -5.633 kJ/mol



2022-Jan-9补充

从2022-Jan-9更新的Multiwfn开始,在产生ORCA输入文件的界面里选择任务类型时可以选择-7,用法同前。此时产生的ORCA输入文件会计算5个能量,依次为:整体的能量、整体基函数下片段1的能量、整体基函数下片段2的能量、片段1的能量、片段2的能量。执行本文文件包里的ORCA_CP.sh就可以对当前目录下所有这种任务的输出文件进行处理,此脚本对每个文件会输出这些信息:(1)5个能量的具体数值 (2)原始的相互作用能 (3)BSSE校正后的相互作用能 (4)BSSE校正能 (5)BSSE校正后的整体的能量。这种任务所做的计算、脚本给出的信息,和Gaussian程序的counterpoise关键词完全相同。

本文文件包里的waterdimer.out是一个例子文件,ORCA_CP.sh对它进行处理的输出如下:

File: waterdimer.out
E(AB) = -152.054002577 Hartree
E(A)_AB = -76.022655799    E(A) = -76.022397044 Hartree
E(B)_AB = -76.024059371    E(B) = -76.022540157 Hartree
Raw interaction energy:
-0.0090654 Hartree
    -5.689 kcal/mol
   -23.801 kJ/mol
BSSE corrected interaction energy:
-0.0072874 Hartree
    -4.573 kcal/mol
   -19.133 kJ/mol
BSSE correction energy:
0.0017780 Hartree
     1.116 kcal/mol
     4.668 kJ/mol
BSSE corrected complex energy:
-152.0522246 Hartree

评分 Rate

参与人数
Participants 4
eV +20 收起 理由
Reason
corei70715 + 5 赞!
wuy069 + 5 谢谢分享
ABetaCarw + 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!

4万

帖子

99

威望

4万

eV
积分
89946

管理员

公社社长+计算化学玩家

 楼主 Author| 发表于 Post on 2020-5-23 06:05:51 | 显示全部楼层 Show all
在本文末尾增加了自动批量计算和输出结合能的脚本
北京科音自然科学研究中心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!

23

帖子

0

威望

2713

eV
积分
2736

Level 5 (御坂)

发表于 Post on 2020-5-23 12:47:31 | 显示全部楼层 Show all
sobereva 发表于 2020-5-23 06:05
在本文末尾增加了自动批量计算和输出结合能的脚本

社长,把这个功能加一个可以定义计算多个片段的功能

本版积分规则 Credits rule

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

GMT+8, 2023-2-6 05:26 , Processed in 0.201167 second(s), 25 queries .

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