平衡态与非平衡态分子溶剂化自由能的计算效率比较

李鹏飞 王美婷 梅晔

摘要：着眼于13个中性氨基酸侧链类似物在水中的溶剂化自由能的计算，来比较两种计算自由能的平衡态动力学模拟和非平衡态动力学模拟方法在高性能计算机上的表现.研究发现，利用非平衡态动力学模拟来计算自由能除了在准确度上和平衡态动力学模拟的计算一致之外，在计算效率和实际所需时间上，非平衡方法计算效率更高，实际所需时间更少.

关键词：溶剂化自由能;氨基酸侧链类似物;平衡态动力学模拟;

非平衡态动力学模拟

中图分类号：0469 文献标志码：A DOI：10.3969/j.issn.1000-5641.2019.01.010

0引言

自由能是一个决定很多化学及生物反应过程的重要物理量.自由能的准确计算是计算生物物理领域中的一个重要方向Hansen和van Gunsteren曾提出在自由能计算时需要解决主要的3个问题.①相空间的采样是否充足;②描述体系的哈密顿量是否准确;③自由能计算方法是否可靠.我们知道，氨基酸侧链类似物的溶剂化自由能及其在不同溶剂之间的配分系数已经被广泛地计算.由于氨基酸侧链类似物的分子都是小体系，而且目前高性能计算机已经广泛普及，所以在对体系的相空间采样上是没有问题的.但是不同的自由能计算模拟方法对于计算机模拟计算的并行要求有所区别.假设在计算机资源充足的前提下，哪种方法并行效率更高实际所需时间更短，这种方法就更为实用.本文中，参数化已经很成熟的分子力学fMolecular Mechanics，MM）哈密顿量被用来描述所研究体系.本文的重点在于比较平衡态动力学模拟方法以及非平衡态动力学模拟方法在高性能计算机上计算氨基酸侧链类似物的溶剂化自由能的准确度、精度以及时间效率.

平衡态自由能计算方法主要包括热力学微扰法（Thermodynamic Perturbation，TP）;热力学积分法（Thermodynamic Integration，TI）[12];Bennett接受比例法（Bennetts Accep-tance Ratio，BAR）以及多态Bennett接受比例法（Multiple Bennett's Acceptance Ratio，MBAR）.这几种方法在计算自由能时具体表现的比较不作为本文的工作重点，此方面可以查阅参考文献.与TI和TP方法相比较，BAR方法具有最小化误差的优点.存在多个模拟窗口时，MBAR方法具有全局最小化誤差的优点.可以严格证明的是，BAR方法是MBAR方法的一种特例.换句话说，只存在两个模拟窗口时，MBAR方法就是BAR方法.本文中，计算分子溶剂化自由能时，为了保证相邻两态之间相空间的充分重叠，平衡态模拟引入了很多中间态，因此本文选取MBAR作为平衡态自由能计算方法.

非平衡态自由能计算方法的思想主要由Jarzynski于1997年首次提出.通过对两态之间的不可逆转变时所消耗功的指数平均建立与两态之间自由能差的等式关系.1998年，Crooks从微观上可逆的马尔科夫系统出发，重新推导了Jarzynski等式.次年，Crooks又发现两态之间非平衡转化的两个不同方向过程功的分布存在一定关系，被称之为“Crooks”关系.实际上，在多条非平衡轨迹模拟之后，再利用Jarzynski等式计算自由能差值时，由于有限的轨迹模拟有时会遇到收敛慢的问题.Shirts于2003年从“Crooks”关系出发重新推导了BAR的公式，可以有效解决这个问题.非平衡态动力学模拟及其在自由能计算方法中工作有很多，这里也不作为本重点，此方面亦可查阅相关参考文献.本文中，Shirts从非平衡角度重新推导BAR的方法被用作非平衡态自由能计算方法.

1方法

自由能是一个状态函数.两态之间的自由能差取决于两者状态，与转化路径无关.本文所使用的平衡动力学和非平衡动力学计算自由能的示意图，详见图1.

1.1平衡态动力学自由能计算方法

如图1（a）所示，在平衡态动力学模拟中，左右两端分别是溶质分子在水中和在气相下的状态.中间增加了一系列的中间态，它们的体系哈密顿量可以表示为H（λ）=（1-λ）H0.0+λH1.0，其中，H0.0和H1.0分别代表着左右两端的状态.中间态的引入是为了增加相邻两个状态窗口之间的相空间重叠程度，以保证自由能计算的可靠性.本文中，这些态的自由能由MBAR[14]方法计算得到.在MBAR方法中，对于K个热力学状态，每个状态动力学采样产生了Ni个无相关的平衡结构，第i个热力学态的自由能A可以写成所有状态自由能的函数，

1.3.1平衡态动力学模拟细节

如图1（a）所示，在平衡态动力学模拟中，每个分子的溶剂化自由能的计算过程分成了11个平衡态动力学模拟窗口，并且利用到了“软核势”，可以平滑地将小分子与水溶剂之间的相互作用关掉.在每一个模拟窗口中，先将体系缓慢加热到298K，然后进行200 ps的预平衡，最后进行了10ns的NPT平衡系综模拟.在这一系列的过程中，非键相互作用的实空间截断在12 A，长程库仑相互作用计算采用PME方法.

1.3.2非平衡态动力学模拟细节

如图1（b）所示，在非平衡态动力学模拟中，每个分子的溶剂化自由能可以由Ⅳ条非平衡态动力学模拟轨迹中功值的指数平均得到.本文对于每个分子的溶剂化自由能计算，分别进行了正向反向各100条非平衡轨迹的模拟.在每条轨迹中，相邻λ以0.2 ps为间隔均匀变化100次，使得λ从0变化到1.因此，每条轨迹的动力学模拟时间是20 ps.在非平衡动力学模拟中的其他条件保持与平衡动力学模拟相同.需要说明的一点是，非平衡动力学模拟的初始结构来自于体系处于初态哈密顿量下的系综.在初态平衡系综模拟中，先将体系进行优化，然后缓慢将体系升温到298K，然后进行200ps的预平衡，最后进行2ns的平衡系综采样，该采样过程中的结构被用作非平衡动力学模拟的初始结构.

2.2平衡态动力学模拟和非平衡态动力学模拟的计算效率比较

假设在计算资源充足的情况下，比较两种方法的计算效率.对于平衡态动力学模拟，由M个模拟窗口来完成，可以同时由M个计算节点来完成.在本文中，每个窗口的动力学模拟时间是10 ns.以Asn体系为例，所使用计算机是华东师范大学第4期IBM超算，使用单节点16核数模拟一个10 ns的平衡轨迹的实际时间大约是5.2 h，平均1 ns需要0.5 h.从计算结果分析来看，对于该体系每个窗口6 ns的模拟结果是已经收敛了的，这样每个窗口只需要3 h.

而对于非平衡态动力学模拟，每条轨迹的计算彼此之间是不需要进行通讯的，换句话说，Ⅳ条轨迹是完全可以独立进行模拟，同时由Ⅳ个计算节点来完成.本文的每条轨迹的动力学模拟时间是20 ps.以Asn体系为例，所使用计算机同样是华东师范大学第4期IBM超算，使用单节点16核数模拟该体系的一条20 ps轨迹的实际时间大约是79 s.当然对于这个体系而言，在非平衡态模拟之前是要在初态下进行平衡系综采样.在本文中，这个体系的初态平衡模拟时间是2ns，同样的计算节点完成该计算实际所需时间0.9 h.由此可以看出，非平衡方法的并行效率更高，计算所需CPU时间和实际时间都远远少于平衡方法.

3结论

（1）平衡态动力学模拟和非平衡态动力学模拟的计算结果完全一致.两种方法都能给出比较准确的计算结果.

（2）从计算效率来看，与平衡方法相比，非平衡方法的计算效率更高，其计算所需cPu时间和实际时间更少.