• 比SMD算溶解自由能更好的溶劑模型uESE的使用

    比SMD算溶解自由能更好的溶劑模型uESE的使用

    文/Sobereva@北京科音 2021-Apr-30


    1 uESE溶劑模型簡介

    SMD是如今非常常用的溶劑模型，在《談談隱式溶劑模型下溶解自由能和體系自由能的計算》（http://www.shanxitv.org/327）里筆者專門介紹過。在Gaussian程序中，SMD的極性部分等同于IEFPCM，而由于其非極性部分考慮較周到，因此溶解自由能計算精度較好。但是SMD算離子體系的溶解自由能誤差遠大于算中性體系的，通常需要用顯式+隱式的雜化溶劑模型才能得到較好結果。

    最近Vyboishchikov等人提出了uESE (universal Easy Solvation Energy Evaluation)溶劑模型，原文見J. Comput. Chem. (2021) DOI: 10.1002/jcc.26531。之前還有xESE溶劑模型，和uESE形式很相似，原文見Phys. Chem. Chem. Phys., 22, 14591 (2020)。uESE/xESE都像SMD一樣考慮了溶質-溶劑相互作用的極性和非極性部分，但uESE/xESE和SMD有很大不同：
    (1)uESE/xESE的極性部分利用的是COSMO模型。通常COSMO是基于電子密度計算的，而是uESE/xESE的COSMO則是基于CM5原子電荷計算的；更具體來說，在求解表面顯著電荷時利用使用CM5電荷計算溶質孔洞表面的靜電勢。
    (2)uESE/xESE模型并不納入SCF迭代過程中，也不影響電子結構，不能計算帶溶劑模型時的受力和Hessian，相當于是一個“后”溶劑模型。uESE/xESE僅適合基于CM5電荷算溶解自由能的目的。由于uESE/xESE中的參數是基于氣相下B3LYP/def2-TZVP級別的結構和波函數計算的CM5原子電荷擬合的，所以自己用uESE/xESE的時候CM5電荷也得在這個級別下計算。
    (3)uESE/xESE包含對溶劑-溶質極性作用部分的經驗校正項。
    (4)uESE/xESE對于不同類型溶劑使用不同的極性部分的校正項和非極性部分形式。溶劑分為四類：水，非水極性質子溶劑，極性非質子溶劑，非極性溶劑。
    (5)uESE/xESE擬合的經驗參數中既有依賴于元素的，也有依賴于溶劑類型的，也有依賴于溶劑本身的。目前xESE只支持水，uESE支持約100種溶劑。極性部分校正項和非極性部分只有H,C,N,O,F,P,S,Cl,Br,I元素的參數（對于其它元素只能給出COSMO極性部分的結果）。而SMD不含依賴原子半徑以外元素的參數。

    uESE和xESE有什么存在意義？主要意義在于算溶解自由能從統計結果上比SMD更好，特別是對于離子體系。根據uESE原文的測試，溶解自由能計算精度（平均絕對誤差，MAE）有下面的關系
    中性體系：xESE≈SM12>=uESE>=SMD
    陽離子體系：uESE>xESE>SMD>=SM12
    陰離子體系：uESE>xESE≈SMD12>>SMD
    即曰，算中性體系溶解自由能，用xESE會比SMD占一點優勢。算陽離子用uESE比用SMD強得多。算陰離子更是絕對不能直接用SMD，不想用雜化溶劑模型的話至少也應當用uESE。根據測試，uESE算中性體系溶解自由能的MAE小于1 kcal/mol，算陰、陽離子溶解自由能的MAE差不多，都小于3 kcal/mol。而SMD算陰離子的MAE則達到8 kcal/mol左右。注：上面的SM12和SMD是同門的，但前者和uESE一樣也是基于CM5原子電荷進行計算，不過極性部分利用的是廣義Born方程形式。

    下圖左側是uESE原文里算大量中性和離子在非水質子溶劑中的溶解自由能（縱坐標）與實驗值（橫坐標）的對比，右圖是SMD算的情況，誤差超過4 kcal/mol的用紅色符號高亮。可見uESE算中性溶質的情況比SMD稍好，而計算離子的情況能好很多。

    下圖是計算在極性非質子溶劑中的溶解自由能的情況。可見此時SMD算離子的情況表現極差，特別是對于溶解自由能非常大的離子

    需要注意的是以上都是統計數據，大家也不要指望算每個中性體系xESE都比SMD好、算每個離子體系uESE都一定比SMD好。筆者寫此文的目的主要是提醒大家算離子體系的溶解自由能的時候，如果懶得用雜化溶劑模型，則至少要用uESE，切勿拿SMD湊合（除非有uESE不支持的元素）。


    2 uESE和xESE程序的使用

    uESE和xESE方法有同名的計算程序，可以在http://iqcc.udg.edu/~vybo/ESE/直接下載，Windows和Linux版可執行文件都提供了，不開源。這兩個程序用的輸入文件相同，uESE可以指定溶劑，而xESE只支持水溶劑。

    uESE和xESE都需要基于命令行使用。uESE在Windows下的使用格式是：uESE.exe [輸入文件路徑] -solvent [溶劑名]。可以用的溶劑見http://iqcc.udg.edu/~vybo/ESE/solvent-list.html。

    uESE/xESE的輸入文件里需要提供各原子的元素在周期表里的序號、原子坐標，以及原子的CM5電荷。筆者開發的Multiwfn波函數分析程序（主頁&免費下載：http://www.shanxitv.org/multiwfn）直接就能計算CM5電荷。為了令用戶能盡可能方便地產生uESE/xESE的輸入文件，從2021-Apr-30更新的Multiwfn開始只要將settings.ini里的uESEinp參數設為1，則Multiwfn計算完CM5電荷后就會問你是否產生uESE/xESE的輸入文件。使用uESE/xESE算溶解自由能用的輸入文件若是Multiwfn產生的，請注意在你發表的文章中提及，并按照Multiwfn啟動時的提示恰當引用Multiwfn。

    Multiwfn計算CM5電荷需要提供含有波函數的信息作為輸入文件，諸如.wfn/.wfx/.fch/.molden/.mwfn等等，產生方式見《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）。不了解Multiwfn的話參看《Multiwfn FAQ》（http://www.shanxitv.org/452）。


    3 使用uESE計算PhO-在乙腈中的溶解自由能實例

    下面就通過一個例子演示uESE程序的使用，是PhO-（苯酚的羥基的質子解離掉后的陰離子）在乙腈中的溶解自由能計算。此例所有相關文件可以在http://www.shanxitv.org/attach/593/file.rar下載。本例Gaussian使用G16 A.03版，Multiwfn是2021-Apr-30更新的3.8(dev)版，操作系統是Win 10 64bit。

    首先用Gaussian對PhO-進行優化。根據uESE原文，應當用B3LYP/def2-TZVP在真空下進行。輸入文件如下

    %chk=C:\PhO-_optfreq.chk
    #p B3LYP/def2TZVP opt freq
    [空行]
    Title Card Required
    [空行]
    -1 1
     C                  0.00000000    0.00000000   -1.80531390
     C                  0.00000000   -1.20821014   -1.10754617
     C                  0.00000000   -1.20820376    0.28716582
     C                  0.00000000    0.00000000    0.98475086
     C                 -0.00000000    1.20820376    0.28716582
     C                 -0.00000000    1.20821014   -1.10754617
     H                  0.00000000    0.00000000   -2.90492418
     H                  0.00000000   -2.16044049   -1.65754504
     H                  0.00000000   -2.16042940    0.83722286
     H                 -0.00000000    2.16042940    0.83722286
     H                 -0.00000000    2.16044049   -1.65754504
     O                  0.00000000    0.00000000    2.41475076

    計算完成后，檢查輸出文件確認沒有虛頻。然后用formchk將PhO-_optfreq.chk轉換為PhO-_optfreq.fch。formchk的使用在《詳談Multiwfn支持的輸入文件類型、產生方法以及相互轉換》（http://www.shanxitv.org/379）里說了。

    將Multiwfn目錄下的settings.ini里的uESEinp設為1，然后啟動Multiwfn，載入PhO-_optfreq.fch，之后依次輸入
    7  //布居分析與原子電荷計算
    16  //計算CM5電荷
    1  //使用內置的球對稱化的自由原子密度
    [按回車]  //產生的uESE輸入文件為當前目錄下的PhO-_optfreq.inp
    然后就可以把Multiwfn窗口關了。

    下載uESE的Windows版可執行文件uESE.exe到某處，把PhO-_optfreq.inp放到與此文件相同的目錄下，然后按Win+R鍵，輸入cmd進入命令行界面，通過cd命令進入此目錄下（這種計算機基本常識不會的話自行Google），運行uESE PhO-_optfreq.inp -solvent acetonitrile。僅需一秒鐘就能運行完，在輸出信息中會看到
     COSMO electrostatic energy =       -66.946 kcal/mol
     Correction term =        12.393 kcal/mol
     Total solvation free energy =       -54.554 kcal/mol

    即曰溶解自由能是-54.554 kcal/mol。uESE/xESE給出的溶解前后都是1M濃度的溶解自由能。從SMD溶劑模型原文的補充材料里可以查到實驗值是-55.1 kcal/mol（也是溶解前后都是1M的情況），可見此例uESE表現得極為理想，誤差僅有0.55 kcal/mol，不過這很大程度上是運氣好。

    下面也用SMD模型算一下溶解自由能，做法在《談談隱式溶劑模型下溶解自由能和體系自由能的計算》（http://www.shanxitv.org/327）里專門講過。基于上面B3LYP/def2-TZVP氣相下優化的結構，在M05-2X/6-31G*級別下算一次氣相單點能和SMD模型表現的乙腈環境下的單點能，結果為
    # M052X/6-31G* scrf(SMD,solvent=acetonitrile)：-306.941786445 Hartree
    # M052X/6-31G*：-306.847792107 Hartree
    因此溶解自由能為627.51*(-306.941786445+306.847792107)=-58.98 kcal/mol。相對于實驗的誤差達到-3.9 kcal/mol，明顯大于uESE的，這體現出了uESE的顯著優勢。

    注意uESE模型在測試、擬合參數的時候對離子體系都是只考慮了帶+1、-1電荷的離子。高價離子與溶劑相互作用顯著強于一價離子，筆者估計用uESE肯定不會得到較好結果，必須用雜化溶劑模型。

    久久精品国产99久久香蕉