在ORCA中做counterpoise校正并計算分子間結合能的例子
在ORCA中做counterpoise校正并計算分子間結合能的例子
文/Sobereva@北京科音
First release: 2020-Mar-13 Last update: 2022-Jan-9
在《談談BSSE校正與Gaussian對它的處理》（http://www.shanxitv.org/46）一文中介紹過BSSE問題與counterpoise校正。Counterpoise校正在高精度弱相互作用計算時通常需要考慮。Gaussian有現成的counterpoise關鍵詞，因此做counterpoise校正很簡單。但是ORCA程序（至少對于撰文時最新的4.2.1版來說）沒有提供現成的counterpoise關鍵詞，只能通過定義鬼原子的方式手工實現，因此實現起來稍微麻煩。不熟悉ORCA者看http://www.shanxitv.org右側ORCA分類里我寫過的大量相關文章。
之前在《詳談Multiwfn產生ORCA量子化學程序的輸入文件的功能》（http://www.shanxitv.org/490）中筆者專門介紹了Multiwfn產生ORCA量子化學程序的輸入文件的功能。為了使得用ORCA做counterpoise任務更方便，筆者2020-Mar-13給Multiwfn的這個功能增加了相應選項，在這里通過一個實例簡單演示一下。最新版Multiwfn可以在http://www.shanxitv.org/multiwfn免費下載，相關知識見《Multiwfn FAQ》（http://www.shanxitv.org/452）。
注意按照下文例子算出來的實際上是復合物結構中單體間相互作用能。而一般意義的結合能在計算時是要考慮單體形成復合物過程中形變造成的能量改變（變形能）。由于下面的例子里單體是剛性的，變形能可以忽略不計，因此下文仍稱作結合能。
本文用的示例體系是下圖這個體系，我們要在DLPNO-CCSD(T)/aug-cc-pVTZ下計算其結合能，并考慮counterpoise校正。
此體系的xyz格式的結構文件，后文中涉及的輸入輸出文件以及要用的shell腳本都可以在這里下載：http://www.shanxitv.org/attach/542/file.rar。并不是必須得用xyz格式，只要是Multiwfn支持的含有結構信息的格式（如pdb、mol2...）都可以作為輸入文件，見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/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
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