Gaussian教程 | 氫鍵與鹵鍵的計算以及BSSE校正

關鍵字:氫鍵 鹵鍵 相互作用能 基組重疊誤差(BSSE) CP校正

一.基本概念

1. 氫鍵(Hydrogen Bond)

當氫原子與一個電負性較大的原子以共價作用結合時,氫原子會表現出較強的正電性,此時,當另一個電負性較大的原子(比如O,N,F原子)靠近這個氫原子時,這兩個原子之間產生的相互作用就稱為氫鍵。氫鍵具有飽和性與方向性。

2. 鹵鍵(Halogen Bonds)

當一個鹵素原子與其他原子以共價作用結合時,由于鹵素原子周圍靜電勢分布的差異,鹵原子本身會表現出雙親性,也就是在垂直于共價鍵的方向,鹵素原子表現出親核性;而在沿著共價鍵軸的方向表現出親電性,這個表現出親電性的區域也稱為“σ-hole”,當“σ-hole”與作為路易斯堿的原子(比如N,O,S)接近時,產生的相互作用就稱為鹵鍵。

3. 相互作用能(Interaction Energy)

從廣義上來說,當任何兩個分子靠近時,假若,形成的復合物的電子能低于兩個單獨分子的電子能加和,那么復合物電子能與兩個單獨分子電子能和之差就可以定義為相互作用能。如果兩個分子主要以氫鍵作用結合,那么就稱為氫鍵相互作用能;同理,如果兩個分子主要以鹵鍵作用結合,就稱為鹵鍵相互作用能。A和B兩個分子間的相互作用能的計算可以統一寫成如下形式:

E(interaction) = E(AB) – E(A) – E(B)

其中,E(interaction)是分子間的相互作用能,E(AB)是分子A和B形成的復合物的能量,E(A)和E(B)是分子A和B獨立存在時的能量。

4. 基組重疊誤差(Basis Set Superposition Error, 簡稱BSSE)

上式從理論上給出了相互作用能的計算公式,可以看出,相互作用能的準確性直接依賴于復合物和單個分子電子能計算的準確性。然而,在實際計算過程中,當計算復合物的電子能時,分子A和B各自的基組會彌散到對方的分子區域,相當于增大了各自的基組,導致計算的復合物能量偏低,最終會使計算的相互作用能比實際更高,這也被稱為基組重疊誤差(BSSE)。

5. 均衡校正(Counterpoise Correction,也稱CP校正法)

為了消除BSSE對相互作用能計算準確性的影響,我們可以采用CP校正法來計算BSSE產生的能量誤差,以得到更準確的相互作用能。

二.利用Gaussian程序計算氫鍵,鹵鍵能并進行BSSE計算(CP校正)的具體步驟

氫鍵和鹵鍵屬于弱相互作用,尤其是當研究的體系中氫鍵或鹵鍵還相對較弱時,此時BSSE帶來的能量誤差很可能與弱相互作用能本身持平,因此若不進行BSSE校正,得到的結果便會很不可靠,Gaussian提供了Counterpoise關鍵字來進行BSSE能量的計算。需要注意的是,BSSE的大小與基組大小密切相關,當基組增大時,BSSE減小。

在本教程中我們將計算HF分子與NH3分子的氫鍵能,以及IF分子與NH3分子的鹵鍵能(如圖1),分別使用PM6,HF/6-31G(d),MP2/6-31G(d),B3LYP/6-31G(d)和B3LYP/6-311++G(d,p)五種方法,并比較BSSE校正前后相互作用能的差異。

氫鍵與鹵鍵相互作用算例

圖1. 算例:氫鍵相互作用(左)與鹵鍵相互作用(右)

計算弱相互作用能的基本流程(包含BSSE校正):

  1. 首先分別優化兩個獨立分子和復合物分子的結構。
  2. 使用優化好的分子作為初始結構,用更大的方法基組對兩個獨立的分子進行單點能計算,并對復合物分子進行BSSE計算。
  3. 結果分析

在這個教程中,為了比較不同方法計算的相互作用能,我們首先使用B3LYP泛函優化3個獨立分子(HF,IF和NH3),以及2個復合物分子(氫鍵復合物和鹵鍵復合物)的結構,這些結構將被使用在后續的能量計算中。作為演示,我們僅給出氫鍵相互作用能計算過程中所涉及的輸入文件,鹵鍵相互作用能的計算過程完全相同,因此只給出最后結果。

1. 計算HF和NH3分子的氫鍵能

第一步:優化HF,NH3和氫鍵復合物的分子結構。

HF分子結構優化的輸入文件:

%mem=1GB
%nprocshared=8
# b3lyp/gen opt freq 

HF

0 1
 F                 -2.14398099    0.55747491    0.69389897
 H                 -1.26398099    0.55747491    0.69389897

H F 0
6-31G*
****

!最后一行為空白行,不可省略

NH3分子結構優化的輸入文件:

%mem=1GB
%nprocshared=8
# b3lyp/gen opt freq 

NH3

0 1
 N                 -0.02122157    0.25666598    0.95462239
 H                  0.31210032   -0.68614710    0.95462239
 H                  0.31211753    0.72806616    1.77111912
 H                  0.31211753    0.72806616    0.13812565

N H 0
6-31G*
****

!最后一行為空白行,不可省略

氫鍵復合物結構優化的輸入文件:

%mem=1GB
%nprocshared=8
# b3lyp/gen opt freq 

NH3-HF

0 1
 N                 -0.02122157    0.25666598    0.95462239
 H                  0.31210032   -0.68614710    0.95462239
 H                  0.31211753    0.72806616    1.77111912
 H                  0.31211753    0.72806616    0.13812565
 F                 -2.14398099    0.55747491    0.69389897
 H                 -1.26398099    0.55747491    0.69389897

N H F 0
6-31G*
****

!最后一行為空白行,不可省略

第二步,用MP2/6-31G(d)方法基組計算HF和NH3分子的單點能以及氫鍵復合物的BSSE

使用第一步優化得到的結構作為初始結構,用MP2/6-31G(d)方法基組計算HF和NH3分子的單點能,以及氫鍵復合物的BSSE。

HF分子單點能計算的輸入文件:

%mem=1GB
%nprocshared=8
# MP2/gen

HF

0 1
 F                  0.00000000    0.00000000    0.09348300
 H                  0.00000000    0.00000000   -0.84134400

H F 0
6-31G(d)
****

!最后一行為空白行,不可省略

NH3分子單點能計算的輸入文件:

%mem=1GB
%nprocshared=8
# MP2/gen

NH3

0 1
 N                  0.00000000    0.00000000    0.11863500
 H                  0.00000000    0.93931000   -0.27681500
 H                 -0.81346600   -0.46965500   -0.27681500
 H                  0.81346600   -0.46965500   -0.27681500

H N 0
6-31G(d)
****

!最后一行為空白行,不可省略

氫鍵復合物BSSE計算的輸入文件:

%mem=1GB
%nprocshared=8
# MP2/gen Counterpoise=2

NH3-HF

0 1 0 1 0 1 
 N(Fragment=1)                 -1.23505300    0.00017600   -0.00001900
 H(Fragment=1)                 -1.61000000   -0.81197600   -0.48735400
 H(Fragment=1)                 -1.61403400    0.82687900   -0.45893100
 H(Fragment=1)                 -1.61219900   -0.01716500    0.94608500
 F(Fragment=2)                  1.44453500    0.00007800   -0.00004800
 H(Fragment=2)                  0.48078400    0.00033300    0.00076900

H F N 0
6-31G(d)
****

!最后一行為空白行,不可省略

其中第3行的Counterpoise=2指示Gaussian程序進行BSSE能量計算,數值等于2表示只計算兩個分子片段間的BSSE能量。

其中第7行表示在計算兩個片段間BSSE能量時需要設置三組電荷和自旋多重度,書寫順序分別為復合物整體、分子片段1與分子片段2。在本算例中,三組電荷和自旋多重度都為0 1,所以第7行為0 1 0 1 0 1。

第8-13行為分子結構描述部分,在BSSE計算中需要為每個原子指定其所屬的片段。由于我們計算兩個片段間BSSE能量,所以需要在分子結構描述部分區分這兩個片段。若某個元素屬于片段1,就在這個元素符號后面寫上(Fragment=1),至于第一個片段設置為1或2并不固定,但必須與前面電荷、自旋多重度保持一致。

第三步:利用GaussView查看CP校正后的分子能量

如圖2所示,用GaussView打開BSSE計算的輸出文件,然后在Results下拉菜單中點擊Summary選項。

查看BSSE校正后的復合物能量

圖2. 查看BSSE校正后的復合物能量

BSSE校正的復合物能量還可以從Gaussian的結果文件里獲取,其特征是含有Counterpoise corrected energy字樣。如下第1行(單位為Hatree)。

 Counterpoise corrected energy =    -156.551181178325
                   BSSE energy =       0.003661058666
              sum of fragments =    -156.530582019229
           complexation energy =     -15.22 kcal/mole (raw)
           complexation energy =     -12.93 kcal/mole (corrected)

我們可以直接用上述BSSE校正后的復合物能量減去兩個獨立的分子能量和,這個數值就是最終的經過BSSE校正的兩個分子的相互作用能(在這里主要是氫鍵能)。用上述方法同樣可以查看兩個獨立分子的單點計算電子能,在這里不做贅述。

根據我們上述的討論,校正后的相互作用能可以寫成如下形式:

E(interaction) = E(AB) + E(BSSE) – E(A) – E(B)

其中E(AB) + E(BSSE)就等于上面得到的BSSE校正后的復合物能量。

NH3計算結果文件最后幾行如下:

 1\1\GINC-C01N01\SP\RMP2-FC\Gen\H3N1\YLIANG2\22-Nov-2019\0\\# MP2/gen\\
 NH3\\0,1\N,0,0.,0.,0.118635\H,0,0.,0.93931,-0.276815\H,0,-0.813466,-0.
 469655,-0.276815\H,0,0.813466,-0.469655,-0.276815\\Version=ES64L-G16Re
 vA.03\State=1-A1\HF=-56.1830004\MP2=-56.3519668\RMSD=2.550e-09\PG=C03V
  [C3(N1),3SGV(H1)]\\@

顯而易見地,優化過后的能量E(MP2)=-56.3519668(Hatree)。

HF計算結果文件最后幾行如下:

 1\1\GINC-C01N01\SP\RMP2-FC\Gen\F1H1\YLIANG2\22-Nov-2019\0\\# MP2/gen\\
 HF\\0,1\F,0,0.,0.,0.093483\H,0,0.,0.,-0.841344\\Version=ES64L-G16RevA.
 03\State=1-SG\HF=-100.0001739\MP2=-100.179425\RMSD=9.310e-09\PG=C*V [C
 *(H1F1)]\\@

則HF優化過后的能量E(MP2)=-100.179425(Hatree)。

對于這個例子,氫鍵相互作用能為:

E(interaction) = E(complex_Corrected) – E(HF) –E(NH3)
= -156.551181178325 – (-100.179425) – (-56.3519668)
= -12.4 kcal/mol


需要注意的是:BSSE計算的輸出文件也會給出一個相互作用能的數值,然而這個數值是復合物能量減去以復合物中單體結構為初始構型計算得到的單點能。而復合物中的單體構型與獨立優化的單體結構是存在差異的,因此,我們需要獨立優化的單體結構作為初始構型進行單點計算,以進行最終的相互作用能的計算。

在這個教程中我們嘗試了5種不同的方法基組組合,分別為PM6;HF/6-31G(d);MP2/6-31G(d);B3LYP/6-31G(d);B3LYP/6-311++G(d,p)。對每一個水平的方法基組重復上述第二步和第三步的過程,我們將得到不同計算水平下的氫鍵能。

需要注意的是,在計算鹵鍵能時,對于碘元素,我們使用贗勢SDD,其余同上,關于基組的使用可以參見教程:《Gaussian教程 | 使用基組和贗勢》

三. 結果分析

在表1中,我們列出了不同方法基組下氫鍵能計算的一些結果:包括BSSE能量和校正后的氫鍵相互作用能。

表1. 不同方法基組下氫鍵能的計算結果(單位kcal/mol)

PM6 HF/6-31G(d) MP2/6-31G(d) B3LYP/6-31G(d) B3LYP/6-311++G(d,p)
BSSE能量 132.5 1.3 2.3 1.8 1.1
校正后的相互作用能 128 -11.0 -12.4 -13.7 -13.4


在表2中,我們列出了不同方法基組下鹵鍵能計算的一些結果:包括BSSE能量和校正后的鹵鍵能。

表2. 不同方法基組下鹵鍵能的計算結果(單位:kcal/mol)

PM6 HF/6-31G(d) MP2/6-31G(d) B3LYP/6-31G(d) B3LYP/6-311++G(d,p)
BSSE能量 190.8 1.3 2.8 2.0 0.8
校正后的相互作用能 168.7 -12.1 -14.1 -17.2 -17.7


從上面的結果可以看出,半經驗方法PM6給出的結果完全不可用,通過B3LYP/6-31G(d)和B3LYP/6-311++G(d,p)計算結果的比較,說明基組的增大能夠降低BSSE能量,這也一致于前人的結論。值得注意的是,盡管HF/6-31G(d)水平也給出了較小的BSSE能量,但不表示這個水平下相互作用能的計算也是準確的。至于在具體的研究中,應該選擇什么樣的方法基組還要視研究體系而定,并多多調研文獻。

四. 注意事項

本文全部的計算在Gaussian 16下完成;另外,本文的算例是在氣相真空中進行計算,實際上溶劑也會對分子間的相互作用能產生影響,因此在研究過程中需要加以考慮。

五. 相關主題

  1. 復合物結構的準備(預測)
  2. CONFLEX教程 | 主客體構象搜索教程》演示了如何用CONFLEX預測兩個或多個分子間的主客體復合物結構。

六. 參考資料

  1. Exploring Chemistry with Electronic Structure Methods, 3rd ed., Gaussian, Inc.: Wallingford, CT, 2015. J. B. Foresman and ? Frisch

七. 聯系我們

購買Gaussian,請聯系我們

奇米777四色影视在线看_国内一点不卡在线播放视频_比较有韵味的熟妇无码