Umbrella Sampling

自由能差是多数物理化学过程的驱动力，例如化学反应、蛋白质与小分子的结合等过程的发生与否和发生速率都由相应的自由能差决定。因此自由能差的准确计算是计算化学和计算生物学领域的重要挑战之一。自由能差可以通过分子模拟方法计算，不过通常的分子模拟方法在处理复杂体系时都会面临着采样效率低下的问题。这是由于在模拟过程中系统可能长期困于一个亚稳态而无法探索相空间更广阔的其他区域。所以怎样能够高效地探索系统在相空间的概率分布全貌是解决这个困境的关键所在。伞形采样（Umbrella Sampling）是一种著名的、成熟的且被广泛接受的增强采样方法。它借助偏置势（bias potential）来改变自由能景观（free energylandscape）从而驱动模拟探索任何想要探索的地方，就像被一把伞保护着一样，不管地形、天气如何都风雨无阻。

——简介——

自由能差的计算是计算科学的一个核心任务，不管热力学性质还是动力学性质都与它息息相关。它的应用范围涉及从固态物理、催化反应、生化过程到理性的药物设计等多个领域。以目前的计算机的算力，一般情况下，我们只能在经典统计力学的框架下近似地获取它。

图1  沿着反应坐标的自由能景观。

在实际中，我们通常关心某个系统沿着一个或多个自由度的反应坐标（reaction coordinate ξ，物理学家更喜欢称其为collective variable）的行为，因此探索沿着反应坐标的自由能景观更为重要。

注：反应坐标是相点x的连续函数Ξ(x)，它可以区分需要研究的两个或多个热力学状态。

根据统计力学，我们知道自由能是与概率密度相关：

这里的P(ξ)是随机变量子集ξ的无条件概率密度函数：

动力学模拟的目的就是采样这个分布，从而数值地计算出我们想要的热力学量和动力学量。但是对于复杂的体系，它通常存在许多亚稳态，这导致在有限的模拟时间内，某个自由能极小值附近相空间被充分地采样，但是像自由能垒（显著高于k_BT）这样的稀有事件（rare event）基本很少被采到。想要获取自由能P(ξ)的全貌，稀有事件的采样是必需的。这时候就需要借助一些增强采样计算来帮忙了。

伞形采样就是这样一种增强采样技术。它的通常做法是首先将反应坐标分为多个窗口（图2），每个窗口都独立地采样，最后不同窗口的采样结果将被合并为完整的自由能轮廓。形象地说，就是一开始只有一个人在漫无目的地探索世界，群山环绕，可能一辈子也只能呆在这个小山村；多年之后，现在技术发达了，有了N架动力十足的无人机，然后它们能够毫不费力地被一口气（不管海拔多高）传送到了世界各地，每架就重点观察它附近的区域，最后将它们的影像整合一下，世界的壮丽景色就尽收眼底。

图2  反应坐标被分为多个窗口，每个窗口仅包含一小段反应坐标。

不过要怎么做到无人机能够准确被传送到各个地点且不会被大风刮走呢？也就是我们要怎么保证每个独立的模拟只重点采它负责的窗口附近的样点呢？其实很简单，只要分别加一个额外的偏置势能将其束缚在某个位置就行（图3），这样它就可以负责重点采样这个区域。

图3  每个窗口添加了偏置势，用于约束采样范围。

每个窗口在外加偏置势后分别采用常规的分子模拟方法进行采样，我们就可以得到每个窗口对应的偏置概率分布（图4）。偏置势选得适当，每个窗口的概率分布应该都是单模的（unimodular），这意味着偏置势的确约束住了采样范围。

图4  外加偏置势后的每个窗口对应的偏置概率密度。

偏置概率分布与原本的无偏概率分布有着明确的函数关系，所以我们可以通过模拟得到的偏置概率分布重新计算无偏概率分布。其实除了需要处理一个由偏置势引起的自由能项，由单个窗口的偏置概率分布就可以准确地恢复全反应坐标范围的无偏概率分布。但是为了获得尽可能准的全局无偏概率分布，最好还是将所有窗口的模拟结果都结合在一起，并且结合过程尽可能采取最终的无偏概率分布统计误差最小的方式。目前常用的结合伞形采样模拟结果的方法有加权柱状统计图分析法（Weighted Histogram Analysis Method, WHAM）和伞形积分（Umbrella Integration， UI）。

以上就是伞形采样的基本概念，详细数学推导将在下一节呈现。

——理论——

无偏与偏置的关系

前面已经介绍了基本概念，那么这里就从如何通过偏置模拟结果恢复无偏自由能开始吧。首先，偏置的哈密顿是这样的：

其中，每个窗口的偏置势B_i是仅依赖于反应坐标Ξ(x)的额外能量项。那么这个偏置的体系对应的概率分布就是这样的：

由于偏置势仅依赖于反应坐标，它在分子的积分中可以作为常数提取出来，

将它代入无偏概率分布的表达式中即可得到：

根据上式得到的无偏概率分布就很容易可以知道无偏自由能了：

式中，F_i=-(1/β)ln<exp[-βB_i(Ξ(x))]>是一个由偏置势引起的类自由能项，它与反应坐标ξ无关。如果我们只取一个窗口，并且这个窗口横跨整个反应坐标范围，那么对于仅关心自由能差的情况而言，F是无关紧要的。

自适应伞形采样

到此为止，我们的偏置势是未知的，也就是说以上推导适用于任何偏置势。要想推导继续下去，我们此时就要引入偏置势的一个具体的表达式。那就需要考虑偏置势怎么取的问题了。我们先想想引入偏置势的目的到底是什么，没错，就是为了增强采样。所以偏置势的引入应该使得偏置概率分布越平坦越好。

此时我们再回过头来看一下偏置概率分布：

其中，C是一个常数。因此，如果B_i(ξ)越接近-A(ξ)，那么偏置概率分布就会越平坦。但是实际上我们压根就不知道自由能的形貌，所以偏置势的选取还是令人头大。

虽然我们不知道自由能景观是什么样的，但是随着分子模拟的进行我们不就可以粗浅地知道了嘛。没错，这就是自适应伞形采样技术（Adaptive Umbrella Sampling）的原理。不知道自由能的形貌也不妨碍我们随便给个偏置势，只要有了这个偏置势（或者没有也行）就可以跑分子模拟了，采样了一段时间后就可以算出一个自由能景观，然后根据它更新偏置势，再接着跑，如此循环往复，偏置势就会越来越接近自由能景观的负数，偏置概率分布也就越来越平，采样就会越来越容易，最后我们就能采到统计误差非常小的分布。

普通的伞形采样

当然，这不是我们在简介中提到的伞形采样的做法。自适应伞形采样以逐渐夷平山峰的方法来探索更广阔的世界，而普通的伞形采样则是以大部队空降的方法看到了世界的全貌。

既然偏置势未知，那也可选取简单的函数试试啊，比如说，简谐势：

只要取得足够大，我们几乎可以将一次模拟困在窗口附近，重点观察每个窗口附近的自由能形貌。但是又不能取得太大，太大了你就会连自己窗口内的部分都看不全。而且WHAM这种叠加技术也要求每个窗口的分布必须有所重叠，所以k_i的选择很重要。

加权柱状统计图分析法

最后我们得到了每个窗口对应的无偏分布，我们要将它们结合在一起，尽量获取一个最佳的分布。先看看加权柱状统计图分析法（WHAM）吧。该方法的做法，正如其名称所示，就是通过各个窗口的分布的加权平均来计算全局分布的：

权值p_i(ξ)以使得P(ξ)的统计误差最小的方式选取，即满足：

再加上Σp_i(ξ) = 1的条件就得到了：

式中，N_i是窗口i采样的总步数。而F_i可以通过下式计算：

这里我们看到F_i又依赖与P(ξ)。因此P(ξ)可以通过不停地迭代上述几个式子至收敛来求解。

伞形积分

另一种结合各个窗口模拟结果的方法是伞形积分（UI）。它通过加权平均平均力（mean force，也就是自由能相对于反应坐标的一阶导数）而非概率分布来绕过了计算F_i的问题。平均力与偏置分布的关系如下：

为了使推导继续下去，此时需要对偏置概率分布做一步近似，即执行累积量展开，并截断到第二项：

如果窗口很小，这种做法是合理的。因此，平均力可以简化为：

然后积分它，就可以得到自由能的形貌：

当然，为了获得更好的自由能景观，我们需要整合所有窗口的结合以获取统计误差最小的解：

最后将这个全局的平均力积分就可以得到全局的自由能景观了。

理论部分到此结束，想看更多的变体版本，请查阅更多的文献吧。

——总结——

加速采样技术有很多种，什么时候应该使用伞状采样的问题无法笼统地回答，这主要取决于特定的系统，一些此领域专家已经比较了一些增强采样方法的适用性，详情可以参考这篇文章的参考文献的参考文献。每种增强采样方法都很有趣，所以呢，如果有时间的话，我打算未来出一系列的增强采样公众号文章。

总之，伞形采样是一种成熟且被广泛接受的计算自由能差的方法，欢迎尝试哦。

参考文献

Kästner,J. Umbrella sampling. Wiley Interdiscip. Rev. Comput. Mol. Sci. 1, 932–942(2011).

点击左下角的"阅读原文"即可查看原文章。

---

作者：林豪禹

审稿：李亦博

编辑：卞薇洁

GoDesign

ID：Molecular_Design_Lab

（ 扫描下方二维码可以订阅哦！）

本文为GoDesign原创编译，如需转载，请在公众号后台留言。