結合自由能計算：MD vs. QM
我在Gromacs結合自由能計算教程（一）中提到過，結合自由能的計算可以使用基於量子化學的方法，也可以使用基於分子動力學模擬（MD）的方法，然而兩者的原理大為不同。量子化學（QM）的方法考慮了實際電子結構，原則上說在計算相互作用焓的時候應該比做了大量近似處理的分子力學精確的多，但缺陷在於其巨大的計算成本；只能對於某一個或者幾個期待的構象進行靜態的研究，難以對構象變化帶來的熵進行計算——量化中熵的計算僅能考慮振轉熵貢獻，是基於hessian矩陣的計算。
相比之下分子力學能夠對體系的構象變化進行充分的採樣，並且可以處理大如蛋白的體系，使用顯式溶劑模型——顯式模型下才能考慮到某些複雜的能量變化，例如疏水空腔的「高能水」釋放，離散溶質效應……
相信很多人和我一樣一直抱有疑問，既然MD在底層計算上精度不夠，QM又難以進行充分的統計學採樣，二者在實際計算中，到底那種更準確？
這裡就要講到一篇文獻Hyrophobic challenge, 文中比較了幾種QM方法和兩種MD方法對於預測疏水中性客體與葫蘆脲的結合自由能與實驗值的差值。
文中所使用的主客體如下圖，其中DAPI為實驗測量鍵合常數所用的熒光探針
對於MD部分，使用牽拉法（傘型採樣）計算自由能
MD1：使用OPC水模型，從R.E.D網站獲得的RESP電荷結合GAFF2力場生成拓撲。
MD2：使用TIP3P水模型，使用AM1-BCC電荷結合GAFF力場生成拓撲。
原則上MD1的精度應該是要高於MD2的，無論是水模型還是電荷、力場大概都要優於MD2。
對於QM部分，自由能為電子能（內能）、振轉熵貢獻、溶劑化自由能三步獨立計算之和。也許是為了比較與MD的水和能計算差異，QM是在真空下基於真空優化的結構做單點和振轉分析的，然後再單獨計算水和能的貢獻（除了QM4考慮了溶劑對構象的影響）。這種做法我個人並不推薦在實踐中這麼做，真空中構象和溶劑中的構象往往會有很大差異，尤其是對於帶電分子這麼做就真是很不明智了；正如之後結果也會提到的，QM4因為考慮了溶劑對構象的影響，實際結果更為準確。如果條件允許，構象優化與結合能、振轉分析都在SMD水溶劑模型下做，不用特意拆開來去算溶劑化自由能，可能會更準確。
幾何構型優化與振轉熵貢獻均在PBEh-3c下完成。作者採用了四種不同的方法計算結合能：
- QM1：色散矯正meta-GGA TPSS-D3ATM/def2-QZVP，水和能在COSMO-RS 12fine模型下
- QM2：色散矯正meta-hybrid PW6B95-D3ATM/def2-QZVP（對於重原子捨去f軌道，對於氫原子捨去g軌道），水和能在COSMO-RS 13模型下
- QM3：色散矯正DLPNO-CCSD(T)/CBS*，水和能在COSMO-RS 13模型下
- QM4：色散矯正DLPNO-CCSD(T)/CBS*，為了考慮溶劑對最優構象的影響，使用XTB+GBSA預優化，然後在PBEh-3c+DCOSMO-RS進一步優化的結構下做單點，水和能在COSMO-RS 13模型下
特別地，對於QM1, QM2，為了改善QM方法對構象採樣不足的弱點，還使用了QMDFF（Stefan自己開發的基於量化計算hessian矩陣生成的分子力場）在優化構象進行模擬退火，得到了幾種能量最低構象，然後在PBEh-3c級別再優化， 一併進行相同級別QM計算，最後進行玻爾茲曼加權，以考慮構象熵貢獻。對於QM3、QM4，由於量化方法的高精度耗時極高，只考慮單個構象能量。
最後的結果如下，縱軸代表理論計算結果，橫軸為實驗結果，理想狀態下所有點都應該落在對角線上。
總體而言，幾種QM方法計算的結合自由能與實驗數據有著更小的誤差，但是和實驗數據的相關性較差；而MD方法總體與實驗數據有著較大的誤差，對於結合能大的體系誤差更大，但與實驗數據有良好的相關性。
這個結果可以說在意料之中，它們的誤差是由其固有的缺陷造成的——如前面所述，MD的底層理論簡陋，QM無法進行充分採樣。
MD1與MD2的差距十分有限，這也可以理解，下面會談到溶劑化能其實各種方法都算的挺準的，而GAFF2對GAFF的改進有限，AM1-BCC雖然略顯粗糙但也有不錯的準確度，綜合下來看就是預測上基本沒什麼不同。
再來看QM的數據，幾種方法共同的一點是，對於簡單碳氫化合物1~5預測基本都在回歸線上，但對於19、20、21、17，前兩者為六元環，後兩者為較長鏈烷烴，此時計算值就與回歸線有較大偏差 （同為長鏈的26和六元環的25誤差也較大，但只在QM4中出現），這是可能是因為在這幾種分子實際結合時的構象較為複雜（苯環看似簡單，實則存在經式緯式構象，並且容易旋轉），量化方法難以充分採樣的缺點在此暴露了出來。另外一點可以作證這個觀點的證據是，採用了構象加權矯正的QM1、QM2，雖然與實驗數據偏差較大，但與實驗數據的相關性好了一些，這和MD的情況類似。
至於水和能的計算方面，COSMO-RS 13（紅線）與MD（藍線）方法都與實驗吻合得很好，看來主要的誤差是在主客體結合能的計算與熵貢獻方面。
總體看來，這一輪較量之中MD是處於下風。但MD就真的不適合計算自由能了嗎？
其實不然，首先要考慮到在這組研究中所使用的主體分子是葫蘆脲，一種較為剛性的分子，它本身的構象變化較少，因此本身比較利於量化計算；並且例子中太多分子都是結構比較簡單的分子，但即便如此還是有個別分子存在構象問題造成的偏差。對於實踐中常常處理的蛋白質空腔以及其客體分子可能會有更大的構象複雜度，這時候量化方法的偏差可能會更大。
另一方面，MD的好處還是在於其較小的計算開銷，能夠處理顯式溶劑、大體系（蛋白質）、長時間採樣，在顯卡加速和多核計算加持下計算速度遠遠優於高精度QM，尤其是對於蛋白質體系，QM幾乎是無法處理的，雖然也有多精度組合方法的QM，但是操作過於複雜，怎樣劃分不同精度處理的界限也是難題。所以我個人認為，對於蛋白質與配體的結合自由能處理應該首先考慮MD方法。
文中還有更多有意思的內容，對此感興趣的請務必讀一下原文：
https://pubs.acs.org.ccindex.cn/doi/full/10.1021/acs.jpcb.7b09175#t1fn5推薦閱讀：
