李鵬飛 王美婷 梅曄

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

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

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

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

0引言

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

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

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

1方法

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

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

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

1.3.1平衡態(tài)動(dòng)力學(xué)模擬細(xì)節(jié)

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

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

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

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

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

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

3結(jié)論

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

（2）從計(jì)算效率來(lái)看，與平衡方法相比，非平衡方法的計(jì)算效率更高，其計(jì)算所需cPu時(shí)間和實(shí)際時(shí)間更少.