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

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

文/Sobereva@北京科音

First release: 2020-Mar-13  Last update: 2020-May-23


在《谈谈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校正。

此体系的xyz格式的结构文件,以及后文中涉及的输入输出文件都可以在这里下载: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,整个片段就被高亮,而且序号也显示出来了,如下所示。

如此这般,把两个片段序号范围都拷出来备用,然后点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

对于当前例子,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