顯著改進SMD溶劑模型描述Br和I元素精度的溶劑模型SMD18的介紹
顯著改進SMD溶劑模型描述Br和I元素精度的溶劑模型SMD18的介紹
文/Sobereva@北京科音 2021-Aug-23
之前筆者在《談談隱式溶劑模型下溶解自由能和體系自由能的計算》（http://www.shanxitv.org/327）中專門說過SMD溶劑模型，這是現階段最適合算溶劑環境下體系能量的隱式溶劑模型。SMD模型需要通過原子球疊加來構建溶質的孔洞用于算溶劑-溶質相互作用的靜電部分，定義原子球需要原子半徑。在SMD原文里，僅對H、C、N、O、F、Si、P、S、Cl、Br專門優化了半徑，對其它元素用的是Bondi范德華半徑來湊合（如果Bondi范德華半徑都沒有，就用2.0埃湊合）。
SMD溶劑模型作者后來的一篇文章Chem. Eur. J., 24, 15983 (2018)是SMD18溶劑模型的原文。這篇文章根據涉及鹵素的體系在乙腈中的反應自由能，發現對于碘用Bondi范德華半徑1.98埃來湊合算的結果很不好，此半徑值明顯太小，于是改為了1991年SM1溶劑模型里對碘優化的半徑2.74埃，這令計算結果與實驗數據相符程度有了極大的提升。后來他們對Br的半徑進行了專門優化，發現用2.60埃做半徑的結果與實驗數據最接近，這比SMD原文里專門給出的Br的半徑3.06埃要小不少，這也令結果大有改進。將SMD里Br和I的半徑分別改為2.60埃和2.74埃之后，就被此文稱為SMD18溶劑模型。
對于溶液中的涉及鹵鍵的形成或斷裂的過程，與Br、I關系密切的化學過程（例如成鍵斷鍵），或者算含有Br和I的體系的溶解自由能，都強烈建議使用SMD18代替SMD模型。
根據文中的測試，SMD18模型非常適合結合M06-2X使用，我也建議用SMD18時都用M06-2X泛函。為了最大程度地誤差抵消，計算溶液中自由能（1M標準態濃度）用的流程也最好和原文一致，具體來說是：對溶質在M06-2X下做優化和振動分析，在這個過程就帶著SMD18溶劑模型，直接取量子化學程序輸出的自由能，然后再加上1.89 kcal/mol的1atm->1M標準態轉變自由能變即可。對Cl、Br、I應當用def2-TZVPD（此時對I來說是贗勢基組，需要結合相應的贗勢），對其它元素用def2-TZVP。如果不清楚怎么在Gaussian中這么用基組，看《給ahlrichs的def2系列基組加彌散的方法》（http://www.shanxitv.org/340）和《詳解Gaussian中混合基組、自定義基組和贗勢基組的輸入》（http://www.shanxitv.org/60）。按照SMD18原文所述，如果存在100cm^-1波數以下振動頻率的情況，應當改用Grimme的準RRHO模型計算熱力學量而非用Gaussian默認用的RRHO模型，這可以通過Shermo程序實現，詳見《使用Shermo結合量子化學程序方便地計算分子的各種熱力學數據》（http://www.shanxitv.org/552），后文有例子。
PS：當然了，用def2-TZVP、def2-TZVPD做opt freq很昂貴，顯然不是可取的做法，但無奈SMD18原文這么干的，為了和SMD18優化的半徑盡可能實現誤差抵消，自己計算時也只能忍著了。如果實在算著太吃力的話，也可以opt freq時候用便宜一些的基組比如def2-TZVP(-f)（這么做對于自由能熱校正量和結構的影響很小），之后再用前述SMD18原文里用的基組算個單點，之后再把電子能量和自由能熱校正量相加。另外，原理上考慮ZPE校正因子更好，但也為了與SMD18原文一致，就不去考慮這點了。
目前的Gaussian不支持SMD18溶劑模型，但只要在SMD溶劑模型計算時改一下Br和I的元素半徑即可輕易實現。這需要寫scrf(SMD,solvent=溶劑名,read)關鍵詞，輸入文件末尾空一行寫modifysph，下面再空一行寫以下內容即可
I 2.74
Br 2.60
下面給一個Gaussian的完整的計算例子輸入文件。輸入文件和Gaussian 16跑的輸出文件可以在http://www.shanxitv.org/attach/613/file.zip里得到。
#P M062X/genecp SCRF(SMD,solvent=acetonitrile,read) opt freq
[空行]
Title
[空行]
0 1
N 2.53946200 -1.07717200 -0.01319200
C 1.74543500 0.00000700 -0.00027500
N 2.53948000 1.07717200 0.01885100
C 3.85512000 -0.67530300 -0.00660900
H 4.66037300 -1.38842900 -0.01940300
C 3.85513000 0.67528200 0.00619900
H 4.66038900 1.38840200 0.01901300
I -0.36464300 0.00000000 0.00004200
C 2.13385200 2.47671900 -0.00046400
H 1.05119700 2.53713200 0.26086500
H 2.74614800 3.02347800 0.75566000
H 2.31269200 2.89578200 -1.02444600
C 2.13383000 -2.47671400 -0.00044100
H 1.05117500 -2.53711700 -0.25836000
H 2.74216300 -3.02400000 -0.75519100
H 2.31661900 -2.89526100 1.02491500
Br -3.43164200 0.00000200 0.00045300
[空行]
@/sob/SMD18.gbs/N
[空行]
@/sob/def2-ECP.txt/N
[空行]
modifysph
[空行]
I 2.74
Br 2.60
[空行]
[空行]
被引入的SMD18.gbs是BSE基組數據庫的def2-TZVP基組文件，但是把Cl、Br、I替換成了def2-TZVPD之后的情況。def2-ECP.txt是def2系列基組標配的Stuttgart贗勢基組的定義。即便你的體系里沒有I，用上面的輸入文件直接替換坐標部分后也可以算，只不過讀取贗勢的設定以及修改I的半徑不生效而已。這兩個文件可以在http://www.shanxitv.org/attach/613/SMD18.gbs和http://www.shanxitv.org/attach/613/def2-ECP.txt下載，上例放到了/sob目錄下。
對于Gaussian 09，特別要注意對元素半徑的修改沒法在opt freq任務中傳遞給freq部分，因此必須opt和freq分別做，這是個bug，在Gaussian 16中得到了修正。
上面的例子的最低頻率只有44.24 cm^-1，因此不建議直接讀取輸出文件里的基于RRHO模型算出來的Sum of electronic and thermal Free Energies，否則低頻對熵的貢獻很不準確，而應當用Shermo程序基于準RRHO模型算。把Shermo的settings.ini里的ilowfreq=后面的值設為2，啟動Shermo并載入上例的輸出文件，可看到結果Sum of electronic energy and thermal correction to G: -3176.6090365 a.u.。將之加上標準態轉換自由能變1.89/627.51=0.00301 a.u.后即得到1M濃度下的乙腈中的自由能-3176.60602 a.u.。
SMD18模型體現出隱式溶劑模型計算能量的準確度對原子半徑是多么的敏感。順帶一提，J. Phys. Chem. A, 123, 9498 (2019)提出的SMD-B模型將SMD定義的原子半徑都改為了Bondi范德華半徑，發現對離子體系（特別是含硫的）的溶解自由能得到一定改進。既可以如上通過modifysph來對各個元素一一修改半徑，也可以在scrf里寫read的同時在輸入文件末尾后面空一行直接加上Radii=Bondi來實現，當然后者省事得多。