平衡態(tài)與非平衡態(tài)分子溶劑化自由能的計算效率比較

李鵬飛 王美婷 梅曄

摘要：著眼于13個中性氨基酸側(cè)鏈類似物在水中的溶劑化自由能的計算，來比較兩種計算自由能的平衡態(tài)動力學(xué)模擬和非平衡態(tài)動力學(xué)模擬方法在高性能計算機上的表現(xiàn).研究發(fā)現(xiàn)，利用非平衡態(tài)動力學(xué)模擬來計算自由能除了在準(zhǔn)確度上和平衡態(tài)動力學(xué)模擬的計算一致之外，在計算效率和實際所需時間上，非平衡方法計算效率更高，實際所需時間更少.

關(guān)鍵詞：溶劑化自由能;氨基酸側(cè)鏈類似物;平衡態(tài)動力學(xué)模擬;

非平衡態(tài)動力學(xué)模擬

中圖分類號：0469 文獻標(biāo)志碼：A DOI：10.3969/j.issn.1000-5641.2019.01.010

0引言

自由能是一個決定很多化學(xué)及生物反應(yīng)過程的重要物理量.自由能的準(zhǔn)確計算是計算生物物理領(lǐng)域中的一個重要方向Hansen和van Gunsteren曾提出在自由能計算時需要解決主要的3個問題.①相空間的采樣是否充足;②描述體系的哈密頓量是否準(zhǔn)確;③自由能計算方法是否可靠.我們知道，氨基酸側(cè)鏈類似物的溶劑化自由能及其在不同溶劑之間的配分系數(shù)已經(jīng)被廣泛地計算.由于氨基酸側(cè)鏈類似物的分子都是小體系，而且目前高性能計算機已經(jīng)廣泛普及，所以在對體系的相空間采樣上是沒有問題的.但是不同的自由能計算模擬方法對于計算機模擬計算的并行要求有所區(qū)別.假設(shè)在計算機資源充足的前提下，哪種方法并行效率更高實際所需時間更短，這種方法就更為實用.本文中，參數(shù)化已經(jīng)很成熟的分子力學(xué)fMolecular Mechanics，MM）哈密頓量被用來描述所研究體系.本文的重點在于比較平衡態(tài)動力學(xué)模擬方法以及非平衡態(tài)動力學(xué)模擬方法在高性能計算機上計算氨基酸側(cè)鏈類似物的溶劑化自由能的準(zhǔn)確度、精度以及時間效率.

平衡態(tài)自由能計算方法主要包括熱力學(xué)微擾法（Thermodynamic Perturbation，TP）;熱力學(xué)積分法（Thermodynamic Integration，TI）[12];Bennett接受比例法（Bennetts Accep-tance Ratio，BAR）以及多態(tài)Bennett接受比例法（Multiple Bennett's Acceptance Ratio，MBAR）.這幾種方法在計算自由能時具體表現(xiàn)的比較不作為本文的工作重點，此方面可以查閱參考文獻.與TI和TP方法相比較，BAR方法具有最小化誤差的優(yōu)點.存在多個模擬窗口時，MBAR方法具有全局最小化誤差的優(yōu)點.可以嚴格證明的是，BAR方法是MBAR方法的一種特例.換句話說，只存在兩個模擬窗口時，MBAR方法就是BAR方法.本文中，計算分子溶劑化自由能時，為了保證相鄰兩態(tài)之間相空間的充分重疊，平衡態(tài)模擬引入了很多中間態(tài)，因此本文選取MBAR作為平衡態(tài)自由能計算方法.

非平衡態(tài)自由能計算方法的思想主要由Jarzynski于1997年首次提出.通過對兩態(tài)之間的不可逆轉(zhuǎn)變時所消耗功的指數(shù)平均建立與兩態(tài)之間自由能差的等式關(guān)系.1998年，Crooks從微觀上可逆的馬爾科夫系統(tǒng)出發(fā)，重新推導(dǎo)了Jarzynski等式.次年，Crooks又發(fā)現(xiàn)兩態(tài)之間非平衡轉(zhuǎn)化的兩個不同方向過程功的分布存在一定關(guān)系，被稱之為“Crooks”關(guān)系.實際上，在多條非平衡軌跡模擬之后，再利用Jarzynski等式計算自由能差值時，由于有限的軌跡模擬有時會遇到收斂慢的問題.Shirts于2003年從“Crooks”關(guān)系出發(fā)重新推導(dǎo)了BAR的公式，可以有效解決這個問題.非平衡態(tài)動力學(xué)模擬及其在自由能計算方法中工作有很多，這里也不作為本重點，此方面亦可查閱相關(guān)參考文獻.本文中，Shirts從非平衡角度重新推導(dǎo)BAR的方法被用作非平衡態(tài)自由能計算方法.

1方法

自由能是一個狀態(tài)函數(shù).兩態(tài)之間的自由能差取決于兩者狀態(tài)，與轉(zhuǎn)化路徑無關(guān).本文所使用的平衡動力學(xué)和非平衡動力學(xué)計算自由能的示意圖，詳見圖1.

1.1平衡態(tài)動力學(xué)自由能計算方法

如圖1（a）所示，在平衡態(tài)動力學(xué)模擬中，左右兩端分別是溶質(zhì)分子在水中和在氣相下的狀態(tài).中間增加了一系列的中間態(tài)，它們的體系哈密頓量可以表示為H（λ）=（1-λ）H0.0+λH1.0，其中，H0.0和H1.0分別代表著左右兩端的狀態(tài).中間態(tài)的引入是為了增加相鄰兩個狀態(tài)窗口之間的相空間重疊程度，以保證自由能計算的可靠性.本文中，這些態(tài)的自由能由MBAR[14]方法計算得到.在MBAR方法中，對于K個熱力學(xué)狀態(tài)，每個狀態(tài)動力學(xué)采樣產(chǎn)生了Ni個無相關(guān)的平衡結(jié)構(gòu)，第i個熱力學(xué)態(tài)的自由能A可以寫成所有狀態(tài)自由能的函數(shù)，

1.3.1平衡態(tài)動力學(xué)模擬細節(jié)

如圖1（a）所示，在平衡態(tài)動力學(xué)模擬中，每個分子的溶劑化自由能的計算過程分成了11個平衡態(tài)動力學(xué)模擬窗口，并且利用到了“軟核勢”，可以平滑地將小分子與水溶劑之間的相互作用關(guān)掉.在每一個模擬窗口中，先將體系緩慢加熱到298K，然后進行200 ps的預(yù)平衡，最后進行了10ns的NPT平衡系綜模擬.在這一系列的過程中，非鍵相互作用的實空間截斷在12 A，長程庫侖相互作用計算采用PME方法.

1.3.2非平衡態(tài)動力學(xué)模擬細節(jié)

如圖1（b）所示，在非平衡態(tài)動力學(xué)模擬中，每個分子的溶劑化自由能可以由Ⅳ條非平衡態(tài)動力學(xué)模擬軌跡中功值的指數(shù)平均得到.本文對于每個分子的溶劑化自由能計算，分別進行了正向反向各100條非平衡軌跡的模擬.在每條軌跡中，相鄰λ以0.2 ps為間隔均勻變化100次，使得λ從0變化到1.因此，每條軌跡的動力學(xué)模擬時間是20 ps.在非平衡動力學(xué)模擬中的其他條件保持與平衡動力學(xué)模擬相同.需要說明的一點是，非平衡動力學(xué)模擬的初始結(jié)構(gòu)來自于體系處于初態(tài)哈密頓量下的系綜.在初態(tài)平衡系綜模擬中，先將體系進行優(yōu)化，然后緩慢將體系升溫到298K，然后進行200ps的預(yù)平衡，最后進行2ns的平衡系綜采樣，該采樣過程中的結(jié)構(gòu)被用作非平衡動力學(xué)模擬的初始結(jié)構(gòu).

2.2平衡態(tài)動力學(xué)模擬和非平衡態(tài)動力學(xué)模擬的計算效率比較

假設(shè)在計算資源充足的情況下，比較兩種方法的計算效率.對于平衡態(tài)動力學(xué)模擬，由M個模擬窗口來完成，可以同時由M個計算節(jié)點來完成.在本文中，每個窗口的動力學(xué)模擬時間是10 ns.以Asn體系為例，所使用計算機是華東師范大學(xué)第4期IBM超算，使用單節(jié)點16核數(shù)模擬一個10 ns的平衡軌跡的實際時間大約是5.2 h，平均1 ns需要0.5 h.從計算結(jié)果分析來看，對于該體系每個窗口6 ns的模擬結(jié)果是已經(jīng)收斂了的，這樣每個窗口只需要3 h.

而對于非平衡態(tài)動力學(xué)模擬，每條軌跡的計算彼此之間是不需要進行通訊的，換句話說，Ⅳ條軌跡是完全可以獨立進行模擬，同時由Ⅳ個計算節(jié)點來完成.本文的每條軌跡的動力學(xué)模擬時間是20 ps.以Asn體系為例，所使用計算機同樣是華東師范大學(xué)第4期IBM超算，使用單節(jié)點16核數(shù)模擬該體系的一條20 ps軌跡的實際時間大約是79 s.當(dāng)然對于這個體系而言，在非平衡態(tài)模擬之前是要在初態(tài)下進行平衡系綜采樣.在本文中，這個體系的初態(tài)平衡模擬時間是2ns，同樣的計算節(jié)點完成該計算實際所需時間0.9 h.由此可以看出，非平衡方法的并行效率更高，計算所需CPU時間和實際時間都遠遠少于平衡方法.

3結(jié)論

（1）平衡態(tài)動力學(xué)模擬和非平衡態(tài)動力學(xué)模擬的計算結(jié)果完全一致.兩種方法都能給出比較準(zhǔn)確的計算結(jié)果.

（2）從計算效率來看，與平衡方法相比，非平衡方法的計算效率更高，其計算所需cPu時間和實際時間更少.