本教程使用的分析腳本文獻：

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

完成了之前的模擬，我們得到了一系列dH/dλ文件，接下來就是進行分析來得到自由能了；分析的方法又有很多種，但總體而言可歸為兩類——自由能微擾（FEP）和熱力學積分法（TI），這一步的理論對於新手來說可能會顯得艱深。在這裡我強烈建議初學者優先學習熱力學積分法，這可以比較直觀地理解並利用dH/dλ文件的意義。

  • 首先從自由能的定義開始，Q為配分函數
  • displaystyle A = -eta^{-1} ln Q=-eta^{-1}ln int e^{-eta U(lambda,vec{q})}dvec{q}
  • 對λ求導，得到
  • displaystyle dA/dlambda = -eta^{-1}frac{d}{dlambda} ln int e^{-eta U(lambda,vec{q})}dvec{q} = -eta^{-1} frac{frac{d}{dlambda}int e^{-eta U(lambda,vec{q})}dvec{q}}{Q}
  • 改寫為
  • displaystyle dA/dlambda = -eta^{-1}frac{-etaint frac{dU(lambda,vec{q})}{dlambda} e^{-eta U(lambda,vec{q})}dvec{q}}{Q} = leftlangle frac{dU(lambda,vec{q})}{dlambda} ight angle_lambda
  • 最後，可以在整個λ範圍內進行積分，以得到最終的TI方程
  • displaystyle Delta A = int_0^1 leftlangle frac{dU(lambda,vec{q})}{dlambda} ight angle_lambda dlambda

當體積功忽略不計時（計算結合自由能的大部分情況），dU=dH，所以dA=dG。

從以上推導可以看出，實際上的自由能就是對dH/dλ進行積分。而實際上由於選取的λ是離散值，最終的積分也是近似的，亦即：

displaystyle Delta G approx sum_{k=1}^K w_k leftlangle frac{dH(lambda,vec{q})}{dlambda} ight angle_k

以上所說的看上去很複雜對不對？不過實際上你完全不懂這些也不要緊，因為已經有人開發了一鍵分析的python腳本alchemical-analysis，只需要配置好腳本後在dHdl文件夾內運行

alchemical-analysis -t 溫度 -s 平衡時間 等參數

即可，自動分析並調用一系列方法（可選）得到誤差和最終自由能

MobleyLab/alchemical-analysis?

github.com圖標

給沒有梯子的讀者下載地址：

鏈接: pan.baidu.com/s/1IfzoDR 提取碼: 421c

如果按照默認的分析參數，以及在mdp中設置了calc-lambda-neighbors = -1，那麼會得到以下一系列方法的分析結果

# Command line was: /bin/alchemical_analysis -t 300 -s 0 -u kcal -w -g

------------ --------------------- --------------------- --------------------- --------------------- --------------------- ---------------------
States TI (kcal/mol) TI-CUBIC (kcal/mol) DEXP (kcal/mol) IEXP (kcal/mol) BAR (kcal/mol) MBAR (kcal/mol)
------------ --------------------- --------------------- --------------------- --------------------- --------------------- ---------------------
0 -- 1 2.006 +- 0.012 2.009 +- 0.024 1.922 +- 0.024 1.915 +- 0.017 1.923 +- 0.011 1.923 +- 0.011
1 -- 2 4.961 +- 0.088 5.493 +- 0.101 6.041 +- 0.078 5.963 +- 0.418 5.691 +- 0.058 5.690 +- 0.058
2 -- 3 3.350 +- 0.094 3.504 +- 0.106 3.144 +- 0.513 1.967 +- 0.192 3.043 +- 0.169 3.043 +- 0.170
3 -- 4 0.033 +- 0.044 -0.263 +- 0.053 -0.251 +- 0.206 -0.389 +- 0.150 -0.108 +- 0.051 -0.108 +- 0.051
4 -- 5 -1.261 +- 0.064 -1.243 +- 0.073 -1.046 +- 0.077 -1.560 +- 0.191 -1.227 +- 0.062 -1.252 +- 0.061
5 -- 6 -1.952 +- 0.059 -2.030 +- 0.068 -1.947 +- 0.136 -1.946 +- 0.097 -1.960 +- 0.049 -1.925 +- 0.045
6 -- 7 -2.120 +- 0.065 -2.086 +- 0.076 -2.298 +- 0.043 -1.372 +- 0.349 -2.244 +- 0.033 -2.287 +- 0.032
7 -- 8 -2.406 +- 0.067 -2.403 +- 0.075 -2.436 +- 0.126 -2.499 +- 0.171 -2.515 +- 0.055 -2.608 +- 0.038
8 -- 9 -1.382 +- 0.013 -1.382 +- 0.014 -1.371 +- 0.021 -1.419 +- 0.019 -1.388 +- 0.013 -1.423 +- 0.010
9 -- 10 -1.469 +- 0.008 -1.476 +- 0.009 -1.426 +- 0.014 -1.500 +- 0.009 -1.493 +- 0.006 -1.469 +- 0.005
10 -- 11 -1.507 +- 0.005 -1.507 +- 0.008 -1.532 +- 0.005 -1.471 +- 0.016 -1.520 +- 0.004 -1.492 +- 0.003
11 -- 12 -0.602 +- 0.003 -0.602 +- 0.003 -0.600 +- 0.004 -0.603 +- 0.004 -0.602 +- 0.003 -0.601 +- 0.001
12 -- 13 -0.905 +- 0.003 -0.908 +- 0.004 -0.913 +- 0.005 -0.898 +- 0.002 -0.900 +- 0.002 -0.906 +- 0.001
13 -- 14 -1.370 +- 0.025 -1.440 +- 0.032 -1.157 +- 0.006 -1.340 +- 0.122 -1.194 +- 0.008 -1.116 +- 0.004
14 -- 15 -2.153 +- 0.025 -2.194 +- 0.028 -2.635 +- 0.096 -1.577 +- 0.069 -1.862 +- 0.014 -1.843 +- 0.005
15 -- 16 -5.884 +- 0.075 -5.713 +- 0.102 -7.271 +- 0.395 -4.572 +- 0.281 -6.370 +- 0.036 -5.308 +- 0.012
16 -- 17 -3.943 +- 0.038 -3.946 +- 0.038 -4.185 +- 0.084 -3.171 +- 0.050 -3.335 +- 0.018 -3.723 +- 0.006
17 -- 18 -4.701 +- 0.025 -4.686 +- 0.026 -5.253 +- 0.021 -4.168 +- 0.169 -5.148 +- 0.011 -5.101 +- 0.009
18 -- 19 -12.421 +- 0.075 -12.289 +- 0.091 -12.869 +- 0.330 -10.787 +- 0.467 -12.284 +- 0.072 -13.453 +- 0.043
19 -- 20 -17.461 +- 0.077 -17.403 +- 0.088 -18.187 +- 0.241 -15.696 +- 0.526 -17.415 +- 0.081 -17.454 +- 0.081
------------ --------------------- --------------------- --------------------- --------------------- --------------------- ---------------------
Coulomb: -47.933 +- 0.190 -47.670 +- 0.213 -51.556 +- 0.583 -41.310 +- 0.790 -47.609 +- 0.117 -47.997 +- 0.105
vdWaals: -3.254 +- 0.265 -2.893 +- 0.264 -2.713 +- 0.595 -3.813 +- 0.658 -3.298 +- 0.213 -3.415 +- 0.234
TOTAL: -51.187 +- 0.326 -50.563 +- 0.339 -54.269 +- 0.833 -45.123 +- 1.028 -50.907 +- 0.243 -51.413 +- 0.256

對於一個穩健的模擬，各種分析方法得到的結果一般不會相差太大

另外，你還會得到一個OMBAR鄰接矩陣

鄰接矩陣反映出某λ值下的模擬與相鄰λ值下的模擬的相空間重疊程度，對於一個穩健的模擬，其鄰接矩陣的主對角線和相鄰兩條對角線上的值都不應該為0。

另外一個文件便是以熱力學積分形式顯示的積分圖：

使用離散值近似積分運算，如上圖所示；其中范德華力部分和庫侖力部分分別積分。

原則上說λ越多積分越準確，但耗時也將大大增加；λ太少，或間距太大，相空間不能很好重疊，會導致積分結果與真實值誤差大。因此選取好適宜的λ間距值很重要。眼尖的讀者這時候可能看到了以上圖片中λ間距值並不等距，這倒並不是說就應該這麼取，實際的λ間距值應該考慮你的體系。

簡單而言，如果某個λ值附近的曲率很大，那麼離散近似帶來的誤差就會很大，這時候應該在這附近多取幾個間距小的λ值以降低誤差；反之，如果曲率較小，那麼可以少取幾個λ值以減少模擬用時。一般而言這種大麴率的情況在范德華力積分中容易出現，在庫侖力中較少，因此庫侖力對應的總λ值數目也較少。你可能需要進行多次模擬才能找到合適的λ值的分布。

另外，一般而言，為了得到可信的數據，一個模擬需要運行3-5個副本，得到的數據才具有較好的統計學可靠性。

至此，我們終於能夠完成自由能的計算了。但是得到的這個數據仍然會存在許多誤差，一方面是由於分子力場自身缺陷導致的系統性誤差，是無法解決的；而另一方面，我們需要做一些應該考慮進去的矯正，這是下一節會談到的內容。


推薦閱讀：
相关文章