Gromacs结合自由能计算教程（一）

以下文章来源于网络，如有侵权,请联系管理员删除!

分子模拟的一大用途就是对体系的平衡态进行模拟采样得到体系的某些热力学数据，其中的自由能计算是分子模拟被应用的较多的一个场景之一；计算自由能的基本思想是研究系统在无限长时间内的能量均值等价于其系综平均值，因而可以通过有限时间的取样来估计系综平均值，进而通过公式求解得到相关自由能数据。

虽然量子化学软件（如Gaussian）也能够处理自由能计算，其原理是把分子的振动，转动，平动等对熵有贡献的项，按照二次项展开，计算它的S与温度的关系，从而可以求出一定温度下的自由能变化。然而，对于动辄上万原子的大体系如蛋白-配体体系而言，哪怕是半经验的方法都极为耗时，更不要提及显式溶剂的计算（一些与溶剂有强烈作用的情况下隐式溶剂模型就不尽人意了）。

自由能计算本质上是处理考察体系从一个热力学态到另一个热力学态发生的能量变化，这可以指其本身的变化（化合价变化）或者是周围环境的变化（从真空到溶剂，或者是从溶剂到配位环境的变化）

因此自由能的计算任务可以解决如以下的问题：

1. 计算化合物在不同环境的分配系数（如计算立普妥在水和辛醇间的分配）；

2. 计算主客体的键合强度（如计算葫芦[8]脲与碘单质键合常数） ——由此衍生出的计算有改变蛋白质或底物的结构以检测它们键合的变化。

然而，一直以来分子模拟计算结合自由能所涉及到的统计热力学知识与一系列复杂的模拟参数设置与繁琐的数据分析，这一道道天然的门槛阻碍了这种方法在非MD专业研究者中的推广。本教程旨在帮助缺乏MD经验的人也能够通过现有模板稍微修改参数，很快上手炼金术自由能计算。

如何计算自由能？

要计算结合自由能，例如说一个主客体识别的结合自由能计算：碘——葫芦脲系统，一种比较直观的思路是所谓的牵拉方法——向已经与葫芦脲络合的碘施加一个谐振势将碘一步步从葫芦脲中慢慢拉扯出来（也就是伞型采样动力学的方法，期间每一步都是要独立模拟的，直到主客体已经离得足够远可被视为二者独自在溶剂盒子中来处理的状态，一系列伞形采样模拟得到平均力势，进而分析得到结合/解离过程的自由能变化。

使用这种方法有两个缺点，一方面是你需要设定一个足够大的模拟盒子（显而易见）；另一方面，尤其是对于处理蛋白质-配体系统时，你的牵拉路径未必是简单的直线，而可能要穿过曲折的通道，这种情况下牵拉路径是极难定义的。

鉴于以上，本教程采取一种更为建议的途径——所谓“炼金术（Alchemical）热力学循环”来计算结合自由能。

这里引入了一个拗口的名词，“炼金术热力学循环”，这并不是我发明创造的名词，老实说第一眼看见我也没搞懂是什么意思，不过弄懂后仔细想想觉得还是蛮恰当的一个比喻。让我们回顾一下古老的炼金术，在今天我们都知道炼金术的天真想法——试图从贫贱金属中提炼出黄金——是不可行的，因为它试图改变原子的种类，这在化学反应中是无法实现的。

但在计算化学中，改变一个原子是再轻松不过的事了，我们轻易地点几下鼠标就能替换掉一个分子中的任意原子或基团。“炼金术热力学循环”正是利用了这个特性——我们可以利用真实世界不存在的状态来计算热力学能，鉴于热力学能是状态函数，不依赖于路径的选取，通过合理地构建热力学态，我们可以找到一种方便计算的途径，而不用考虑它本身是否真实。

举个例子，例如计算甲烷水合自由能，利用炼金术途径，我们不需要将甲烷分子一步步从水里拉到真空中，我们只需要慢慢地关闭（或打开）甲烷与水的相互作用（静电力、L-J势），也可以得到自由能。

拿蛋白质与配体的结合自由能来分析，一个讲得比较清楚的图示如下

在这个热力学循环中，我们需要模拟的系统由它们周围的黑色圆圈表示，限制势的存在由回形针指示，白色配体意味着它不与环境相互作用（但仍存在分子内作用），蓝色表示显示水溶剂环境。注意在这里引入了限制势，限制势的引入是为了防止配体在蛋白质中乱跑（尤其是关闭了大部分配体与环境的相互作用时）、阻碍收敛，只要最后在能量加上校正项就行。

这个计算将包含两个模拟：配体在水中耦合（或解耦）/配体在蛋白中耦合（或解耦）。看上去我们似乎还需要计算C到D的自由能变化，但实际上C与D的状态是完全相同的，反正此时的配体都不与周围环境有任何相互作用了，换个环境对自由能改变为0。

需要注意的是在关闭相互作用时，我们总是先关闭静电作用，再关闭非静电项；如果反过来的话，那就会出现这样的情况——正负电荷的原子可能会靠的无限近，这显然是不合理的。反之，如果使用的是逐步打开相互作用的途径，那就应该先打开非静电项，再打开静电项。

对每一步进行积分，即可得到总过程的自由能变化，这就是热力学积分计算自由能的想法了。具体的公式以及实际操作将会在下一节讲解，敬请期待。