背景

分子模擬的一大用途就是對體系的平衡態進行模擬採樣得到體系的某些熱力學數據，其中的自由能計算是分子模擬被應用的較多的一個場景之一；計算自由能的基本思想是研究系統在無限長時間內的能量均值等價於其系綜平均值，因而可以通過有限時間的取樣來估計系綜平均值，進而通過公式求解得到相關自由能數據。

雖然量子化學軟體（如Gaussian）也能夠處理自由能計算，其原理是把分子的振動，轉動，平動等對熵有貢獻的項，按照二次項展開，計算它的S與溫度的關係，從而可以求出一定溫度下的自由能變化。然而，對於動輒上萬原子的大體系如蛋白-配體體系而言，哪怕是半經驗的方法都極為耗時，更不要提及顯式溶劑的計算（一些與溶劑有強烈作用的情況下隱式溶劑模型就不盡人意了）。

自由能計算本質上是處理考察體系從一個熱力學態到另一個熱力學態發生的能量變化，這可以指其本身的變化（化合價變化）或者是周圍環境的變化（從真空到溶劑，或者是從溶劑到配位環境的變化）

因此自由能的計算任務可以解決如以下的問題：

1. 計算化合物在不同環境的分配係數（如計算立普妥在水和辛醇間的分配）；
2. 計算主客體的鍵合強度（如計算葫蘆[8]脲與碘單質鍵合常數） ——由此衍生出的計算有改變蛋白質或底物的結構以檢測它們鍵合的變化。

……

然而，一直以來分子模擬計算結合自由能所涉及到的統計熱力學知識與一系列複雜的模擬參數設置與繁瑣的數據分析，這一道道天然的門檻阻礙了這種方法在非MD專業研究者中的推廣。本教程旨在幫助缺乏MD經驗的人也能夠通過現有模板稍微修改參數，很快上手鍊金術自由能計算。

在這個教程中，我將關注於結合自由能的計算，以我的課題中處理的系統為例，逐步講解Gromacs的結合自由能計算實現方法。所涉及內容較多，篇幅大，將分成幾節講述。

在此之前強烈推薦閱讀以下幾篇參考鏈接，對理解本篇內容的原理很有幫助：

https://jerkwin.github.io/9999/10/31/GROMACS%E4%B8%AD%E6%96%87%E6%95%99%E7%A8%8B/#6-%E8%87%AA%E7%94%B1%E8%83%BD%E8%AE%A1%E7%AE%97-%E6%B0%B4%E4%B8%AD%E7%9A%84%E7%94%B2%E7%83%B7gmx-50

（李老師翻譯的甲烷水合自由能計算教程）

http://www.alchemistry.org/wiki/Absolute_Binding_Free_Energy_-_Gromacs_2016

（本教程使用的腳本來自於以上網站，對原作者的共享精神示以敬意）

J Comput Aided Mol Des (2015) 29: 397. https://doi.org/10.1007/s10822-015-9840-9

（本教程使用的分析腳本說明）

---

在此之前，還是簡要講解一些原理上的問題。

如何計算自由能？

要計算結合自由能，例如說一個主客體識別的結合自由能計算：碘——葫蘆脲系統，一種比較直觀的思路是所謂的牽拉方法——向已經與葫蘆脲絡合的碘施加一個諧振勢將碘一步步從葫蘆脲中慢慢拉扯出來（也就是傘型採樣動力學的方法，期間每一步都是要獨立模擬的，可以參閱https://jerkwin.github.io/GMX/GMXtut-3/），直到主客體已經離得足夠遠可被視為二者獨自在溶劑盒子中來處理的狀態，一系列傘形採樣模擬得到平均力勢，進而分析得到結合/解離過程的自由能變化。

使用這種方法有兩個缺點，一方面是你需要設定一個足夠大的模擬盒子（顯而易見）；另一方面，尤其是對於處理蛋白質-配體系統時，你的牽拉路徑未必是簡單的直線，而可能要穿過曲折的通道，這種情況下牽拉路徑是極難定義的。

鑒於以上，本教程採取一種更為建議的途徑——所謂「鍊金術（Alchemical）熱力學循環」來計算結合自由能。

這裡引入了一個拗口的名詞，「鍊金術熱力學循環」，這並不是我發明創造的名詞，老實說第一眼看見我也沒搞懂是什麼意思，不過弄懂後仔細想想覺得還是蠻恰當的一個比喻。讓我們回顧一下古老的鍊金術，在今天我們都知道鍊金術的天真想法——試圖從貧賤金屬中提煉出黃金——是不可行的，因為它試圖改變原子的種類，這在化學反應中是無法實現的。

但在計算化學中，改變一個原子是再輕鬆不過的事了，我們輕易地點幾下滑鼠就能替換掉一個分子中的任意原子或基團。「鍊金術熱力學循環」正是利用了這個特性——我們可以利用真實世界不存在的狀態來計算熱力學能，鑒於熱力學能是狀態函數，不依賴於路徑的選取，通過合理地構建熱力學態，我們可以找到一種方便計算的途徑，而不用考慮它本身是否真實。

舉個例子，例如計算甲烷水合自由能，利用鍊金術途徑，我們不需要將甲烷分子一步步從水裡拉到真空中，我們只需要慢慢地關閉（或打開）甲烷與水的相互作用（靜電力、L-J勢），也可以得到自由能。

拿蛋白質與配體的結合自由能來分析，一個講得比較清楚的圖示如下

在這個熱力學循環中，我們需要模擬的系統由它們周圍的黑色圓圈表示，限制勢的存在由回形針指示，白色配體意味著它不與環境相互作用（但仍存在分子內作用），藍色表示顯示水溶劑環境。注意在這裡引入了限制勢，限制勢的引入是為了防止配體在蛋白質中亂跑（尤其是關閉了大部分配體與環境的相互作用時）、阻礙收斂，只要最後在能量加上校正項就行。

這個計算將包含兩個模擬：配體在水中耦合（或解耦）/配體在蛋白中耦合（或解耦）。看上去我們似乎還需要計算C到D的自由能變化，但實際上C與D的狀態是完全相同的，反正此時的配體都不與周圍環境有任何相互作用了，換個環境對自由能改變為0。

需要注意的是在關閉相互作用時，我們總是先關閉靜電作用，再關閉非靜電項；如果反過來的話，那就會出現這樣的情況——正負電荷的原子可能會靠的無限近，這顯然是不合理的。反之，如果使用的是逐步打開相互作用的途徑，那就應該先打開非靜電項，再打開靜電項。

對每一步進行積分，即可得到總過程的自由能變化，這就是熱力學積分計算自由能的想法了。具體的公式以及實際操作將會在下一節講解。

推薦閱讀：

相关文章