自由能微扰

玻尔百科

核心要点

自由能微扰（FEP）基于 Zwanzig 方程，通过模拟初始态并对到受扰动的终态的能量变化的指数进行平均，从而计算自由能差。
FEP 的成功应用需要克服采样重叠不佳的问题，通常将变换过程划分为多个更小的步骤（“窗口”），并使用软核势来防止模拟失稳。
热力学循环是一种强大的策略，它使 FEP 能够计算物理性质的相对变化，例如蛋白质突变对药物结合亲和力的影响。
FEP 是一种应用广泛的通用方法，从预测生物化学中的耐药性和 pKa 偏移，到计算材料中的缺陷能，甚至量子同位素效应。

引言

我们如何预测一种新药能否有效结合其靶点？或者，一个单一的突变如何导致抗生素耐药性？科学中的这些基本问题都由自由能决定，自由能是一个涵盖了系统所有可能微观状态的物理量。然而，直接计算一个复杂系统的绝对自由能，在计算上是不可能的。本文将探讨自由能微扰（FEP），这是一种强大的计算方法，它通过专注于计算两个状态之间的自由能差，巧妙地回避了这一限制。它通过定义一条非物理的、“炼金术”般的路径，在数学上将一个状态转变为另一个状态来实现这一目标。

本文将引导您进入计算炼金术的世界。首先，在“原理与机制”一章中，我们将深入探讨 FEP 的理论核心，推导其基础——Zwanzig 方程，并探索使这些计算成为可能的实际挑战与解决方案。随后，“应用与跨学科联系”一章将展示如何利用热力学循环巧妙地应用 FEP，以解决药物设计、生物化学和材料科学中的关键问题。

原理与机制

炼金术士的梦想：计算自由能的变化

在物理与化学的宏大舞台上，一些最重要的问题都与变化有关。当药物与蛋白质结合时，能量差异是多少？一种材料的某种晶体结构比另一种稳定多少？这些问题不仅仅关乎能量，更关乎自由能，一个更精妙、更强大的概念。能量告诉你系统在单一、冻结构象下的情况，而自由能（ $F$ 代表恒定体积下的亥姆霍兹自由能，或 $G$ 代表恒定压力下的吉布斯自由能）则讲述了完整的故事。它考虑了系统可能存在的所有状态，其原子的所有摆动、振动和旋转，并按其概率加权。自由能是真正主导自发性和平衡的物理量。

困难在于，自由能是出了名地难以直接计算。它通过配分函数 $Z$ 的对数来定义，而配分函数是一个对所有可能微观状态的概率进行的令人生畏的求和。对于一个复杂系统来说，计算这个总和在计算上是不可能的。这就像试图数清世界上所有海滩上的每一粒沙子。

那么，我们该如何进行呢？我们“作弊”。或者说，我们使用一个巧妙的、感觉像作弊的技巧。我们不计算两个状态（比如状态 A 和状态 B）的绝对自由能，而是关注其差值 $\Delta F = F_B - F_A$ 。为此，我们想象一条连接这两个状态的“炼金术”路径。想象一个可以转动的旋钮，上面标有一个从 0 到 1 的参数 $\lambda$ 。当 $\lambda=0$ 时，系统处于状态 A。当你转动旋钮时，系统的性质——甚至支配其原子的法则——会平滑地转变。当你到达 $\lambda=1$ 时，系统就处于状态 B。这是一个非物理的、纯数学的构造，但它却是解开问题的关键。问题于是变成了：我们能否找到一种方法，通过监测这条神奇路径上发生的事情来计算总的自由能变化？

Zwanzig 方程：连接两个世界的桥梁

答案在于一个优美而又出奇简单的关系式，即 Zwanzig 方程，它是自由能微扰（FEP）方法的基石。它在系统的微观细节和宏观的自由能差之间架起了一座精确的理论桥梁。

让我们看看这座桥是如何搭建的。我们从基本定义出发：末态（我们称之为状态 1）和初态（状态 0）之间的自由能差与它们配分函数的比值有关：

\Delta F_{0 \to 1} = F_1 - F_0 = -k_B T \ln\left(\frac{Z_1}{Z_0}\right)

这里， $k_B$ 是玻尔兹曼常数， $T$ 是温度。状态 $i$ 的配分函数是 $Z_i = \int \exp(-\beta U_i(\mathbf{x})) d\mathbf{x}$ ，其中 $\beta = 1/(k_B T)$ ， $U_i(\mathbf{x})$ 是给定原子构象 $\mathbf{x}$ 的势能。

这个比值是棘手的部分。但请看这里。我们可以写出 $Z_1$ ，然后巧妙地乘以一：

Z_1 = \int \exp(-\beta U_1) d\mathbf{x} = \int \exp(-\beta U_1) \frac{\exp(-\beta U_0)}{\exp(-\beta U_0)} d\mathbf{x}

重新整理指数中的项：

Z_1 = \int \exp(-\beta (U_1 - U_0)) \exp(-\beta U_0) d\mathbf{x}

现在，我们除以 $Z_0$ ：

\frac{Z_1}{Z_0} = \frac{\int \exp(-\beta (U_1 - U_0)) \exp(-\beta U_0) d\mathbf{x}}{Z_0} = \int \exp(-\beta (U_1 - U_0)) \left(\frac{\exp(-\beta U_0)}{Z_0}\right) d\mathbf{x}

括号中的项 $\exp(-\beta U_0)/Z_0$ 正是当系统处于状态 0 时，在构象 $\mathbf{x}$ 处找到它的概率密度。因此，整个积分就是物理量 $\exp(-\beta (U_1 - U_0))$ 在从状态 0 采样的所有构象上的系综平均。我们用 $\langle \dots \rangle_0$ 来表示这个平均值。

将这个优雅的结果代回我们关于 $\Delta F$ 的表达式，我们就得到了 Zwanzig 方程：

\Delta F_{0 \to 1} = -k_B T \ln \left\langle \exp(-\beta [U_1(\mathbf{x}) - U_0(\mathbf{x})]) \right\rangle_0

这个方程意义深远。它告诉我们，要找到两个世界之间的自由能差，我们根本不需要探索第二个世界！我们只需要对初始世界（状态 0）进行模拟，对于我们采样的每一个构象，我们计算如果我们突然切换到第二个世界的法则，这个构象将会有的能量差 $\Delta U = U_1 - U_0$ 。然后，我们计算这个“微扰”能量差的指数平均值。这是一种强有力的重要性采样形式，我们利用一个状态的分布来测量另一个状态的性质。对于恒定压力的系统（在化学和生物学中很常见），一个类似的方程对于吉布斯自由能差 $\Delta G$ 也成立。

一个简单的例子：谐振子

为了看到 Zwanzig 方程的实际应用，让我们考虑物理学中最简单的系统之一：一个连接在弹簧上的粒子，由谐振子势 $U(x) = \frac{1}{2} c x^2$ 描述。想象我们的“炼金术”变化是使弹簧变得更硬，将其力常数从 $c_0$ 变为 $c_1$ 。这个变化的自由能代价是多少？

使用 Zwanzig 方程，我们可以精确地解决这个问题。我们会用初始弹簧常数 $c_0$ 来模拟系统。对于粒子访问的每个位置 $x$ ，我们会计算 $\Delta U = \frac{1}{2}c_1 x^2 - \frac{1}{2}c_0 x^2$ 。然后我们会计算 $\exp(-\beta \Delta U)$ 的平均值。对于这个简单的例子，积分可以被解析求解。最终的结果非常简洁：

\Delta F = \frac{1}{2} k_B T \ln\left(\frac{c_1}{c_0}\right)

这个结果完全符合直觉。使弹簧变硬（如果 $c_1 > c_0$ ）需要付出自由能，并且这个变化与温度成正比。在绝对零度下，自由能的变化为零（如果我们忽略量子效应），因为粒子只会停在势阱的底部。在更高的温度下，粒子会探索更广的位置范围，用更硬的弹簧限制它所带来的熵代价变得更加显著。

微扰的陷阱：当世界不重叠时

Zwanzig 方程在理论上是精确的，但在实践中，它隐藏着一个险恶的陷阱。该方法依赖于对状态 0 的模拟能为状态 1 的重要构象提供良好的采样。这就是相空间重叠的概念。如果两个状态非常不同，它们的重要构象可能根本不重叠，FEP 可能会惨败。

想象一下，你生活在一个永远阳光明媚的热带气候中（状态 0），你想估算西伯利亚人穿的冬衣的平均厚度（状态 1）。你可以调查你的邻居，问他们：“如果气温突然降到 $-40^\circ$ C，你会穿上什么外套？”你的大多数邻居都穿着短裤和 T 恤，他们完全没有参考框架。他们的回答将毫无意义。然而，你可能会找到一个刚从雅库茨克搬来的人。他拥有一件合适的毛皮大衣。当你对结果进行平均时，他那唯一的、合理的回答与其他所有人的回答相比会显得如此极端，以至于它将完全主导你的整个估算。你的最终平均值将极度不准确，并且如果你碰巧找到第二个西伯利亚移民，这个平均值会发生剧烈变化。

这正是 FEP 在重叠不佳时遇到的问题。对于状态 1 最重要的构象（例如，药物的低能量结合姿态）在状态 0 的模拟（药物未结合，漂浮在水中）中极为罕见。当状态 0 的模拟因统计上的偶然性而采样到这些罕见的、类似状态 1 的构象之一时，能量差 $\Delta U = U_1 - U_0$ 将是一个很大的负数。权重因子 $\exp(-\beta \Delta U)$ 将会是一个天文数字。整个平均值被这些罕见但权重极高的事件所主导，导致巨大的统计误差和非常不可靠的估计。这种高方差在形式上是由权重分布的二阶矩驱动的，它对这些罕见事件更为敏感。

这个问题的一个明显症状是滞后现象：计算出的正向自由能变化 $\Delta G_{A \to B}$ 不等于反向变化的负值 $-\Delta G_{B \to A}$ 。这表明模拟没有正确收敛，结果不可信。

实用的魔法：让 FEP 发挥作用

幸运的是，我们有一系列技巧来驯服重叠不佳这头猛兽。指导原则很简单：如果两个状态之间的跳跃太大，就采取更小的步长。

分层（窗口化）： 我们不试图一次性完成从 $\lambda=0$ 到 $\lambda=1$ 的整个炼金术转变，而是将路径分解成许多小的、中间的步骤或“窗口”。例如，我们可以使用 $\lambda = \{0.0, 0.1, 0.2, \dots, 1.0\}$ 。然后我们计算每一步小的自由能变化 $\Delta F_{\lambda_i \to \lambda_{i+1}}$ ，此时相空间重叠要好得多。总的自由能变化就是所有窗口变化的总和。我们甚至可以更聪明地在系统变化最快的区域放置更多的窗口，将我们的计算预算分配到最需要的地方。

软核势： 当我们“创建”或“湮灭”一个原子时，会出现一个特别棘手的问题，这在计算溶剂化或结合自由能时是常见操作。想象一个原子在 $\lambda=0$ 时从真空中出现。当我们慢慢开启它的相互作用时，如果环境中的另一个原子碰巧进入了同一空间，会发生什么？使用像 Lennard-Jones 势这样的标准势，排斥能以 $1/r^{12}$ 的形式缩放，当原子间距离 $r$ 趋于零时，这将变为无穷大。这种“端点灾难”会产生无穷大的力，并使模拟崩溃。

解决方案是使用软核势。可以这样想：原子不是以一个坚硬、不可穿透的球体形式出现，而是首先以一个其他原子可以无碍穿过的“幽灵”形式出现。随着 $\lambda$ 的增加，这个幽灵慢慢固化成一个正常的原子。在数学上，这些势被修改，使得即使在任何中间 $\lambda$ 值的零距离处，能量也保持有限，从而防止了灾难性的发散，并允许模拟平稳进行。

微扰之外：其他工具一瞥

FEP 是一个强大的工具，但它不是计算炼金术士武库中唯一的工具。另外两种主要方法值得了解。

热力学积分（TI）： TI 不计算离散窗口之间的自由能差，而是采用一种基于微积分的方法。它在几个 $\lambda$ 点计算自由能曲线的斜率 $\langle \partial U / \partial \lambda \rangle_\lambda$ 。然后通过将此斜率从 $\lambda=0$ 积分到 $\lambda=1$ 来找到总的自由能差。这类似于通过积分路径的陡峭程度来计算徒步旅行中总的海拔变化。

Bennett 接受率（BAR）： BAR 方法是一种在统计上更复杂的方法。它认识到运行一个正向模拟（ $A \to B$ ）和一个反向模拟（ $B \to A$ ）可以为你提供两组数据来估计同一个量。BAR 不是分开处理它们，而是以一种统计上最优的方式结合来自两个方向的数据。通过求解一个自洽方程，它为 $\Delta G$ 产生一个单一的、低方差的估计，这比单向 FEP 可靠得多，尤其是在重叠不佳的情况下。因此，BAR 及其多态泛化形式（MBAR）通常被认为是自由能计算的“金标准”。

通过这些原理和机制，曾经是炼金术士的梦想——将一种物质嬗变为另一种——已经成为一种常规的、尽管具有挑战性的计算工具。通过使用 FEP、TI 和 BAR 等方法小心地驾驭数学领域，我们可以准确地预测驱动自然界基本过程的自由能变化。

应用与跨学科联系

在掌握了自由能微扰的原理之后，我们可能会觉得它只是统计力学中一个优美但相当抽象的部分。但事实远非如此。当我们将这套理论机制应用于现实世界中那些混乱、复杂而又引人入胜的问题时，真正的魔力才开始显现。自由能微扰不仅仅是一个方程；它是一台计算显微镜，一种“计算炼金术”，让我们能够窥探原子世界，并预测那些塑造我们生活的宏观性质，从药物的疗效到新材料的特性。这是一场想象力的旅程，我们巧妙地设计出非物理的“炼金术”路径，来解决非常真实的物理问题。

真理的试金石

在我们开始探索未知之前，科学家的职责是在他们完全理解的事物上测试他们的工具。对于自由能计算而言，不起眼的谐振子——物理学家用来模拟任何振动物体（从固体中的原子到弹簧上的质量）的模型——提供了一个完美的试验场。我们可以用完美的数学精度计算两个不同刚度（比如 $c_0$ 和 $c_1$ ）的弹簧之间的自由能差。精确的答案出奇地简单： $\Delta F = \frac{1}{2}k_B T \ln(c_1/c_0)$ 。

通过尝试用自由能微扰重现这个已知答案，我们了解了该方法的实际陷阱。想象一下，试图通过从一个“软弹簧”状态采样构象来计算“硬弹簧”状态的自由能。软弹簧探索的位置范围很广，但硬弹簧主要局限于中心附近的一个狭窄区域。我们的 FEP 计算依赖于软弹簧偶然采样到那些对硬弹簧来说很重要的、被压缩的稀有构象。如果刚度 $c_1$ 和 $c_0$ 相差太大，它们重要空间区域之间的“重叠”几乎为零。我们的计算将被稀有事件主导，并且无法收敛。这就像试图通过研究一只老鼠来了解一头大象，寄希望于老鼠会漫步到大象的脚印里。

这个简单的例子教给我们一个深刻的教训：FEP 中的“微扰”必须是温和的。在实践中，这意味着我们不能直接从状态 A 跳到状态 B。相反，我们必须构建一条由许多小的、中间步骤组成的平缓路径，用参数 $\lambda$ 索引，并将每一步小的自由能变化相加。为了防止原子在这条路径上出现或消失时相互碰撞，我们使用巧妙的“软核”势，这些势在短距离内表现良好。这种谨慎的、循序渐进的方法，正是使我们的计算炼金术成为可能和可靠的原因。

循环的艺术：为生物学破解热力学

也许自由能微扰最引人注目的应用是在生物化学和药物设计领域。核心挑战通常是预测一个微小的变化——比如蛋白质中的一个突变，或对药物分子的一个修饰——将如何影响其功能。直接模拟药物与蛋白质结合以计算结合自由能 $\Delta G_{\text{bind}}$ ，在计算上是极其庞大的。但我们通常不关心绝对的结合能，而是关心结合能的变化 $\Delta\Delta G_{\text{bind}}$ 。这正是该方法的天才之处，通过使用热力学循环而大放异彩。

想象一下，我们想知道酶中的一个突变如何影响其与抗生素的结合能力。这是理解和对抗抗生素耐药性的核心问题。假设我们有野生型酶 ( $P_{WT}$ ) 和一个突变体 ( $P_{MUT}$ )。我们想知道它们与抗生素 ( $L$ ) 的结合亲和力之差，即 $\Delta\Delta G_{\text{bind}} = \Delta G_{\text{bind}}^{\text{MUT}} - \Delta G_{\text{bind}}^{\text{WT}}$ 。直接计算任何一个 $\Delta G_{\text{bind}}$ 都很困难。

但由于自由能是一个状态函数——意味着变化只取决于起点和终点，而与路径无关——我们可以构建一个巧妙的弯路。考虑以下四种状态：

未结合的野生型蛋白 ( $P_{WT}$ )
已结合的野生型蛋白 ( $P_{WT}L$ )
未结合的突变蛋白 ( $P_{MUT}$ )
已结合的突变蛋白 ( $P_{MUT}L$ )

物理过程是水平路径：配体与野生型和突变体的结合。炼金术的、非物理的过程是垂直路径：在未结合（apo）状态和已结合（holo）状态下，神奇地将野生型蛋白转变为突变体。围绕这个闭合回路的总自由能变化必须为零。这给了我们一个惊人地简单的方程：

$\Delta G_{\text{bind}}^{\text{WT}} + \Delta G_{\text{mut, holo}} = \Delta G_{\text{mut, apo}} + \Delta G_{\text{bind}}^{\text{MUT}}$

重新整理这个方程，我们便得到了我们想要的东西：

$\Delta\Delta G_{\text{bind}} = \Delta G_{\text{bind}}^{\text{MUT}} - \Delta G_{\text{bind}}^{\text{WT}} = \Delta G_{\text{mut, holo}} - \Delta G_{\text{mut, apo}}$

这是一个优美的结果。我们通过计算两个非物理的炼金术转变的自由能，找到了物理结合能的变化！如果计算出的 $\Delta\Delta G_{\text{bind}}$ 是正值，意味着突变体中的结合变弱了，这可能是耐药性的一种机制。如果是负值，则突变增强了结合，这或许指向一个更有效的药物设计方向。这一思想是现代计算药物发现的基石。

同样的“循环技巧”使我们能够解决其他基本的生物化学问题。例如，氨基酸残基的酸度，即其 $\mathrm{p}K_{\mathrm{a}}$ 值，会受到其在折叠蛋白质中局部环境的显著影响。为了计算这种偏移，我们再次使用一个循环。我们计算在蛋白质内部使该残基去质子化的炼金术自由能，然后计算在体相水中对一个小的模型化合物（如天冬氨酸残基对应的乙酸）进行相同转变的自由能。这两个炼金术计算之间的差异给了我们去质子化自由能的差异，从中我们可以找到 $\mathrm{p}K_{\mathrm{a}}$ 的偏移。这种方法非常强大，因为它巧妙地回避了计算单个质子的绝对溶剂化自由能这个 notoriously difficult 的问题。当我们研究的反应，如去质子化，涉及共价键的断裂和形成时，我们甚至可以将 FEP 与量子力学结合在一种混合的 QM/MM 方案中，用量子力学方法处理反应部分，而蛋白质环境则保持经典力学处理。

从生物学到材料科学：FEP 的更广阔天地

这种“炼金术循环”思维的力量远远超出了生物学。它是统计力学的一个通用工具。考虑化学中的一个基本性质：溶解度。一个给定的分子会溶于水吗？如果会，溶解度是多少？这取决于标准溶解自由能 $\Delta G_{\text{sol}}^{\circ}$ ，即一个分子从其纯晶相转移到溶液中的自由能变化。

再一次，热力学循环为我们提供了帮助。我们可以间接地计算 $\Delta G_{\text{sol}}^{\circ}$ ，通过想象我们首先将晶体中的分子变成气体（升华， $\Delta G_{\text{sub}}^{\circ}$ ），然后将气体分子溶解到水中（水合， $\Delta G_{\text{hyd}}^{\circ}$ ）。我们可以用 FEP 通过在水中“炼金术”般地“消失”分子来计算 $\Delta G_{\text{hyd}}^{\circ}$ 。或者，我们可以直接通过将晶体中的一个分子炼金术般地转变为一个不相互作用的“幽灵”粒子，并对水中的一个分子做同样的操作来实现。这两个炼金术过程的自由能之差就给出了我们想要的 $\Delta G_{\text{sol}}^{\circ}$ 。这些方法在从药理学到环境科学等领域都具有不可估量的价值。

这种思维方式甚至让我们能够理解固体材料的结构。完美的晶体是一个有用的理想化模型，但真实的材料是由它们的缺陷——空位、杂质和其他缺陷——来定义的。FEP 允许我们计算在晶格中创建一个单一缺陷的自由能代价。但在这里，我们遇到了统计力学中另一个微妙而优美的思想。如果我们计算在某个特定原子位点创建一个缺陷的自由能，我们还没有完成。如果晶体有 $N$ 个相同的、对称等价的位点可以形成这个缺陷，那么带有一个“缺陷”的系统的真实自由能必须更低。为什么？因为系统有更高的熵——选择在哪里放置缺陷的“选择的熵”。这导致了一个对称性校正项 $-k_B T \ln N$ ，我们必须将其加到我们的单位点 FEP 结果中。这是一个引人注目的提醒：自由能不仅仅关乎能量，还关乎计算事物存在的不同方式的数量。

量子前沿：权衡同位素效应

为了真正领略自由能微扰的惊人范围，考虑最后一个令人脑洞大开的应用：动力学同位素效应（KIE）。当一个氢原子被其更重的同位素氘取代时，涉及该原子的化学反应速率通常会发生变化。这是一种纯粹的量子力学效应，由零点振动能和隧穿效应的差异驱动。

我们经典的 FEP 框架如何可能捕捉到这一点呢？通过将其与一种称为路径积分分子动力学（PIMD）的技术相结合，该技术将每个量子粒子表示为一个由经典“珠子”组成的环，我们可以计算量子统计性质。现在，是展现炼金术天才一击的时刻：我们定义一条路径，不是将一种元素变为另一种，而是缓慢地将一个原子的质量从氢的质量变为氘的质量。我们进行两次这种质量嬗变的 FEP 计算：一次是针对反应物，另一次是针对反应的过渡态。这些量子自由能变化的差异给了我们氢和氘反应活化能垒的差异，从中我们可以计算出 KIE。

请仔细思考一下。我们正在使用“计算炼金术”在计算机内部改变自然界的一个基本常数——原子核的质量——来预测化学反应速率上的一个微妙的量子效应。这是对统计力学原理的力量和美的终极证明，表明只要有一条巧妙的路径和对基本原理的深刻理解，几乎没有我们不敢 tackling 的问题。