自由能计算，PMF，伞形采样，WHAM

一个过程的自由能的分子模拟计算被称为分子模拟领域的四个主要问题之一. 可能是因为自由能的概念更加令人困惑. 下面说说我的理解，每个人都在打砖块.

要谈论自由能，首先要讨论其“能量”与其兄弟之间的区别. 每个做过模拟的人都知道，一个模拟系统，无论是NPT还是NVT，都具有总能量，即内部能量. 内部能量分为动能和势能两种. 势能项可以分解为许多子项，即键能，非键相互作用能等，这一点很清楚. 无论模拟是否平衡，总能量是否守恒，以及Plot的总能量，均清晰可见. 但是什么是自由能？

自由能是在物理化学中教授的，但是教科书使用一堆数学公式严格包装了物理意义，使人们很难看到其物理性质. 我们必须首先阐明，自由能是自发过程的标准: 自发能量减少的过程可以自发进行. 一个可以自发进行的过程必须减少其自由能. 这是热力学的第二定律.

以两个分子的结合过程为例，结合能和结合自由能有什么区别？ “束缚能量”很容易理解. 真空中的两个分子相距很远，干燥的物料彼此接触. 系统的总能量减少. 释放的能量是“结合能”，它是分子之间的吸引力的表现（由分子之间的永久偶极或瞬态偶极引起）. 结合能是内部能量的变化量. 但这显然不是“与自由能结合”. 什么是“结合自由能”？那些已经学会了物化的人将立即想到自由能无序定义中的“熵”. 是的，在联合之前，这两个分子喜欢去哪里，拥有多少自由，现在他们在一起，必须在一起，并且无序程度明显降低了. 如果反应在溶液中进行，则不仅减少了结合中涉及的分子的无序性，而且还影响了周围的环境分子. 让我将这两个分子保留在一定距离内的水中. 假设这两个分子是苯，它们是疏水的，那么这两个人很快就会在水中聚在一起. 这个过程显然是一个自发的过程，能量在前后变化，系统的总能量仍然减少. 尽管这种能量包括环境能量的变化，但是这种能量变化与自由能结合吗？不，它仍然是内部能量的变化，减少的能量以热量的形式释放. 但是，在结合前后，这两种苯对周围水分子的排列有影响，导致熵变. 自由能变化=内部能变化+熵能变化，

ΔA＝ΔU-TΔS. 为什么这里有减号？在苯和苯混合之前分子数密度公式，先将黑白混合，并且熵必须很大. 现在黑色是黑色，白色是白色分子数密度公式，熵明显减小了. 上帝的老人显然不喜欢秩序，至少不喜欢这种对平衡的态度，因此熵降低了，但是相应的熵能量却必须增加，这是非常不利的. 因此，TΔS在公式中为负. 但是您会问，熵S是如何如此神秘的东西，如何测量和模拟它？不要先忙在这里，我们首先总结结合自由能ΔA和结合能ΔU之间的差异: ΔA=ΔU-TΔS. 也就是说，结合能ΔU是结合自由能ΔA的一部分，另一部分是熵能TΔS的变化.

以上说明适用于NVT系统，因为体积不变，所以没有抵抗外部压力的功PΔV. 如果是NPT系统，则必须更改盒子的体积. 此时，该自由能称为吉布斯自由能，另外还有一个外部功PΔV，即，ΔG＝ΔU+PΔV-TΔUS. 如果将前两个变化合并在一起，则称为焓变化ΔH=ΔU+PΔV. 这是题外话. 那么以前的NVT系统的自由能变化是什么？是的，亥姆霍兹自由能. 以下两个项目之间没有区别，它们统称为自由能. 具体的自由能取决于您使用的是NVT还是NPT.

在继续之前，高斯定量版本的兄弟问: 两个分子的自由能如何在前面的真空中结合？实际上，您的计算全部为0K. 热力学第三定律说，熵在0K时为0，并且熵能没有变化. 因此，您的组合自由能等于组合能. 这里应该补充一点，高斯还可以计算除0k以外的两个分子的结合自由能. 该如何计算？就是扩大分子的振动，旋转，平移等对熵有贡献的项，根据二次项，计算其S与温度之间的关系，以便在一定温度下的自由能变化为获得. 它很远，然后回到我们最关心的问题: 如何模拟过程中自由能的变化.

但是，在我开始之前，我必须举一个例子. 这个例子可以帮助我们对系统的散热/工作有一个清晰的了解. 这并不意味着总的自由能会发生变化: 想象一下，将盐撒到一杯水中，这时盐迅速溶解了. 一些盐溶解以释放热量，而另一些盐溶解以吸收热量. 吸热意味着系统的能量将增加. 请问，为什么溶解过程会自发吸收热量？

答案很明确: 可以自发进行的过程表明，系统的总自由能减少了. 自由能分为内能和熵能两部分. 尽管吸热能使内能升高，但是只要能够充分降低熵，最终的自由能仍会降低，因此可以自发进行. 这种熵降低的根源是什么？显然，原始阴离子和阳离子仅在一种状态下结合在一起. 溶解后，它们可以游动，活动空间大大增加，这意味着状态总数大大增加，熵也大大增加. . . 但是应该注意的是，尽管离子的熵增加了，但是水分子呢？由于部分水分子被离子束缚在其壳中，因此活性降低，因此某些水分子的熵可以略微增加. 但是，它不会影响整体情况.

现在，您应该已经完全了解自由能，内部能和熵能之间的区别. 下面切入主题，如何使用分子模拟来模拟过程的自由能变化？例如，底物分子和蛋白质分子的组合？还是沸石吸收了H2？还是一个小分子穿过双层膜？这些都是变化. 任何更改过程都可以使用一种称为反应坐标的东西将过程从头到尾连接起来. 例如，如果小分子通过某个通道，则该通道的中心轴可以定义为反应坐标. 该示例相对简单，并且很难为某些复杂过程定义反应坐标. 例如，您是否给出了在沸石表面吸附许多H2的反应坐标？您只能关注一个H2，但是很难确定要吸附的方位角.

如前所述，自由能中有一种叫做熵的东西. 可以模拟这个东西吗？不行熵的定义称为“系统所有可访问的微观状态的数量的对数”，S =klnΩ，Ω，无法计算. 此外，我们都知道分子模拟可以一步一步地计算出样品. 如果样本对应于一个状态，则意味着您可以通过执行一些步骤来拥有一些状态. . . 因此，分子模拟不能模拟熵. 严格来说，应该说“熵不是系统微观状态的平均热力学量”，或者“熵是与相空间的体积有关的量”. 只有那些具有相应微观热力学量的宏观特性才能被模拟和平均，例如体积（例如内部能量）（例如压力）.

好的，您可以模拟内部能量，但是不能模拟熵. 您如何计算自由能？ ？ ？

兄弟，在讨论这个问题之前，让我们先看看实验如何确定反应前后自由能的差异. 记住，您如何测量？ . . 测量浓度的变化. . 通过反应平衡后的浓度分布，可以获得平衡常数（产物浓度m /反应物浓度n，m，n是平衡系数，您知道），平衡常数为K，则∆G°= -RT ln K . . . 但是这种魔力是如何改变的？关键在于该平衡常数，该常数反映了系统中分子的选择: 是选择反应态还是产物态. . 如果需要产品状态，则K值较大；如果渴望反应物质状态，则K值小. . . 这一向往体现了对分子的直觉选择: “无论酷还是呆在哪里”，以及“任何自由都应尽可能低”的自由意志. . . 因此，最终浓度反映了由于自由能下降导致的两个状态之间的分布状态的差异，即概率密度！ ！如果您能理解，恭喜您，这很容易理解. 为了模拟自由能，分子模拟将从分布概率开始！ ！ ！

分子模拟如何模拟自由能的变化？简而言之，如果系统具有两个状态A和B，则让系统运行. . . 快跑. . 运行统计后，我在A状态下花费了20％的时间，在B状态下花费了P2 = 80％的时间. 概率为PA = 0.2，PB = 0.8，那么，平衡常数K = P2 / P1 = 4，AB两个状态的自由能差为∆G°= -RT lnK = -0.6ln4≅0.8 kcal / mol （T〜300K，kT约为0.6kcal / mol，这就是模拟自由能变化的分子模拟原理.

要添加2点.

首先，真实系统不能只有两个状态A和B. 例如，沿着反应坐标，当两个分子彼此接近时，存在无数个状态. 如果您离散化反应坐标并将其切成小块，则状态数取决于您如何切割它们. 您可以放开系统，让它在每个段上自由跳转. 最后，您算出每个段上的舞步数. 最终，您可以获得概率密度分布，最后使用公式将其转换为每个线段（状态），其上的自由能被连接成一条起伏的曲线，称为PMF，全称是平均力的势能. 它是传说中的自由能势能面. 是的，能量变化的曲线或曲面称为势能曲面；自由能势能面称为PMF. 为什么要取一个如此无休止的名字呢？当这两个分子接近时，整个系统中的所有其他分子都会对其产生影响（力），最终结果是平均效果，因此称为电位

平均力. 无论如何，它真的在Wiki上得到了解释. 将此限制为两个分子相互接近的过程. 但是我认为，只要您可以定义反应坐标，就可以将其背后的思想应用于任何更改过程. 然后，坐标上每个点的概率也反映了所有分子对该点的平均影响. 对吧？

第二，概率密度如何转换为PMF值？回忆一个经典的例子，它是模拟水中两个甲烷分子的组合自由能势能面: 释放这两个分子使其在水中运行，最后计算它们之间的径向分布函数g®，也称为RDF， RDF恰好是两个距离不同的分子的概率密度. 然后根据以下公式计算每个点的PMF值:

W®= -kTln

W®是一条优美的曲线，它沿反应坐标（即两者之间的距离r）结合了自由能. 通常，您会在一定距离处找到一个最小点，而该点是它们之间最快乐的距离. 上式是所谓的柯克伍德公式. 现在您应该对正在发生的事情有一个模糊的想法: 自由能是耗尽的可能性. 这个概率与它有什么关系？反应坐标势能面U®. 为什么？因为仿真系统根据某个概率访问一个点，所以该概率取决于该点的势能U. 势能越高，进入的可能性越低. 您可能立即想到了Boltzmann分布，对，对，

P = kexp（-U / kT）

请注意，此概率P是整个系统访问某个能量状态的概率. 我们要求的概率是某个反应坐标下的概率. 系统分为反应部分和环境部分，关注的是反应一些不同状态之间的相对概率. 也就是说，假设r代表反应坐标

P®= k exp（-U®/ kT）

通过该公式，我们现在已成功地将内部能势能面U®转换为概率P®，然后将该概率转换为自由能势能面W®，包括来回熵！ ！ ！太神奇了吗？关键在于玻尔兹曼分布. 分布函数是势能面和自由能势能面之间的链接. “补充说明: 关键是玻耳兹曼公式中的kT. 以热驱动的形式，kT将分子从其喜欢的低势能状态推到了它不喜欢的高势能状态，从而导致分子数增加相空间（相空间的体积增加）中可访问状态的数量；温度越高，kT驱动的状态数量就越多. 这就是熵的包含方式.

好的，您现在会说很简单，只要我们知道沿反应坐标的分布概率，就可以找到目标反应系统的自由能势能面. 那么如何找到响应坐标分布的概率呢？最直接的想法之一是将反应坐标分为多个点，然后在每个点进行模拟，然后计算每个点的概率密度Pi，最后将其转换为自由能势能面. 现在，关键是找到每个点的分布概率.

在前面的两个甲烷示例中，我们逃跑了. 现在，我想告诉您不幸的消息. 由于分子模拟本身的缺陷，这种“放手”方法不适用于大多数变化.

范玮: 发生了什么事？这种“开放运行”方法要求系统在反应坐标上的每个点都进行良好采样. 然后我问您，您的模拟会自动访问所有点吗？通常不会. 为什么？如玻耳兹曼公式所示，模拟系统到达一个点的概率与该点系统总能量E的指数幂成反比！！！不会疼的我们的仿真系统将大部分时间诚实地保持在低能耗的配置中，而没有进取精神. 结果是欠采样. 我们如何进行统计？这时，有些兄弟感到困惑. 高能量配置，模拟系统无法爬升，实际系统能否爬升？ . . . 是的，本构能量很高，并不是说它永远无法达到，但是概率很低，但是只要时间足够长，系统就能爬升并停留一段时间，这就是为什么有些反应会发生的原因1-2天，即中间必须要穿越的能量障碍很高，即使穿越的可能性也很低. . . 升高温度会增加交叉的可能性. . . 所有功率都取决于一个小的热运动kT = 0.6kcal / mol，从而创造了跳入龙门的鲤鱼奇迹.

您可以模拟1-2天吗？不，可悲的研究僧侣只能运行ns的仿真，因此系统无法访问反应坐标的高能垒部分. 他们应该毕业吗？这时，伞式采样就诞生了. . . 那些压迫青年的人有福了. .

由G.M. Torrie，来自加拿大的J.P. Valleau，1977年提出了一种巧妙的方法: 您的系统不愿意停留在我的反应坐标系的高能段中，不是吗？兄弟，我用绳子把你绑在那里. 具体方法是向反应物施加弹簧电势w = k *（x-x0）^ 2. 一旦反应物坐标x不在我设定的x0处，势能就会上升，您必须返回. Brother将反应坐标切成许多小段（科学名称是window），并在每个小段上像这样跑了一段时间. 这称为有偏采样，即有偏采样. 暴力方式迫使系统沿着反应坐标移动. 但是我们可以接受失真的结果吗？不. 如何将有偏抽样的结果转换回无偏抽样的概率密度？这需要一个称为WHAM的“三步布局”过程.

本文来自电脑杂谈，转载请注明本文网址：
https://www.pc-fly.com/a/dianqi/article-255615-1.html