通過SMD溶劑模型描述離子液體溶劑環境的方法
通過SMD溶劑模型描述離子液體溶劑環境的方法
文/Sobereva @北京科音 2018-Aug-7
離子液體已被廣泛作為溶劑。常有人問離子液體溶劑環境怎么在量子化學計算中體現，本文就專門說一下。閱讀本文之前務必先閱讀《談談隱式溶劑模型下溶解自由能和體系自由能的計算》（http://www.shanxitv.org/327）了解基礎知識。
1 原理
Truhlar搞出來的SMD溶劑模型是一種比較普適的隱式溶劑模型。原理上，只要你給了溶劑的以下參數，就可以通過SMD溶劑模型描述溶劑環境：
● 靜態介電常數eps(298K)
● 動態介電常數epsinf或折射率（前者是后者的平方）
● Abraham氫鍵酸度：來自Abraham的一些文章
● Abraham氫鍵堿度：來自Abraham的一些文章
● 氣液界面表面張力(298K)
● 芳香度：溶劑分子中芳香碳占重原子（即非氫原子）的比例
● 鹵素度：溶劑分子中鹵素占重原子（即非氫原子）的比例
對于SMD溶劑模型而言，離子液體作為溶劑并沒有任何特殊性，只要提供了離子液體的以上參數，就可以用SMD恰當地描述離子液體溶劑環境。這種做法的可靠性在SMD原作者的一篇文章J. Phys. Chem. B, 116, 9122 (2012)中得到了充分驗證，文中測試了一批常見的離子液體體系。從測試結果看，在計算溶解自由能的計算精度上，離子液體和普通分子作為溶劑的情況差不多，盡管對于有的離子液體誤差比較大，比如[EMIM][PF6]。
用SMD描述離子液體的一個困難之處是很多離子液體的SMD參數找不全。如果你感興趣的只是溶劑-溶質之間的極性部分作用還好，只需要定義eps，有些時候也需要定義epsinf（或折射率），這就夠了，而其它參數不會被利用；而如果還需要考慮溶劑的非極性部分貢獻，那就很麻煩，表面張力可能還不難找，但Abraham氫鍵酸度和堿度沒地方搞去。
為了解決離子液體SMD參數往往找不全的問題，JPCB那篇文章里提出一個模型叫SMD-GIL，其中GIL是generic ionic liquid的縮寫。這個模型并不是一個新的溶劑模型的形式，它的意思只不過是把除了芳香度和鹵素度以外的SMD參數都直接設為許多離子液體實驗上已知的參數的平均值。SMD-GIL的參數值在JPCB文中給了：
靜態介電常數=11.5
折射率=1.43（折合動態介電常數=1.43^2=2.0449）
Abraham氫鍵酸度=0.229
Abraham氫鍵堿度=0.265
表面張力=61.24 cal/mol/A^2
而離子液體的芳香度和鹵素度自己根據其化學組成手動一算就知道。
根據JPCB文中的測試，使用SMD-GIL參數計算和使用實際離子液體的實驗參數做SMD計算得到的結果精度半斤八兩，互有勝負。因此，通過SMD溶劑模型描述你用的離子液體，有實驗值的參數可以直接用實驗值，而沒有的就用SMD-GIL給出的離子液體的平均參數代替即可，結果是合理的；或者如果你懶得找實驗參數，則所有參數都用SMD-GIL的也完全沒問題。JPCB文中也給出了一些常見離子液體的實驗參數，見文中表1、3。
2 實例：計算乙醇在[BMIM][NTf2]中的溶解自由能
[BMIM][NTf2]的結構如下
乙醇在[BMIM][NTf2]中的常溫下的溶解自由能實驗值為-3.76 kcal/mol（來自那篇JPCB文章補充材料Table S1.6，對應氣相和溶劑下都是1M標準態濃度）。我們看看利用SMD-GIL的參數算出來的結果是多少。如《談談隱式溶劑模型下溶解自由能和體系自由能的計算》所述，我們先對結構優化，然后用M052X/6-31G*在溶劑模型下和氣相下分別算個單點，求差就是常溫下的溶解自由能（對應氣相和溶劑下都是1M標準態濃度）。
SMD幾乎是目前最流行的隱式溶劑模型，因此很多主流程序都已經支持SMD，比如Gaussian、NWChem、ORCA、GAMESS-US等（而非主流的Dmol3和ADF之流尚不支持）。這里我們用量化計算最常用的Gaussian作為例子，結構已經在B3LYP/TZVP下優化過（用什么級別優化無所謂，只要是合理級別即可，對溶解自由能計算結果影響很小）。
單點任務的輸入文件如下：
# M052X/6-31G*B3LYP/TZVP opted
0 1
C 1.17643900 -0.39950500 0.00000000
H 1.15236000 -1.03828400 0.88463100
H 2.11738900 0.15557800 0.00000000
H 1.15236000 -1.03828400 -0.88463100
C 0.00000000 0.55417200 0.00000000
H 0.04041900 1.20079700 -0.88644000
H 0.04041900 1.20079700 0.88644000
O -1.20073500 -0.22291400 0.00000000
H -1.95570500 0.37470800 0.00000000
溶劑模型下的單點任務的輸入文件如下。這里我們自己完全自定義了一個新溶劑，除了芳香度和鹵素度以外的溶劑參數就是上面提到的SMD-GIL的平均離子液體溶劑參數。
# M052X/6-31G* SCRF(SMD,read,solvent=generic)
B3LYP/TZVP opted
0 1
C 1.17643900 -0.39950500 0.00000000
H 1.15236000 -1.03828400 0.88463100
H 2.11738900 0.15557800 0.00000000
H 1.15236000 -1.03828400 -0.88463100
C 0.00000000 0.55417200 0.00000000
H 0.04041900 1.20079700 -0.88644000
H 0.04041900 1.20079700 0.88644000
O -1.20073500 -0.22291400 0.00000000
H -1.95570500 0.37470800 0.00000000
eps=11.5
epsinf=2.0449
HBondAcidity=0.229
HBondBasicity=0.265
SurfaceTensionAtInterface=61.24
CarbonAromaticity=0.12
ElectronegativeHalogenicity=0.24
上面文件中設的CarbonAromaticity就是芳香度，ElectronegativeHalogenicity是鹵素度。[BMIM][NTf2]的化學組成為C10H15F6N3O4S2，重原子數（非氫原子數）為25，此體系中咪唑環上的碳屬于芳香碳，數目為3，因此芳香度為3/25=0.12。體系中鹵素有6個，因此6/25=0.24就是鹵素度。
對兩個文件用Gaussian16 A.03版進行計算，氣相單點能為-154.9998251，溶劑下單點能為-155.0061949，因此溶解自由能為627.51*(-155.0061949+154.9998251)=-3.997 kcal/mol，和實驗值-3.76 kcal/mol相符不錯。