Gromacs結合自由能計算教程（二）
如何計算自由能？
續上回講到的，上一章簡要介紹了自由能計算的相關背景與鍊金術熱力學途徑。這一節將會講解計算鍊金術自由能中每一步熱力學過程所涉及的自由能計算方法。
首先回顧一下鍊金術熱力學的途徑，我們會涉及到一個講研究對象從其環境中解耦和的過程，在這個過程中我們會先後且逐步地關閉與環境的靜電作用與范德華作用，保證相鄰模擬間儘可能大的相空間重疊，降低自由能計算的不準確性。再次強調這個過程是分開進行的，並且靜電作用總是先於范德華力解耦以避免正負電荷過於近的不合理接觸（范德華力的交換互斥項可以避免這種不合理接觸）。
另外，正如你此時可能會感到疑問的一點是，我在這裡所述的范德華力與我們大學課本中所講述的似乎略有出入。
在大學課本中的范德華力是基本上是這麼說的：
「范德華力的主要來源有三種機制：
取向力：極性分子之間的相互作用力誘導力：極性分子極化非極性分子，產生誘導偶極矩，並相互作用
色散力：一對非極性分子由於本身電子的概率運動，可以相互配合產生一對方向相反的瞬時偶極矩，這一對瞬時偶極矩相互作用。這種機制是非極性分子中范德華力的主要來源，1930年由F.W.倫敦首先根據量子力學原理給出解釋「 ——from Wikipedia
到目前為止我所說的范德華力實際指的大約（說是大約是因為這樣說也不準確）只有第三項——色散力，嚴謹的表達是分子動力學軟體用來描述原子間作用的所謂L-J勢，包含了近距離的交換互斥項與遠距離的吸引項。搞清楚分子動力學中的范德華力這個基礎概念十分重要，我會在之後的文章中單獨論述。事實上許多中文教材在范德華力的描述相當old school，存在著自相矛盾與糾纏不清的地方。關於這一點，強烈推薦閱讀維基上范德華力的英文詞條，講得清晰明了得多：
Van der Waals force
回到鍊金術自由能計算。我們在這裡引入了一個變數λ來線性地描述解耦過程，Gromacs提供了一系列關鍵字如「couple-lambda」等來實現相關操作（前提是使用關鍵字「free-energy = yes」）。在熱力學積分法中，我們計算自由能的方法實際上是使用?U/?λ對λ進行積分，λ實際取的是離散值，所以我們得到的自由能實際也是近似值；原則上如果λ無限細分，所得的自由能就可以無限逼近真實值。
如上圖就是先後解耦庫倫作用與范德華力整個過程的熱力學積分示意。
關於熱力學積分的更具體的物理原理可以參閱，裡面也介紹了熱力學微擾法（FEP）的推導，之後我們也會提到
自由能微擾和熱力學積分方法 - 分子模擬 - 計算化學公社事實上按照這個鍊金術途徑處理自由能的方法有很多種，在之後的分析結果會進行進一步的討論。
值得一提的是，為了保證?A/?λ函數的連續性，在關閉范德華力時我們總是會使用「軟核勢能」來描述原子，著可以實現真正意義上的「逐步關閉原子相互作用」，關於軟核勢能的理論文章可以參考
10.1002/jcc.21909前戲這麼長，估計大家都已經開始失去耐心了。現在讓我們開始上手自由能計算在Gromacs的操作流程吧。
Gromacs中關於自由能的參數設置
自由能計算在Gromacs中主要涉及以下額外關鍵詞
free-energy = yes
couple-moltype = ligand
couple-lambda0 = none
couple-lambda1 = vdw-q
couple-intramol = no
separate-dhdl-file = yes
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
init-lambda-state = 0
coul-lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.00 0.0 0.00 0.00 0.0 0.1 0.2 0.4 0.5 0.6 0.8 1.0
vdw-lambdas = 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.85 0.9 0.95 0.97 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
nstdhdl = 100
calc-lambda-neighbors = -1
Couple-moltype是你要耦合（或解耦）的分子（基團）名稱
Couple-lambda0是你期望的起點狀態，Couple-lambda1是終點狀態；在給出的例子中，Couple-lambda0=none，Couple-lambda1=vdw-q意思是在初始狀態不與環境相互作用，在最終狀態完全打開范德華力與靜電作用，所以著定義的是一個耦合的過程。如果是解耦過程應該反過來寫。
couple-intramol =no 指的是在這個過程中不處理分子內去耦合/耦合的問題
sc-alpha, sc-power, sc-sigma三項是關於軟核作用勢的設定，需要參考軟核作用勢的相關文獻。不知道軟核作用勢的情況下就用這套值也沒太大問題。
init-lambda-state 指的是當前模擬的λ值。事實上在一個計算中會涉及十幾個或更多λ下的計算。可以使用腳本mk_mdp.sh 一鍵生成不同lambda的mdp文件（將在下一節一併給出壓縮包），給不同過程調用。
vdw-lambdas與coul-lambda對應每個λ下的范德華力與靜電力的耦合量，0代表沒有相互作用，1代表完全打開。這些值的設定有一定技巧性，需要根據模擬體系進行調整，之後會進一步討論。在這裡是一個耦合過程，所以先打開了范德華作用，再打開靜電作用。
nstdhdl = 100指的是每100幀保存一次?H/?λ數據
calc-lambda-neighbors = -1指的是計算並輸出所有lambda值的ΔH供MBAR分析。通常情況下，在設置了init-lambda-state時，calc-lambda-neighbors = 1指的是計算當前lambda下前後一個lambda值對應的ΔH，對於普通BAR分析已經足夠了，而對於MBAR分析，需要使用-1，表示計算所有值。關於MBAR分析之後會提到
推薦閱讀：
