自由能计算-讲解

自由能的计算在计算科学领域是一个非常重要的研究课题, 所有物理、化学以及生物过程自发进行达到的程度都是由初态(反应物)与终态(产物)之间的自由能变化决定的。

对于一个反应来说, 反应的自由能差对应于该反应的热力学上的可逆功。

目前, 自由能的计算方法有多种类型, 经典的方法有 Umbrella Sampling、Free Energy Perturbation和 Thermodynamic Integration等, 这类方法原理上比较严谨, 

计算结果也比较精确, 但是需要大量的数据采集, 所需的计算资源很多, 在实际应用上受到体系大小以及计算资源的限制。

最新发展出来的自由能计算方法有 Jarzynski Equality、Metadynamics等,

Jarzynski Equality 由大量不可逆功的系统平均值来求解反应过程的自由能数据, 但是在具体使用时, 此方法的准确度受研究体系的影响很大

Metadynamics 则是不断地对研究体系施加额外的高斯型排斥势能函数, 使目标分子能够脱离自由能的势阱, 访问到高自由能区域, 

进而通过施加的排斥势能函数来求解自由能数据, 在求解自由能二维势能面上有很大的优势.

在生物学领域中, 生物大分子之间的相互作用、蛋白质折叠、离子通道以及生物分子与材料的相互作用等都是很重要的研究课题。

在这些研究中, 对生物大分子结构与功能, 和对新设计的生物分子或生物材料的评价是很重要, 都可以通过对其能量学的考察得到揭示。

然而, 由于上述各种自由能计算方法的原理不同, 它们在具体操作的步骤, 计算自由能的效率, 以及对研究体系的适用性等方面都有各自的特点, 

因此在计算体系自由能数据的时候, 需要根据研究体系的不同, 计算资源的多少来选定特定的计算方法. 下面我们对各个方法的原理、计算效率以及在生物大分子体系中的适用性等方面进行阐述, 

使读者对其能有一个较为全面的了解。

Umbrella Sampling

根据统计热力学的基本原理可知, 系统的自由能与研究体系的概率分布有关. 对于系统从状态 0 转变为状态 1 的过程, 其自由能差值可以由研究系统在两个状态下的概率分布得到:

因此对于一些特殊体系的自由能数据, 我们可以统计出体系在所要研究两种状态的概率分布之比, 进而使用公式(1)来计算相关的自由能数据. 在实际应用中, 人们往往对自由能数据随着给定的反应坐标 ξ 的变化趋势比较关注, 这种自由能随着某个反应坐标变化的曲线称为平均力势(potential of mean force, PMF). 例如, 计算体相水中两个水分子分离的自由能曲线时, 可以通过分子动力学模拟, 得到两个水分子直接的径向分布函数 (radial distribution function, RDFξ), 然后即可通过公式(2)得到相关的PMF 数据:

其中 ξ 为反应坐标, 即水分子分离的距离, RDF0 为选择的自由能参考点的数据. 因此对于概率分布求解简单的体系, 可以直接通过概率分布求解相关自由能数据. 然而在实际体系中, 由于受模拟时间的限制, 研究体系(目标分子)绝大多数的时间都处在自由能低的区域, 一般的动力学模拟很难采集到其在自由能高的区域的分布, 无法利用公式(1)求解相关自由能数据. Umbrella Sampling 就是为了解决高自由能区域的取样问题而提出来的.

1977 年, Umbrella Sampling 最早由 Torrie 和Valleau 提出. 在这个方法中, 为了解决普通分子动力学在高自由能区域的取样问题, 整个反应路径被分为多个区域, 这些区域被称为 Umbrella Sampling 的窗口. 在每一个窗口 i 中, 一个额外的势能项被添加到研究体系的哈密顿量中:

其中 x 和 p 分别表示研究体系的位置和动量, 而额外加的项为一个静态的谐振势, 它是反应坐标r= r(x)的函数, 对于给定的中心位置 r0, 则有:

谐振势 hi(r)的作用是把研究体系的位置限制在中心位置 r0 附近, 从而加强研究体系在这段区域的取样. 整个反应路径被分割成许多段很小的区域, 这样的区域在 Umbrella Sampling 中称为窗口. 对每一个窗口, 系统都被谐振势限制在对应的位置 r0. 在对应的谐振势的限制下, 通过分子动力学模拟, 我们就能够得到研究体系在反应路径上的一系列的有偏差的概率分布 ρi(r). 需要指出的是, 这里相邻的两个窗口的概率分布之间要有足够重叠的部分, 才能得到准确的自由能曲线. 然后这一系列的概率分布 ρi(r)再通过权重直方分析法(weighted histogram analysis method, WHAM)组合起来, 去除掉额外添加 的谐振势的影响, 从而得到研究体系沿着反应坐 标 ξ 的准确的概率分布, 进而求解得相关的自由能曲线. 

由于 Umbrella Sampling 计算自由能具有严格的理论依据, 因此可以得到非常准确的自由能数据. 在生物体系中, 对于 Umbrella Sampling 来说, 溶于水溶液的生物大分子体系是非常大的 , 所需的大量的计算资源是 Umbrella Sampling 在该研究领域推广的最大阻力. 此外, Umbrella Sampling 还要求相邻窗口之间拥有足够的重叠部分. 对于很多体系而言, 某些自由能较高的区域需要很多次的尝试, 才能得到较好的取样, 这不仅增加了所需的计算资源, 拖延了计算周期, 更降低了操作的便捷性. 因而, 在生物大分子领域, 相对于在实际研究体系中的应用, Umbrella Sampling 更多地应用于对新提出的自由能计算方法准确性的检验 上.

Free Energy Perturbation

Free Energy Perturbation 由 Zwanzig 于 1954 年提出将研究体系的两个状态(初态和终态)的自由能差与其能量差值的系统平均联系到了一起:

这样研究体系初态和终态之间的自由能差值就能由干扰势能的系统平均值计算得到. 这种方法的有效性已经被很多体系所证实, 尤其是波动服从高斯概率分布体系。上述公式只能用来计算初态和终态之间只有较小差异的研究体系, 即干扰势能V很小的体系. 而在实际应用过程中, 所计算自由能差的初态和终态往往有比较大的差异, 这就使得上述公式不能直接应用到研究体系中. 解决方法是在初态和终态之间加入一系列的中间态: 引入耦合变量λ, 通过微调λ人为制造出一系列只有微小差异的中间态, 每一个中间态的能量为初态和终态的耦合:

然后通过公式(7)求解相邻两个中间态的自由能差值∆Gi, 最后将所有中间态之间的自由能差值加和, 就得到初态和终态之间的自由能差:

Free Energy Perturbation 计算自由能的基本思想就是从一个已知的状态出发, 通过一系列的微小变化, 使之达到终态, 在每一个中间态上进行分子动力学模拟, 得到体系的势能变化, 最终求解得到初态和终态的自由能差. 根据λ的间隔大小, 可以将 Free Energy Perturbation 分为三种方法: 固定步长, 慢增长以及动态步长. 在固定步长方法中, 相邻两个中间态的间隔是固定的. 在慢增长方法中, 中间态的数据趋近于无穷多个, 相邻两个中间态之间的差异接近于零. 在动态步长方法中, 每一步微扰的步长是可变的, 系统根据上一次微扰的自由能数据变化大小确定下一次微扰的步长.根据热力学原理构建热力学循环以后, Free Energy Perturbation 就可以用来计算相关反应过程的自由能变化, 对于一般的研究体系, 需要使用 Free Energy Perturbation 进行两次自由能计算, 再根据构建好的热力学循环, 得到最终的自由能变化. 然而, 由于需要对一系列的中间态进行取样, Free Energy Perturbation 计算自由能的过程就需要消耗大量的计算资源, 在计算效率上与其他方法有较大的差距. 由于计算过程只涉及到初态和终态, 中间的状态是人为设置的, 所以最终结果只有一个自由能变化的数据, 不像其他方法那样具有反应坐标以及对应的自由能数据变化趋势. 因此对于小分子体系来说, 使用Free Energy Perturbation 方法可以得到较为准确的自由能数据, 但是涉及到生物大分子体系时, 此方法会给计算带来非常大的压力, 可以截取关键片段进行相关计算, 或者使用其他计算效率相对较高的方法.

Thermodynamic Integration

Thermodynamic Integration 的基本计算公式是由Kirwood 提出的. 在这个方法中, 与 Free Energy Perturbation 比较类似的是, 研究系统的两个状态的自由能变化也由一个耦合变量λ连接, 初态时λ= 0, 终态时λ= 1. 通过连续改变λ, 得到一系列介于初态和终态的中间态, 通过对这些中间态的取样, 最终计算得到自由能变化数据. 以 H 表示研究体系在某个状态λ下的能量, 则其具体的自由能计算公式为:

在实际操作中, 自由能能量的动力学部分并没有被包含在这里面, 这是因为在计算中, 可以根据相关基本原理明确地得到动力学部分的解析解, 同时在利用热力学循环求解研究体系自由能变化时, 动力学的贡献会相互抵消. 在 Thermodynamic Integration 方法中, 计算自由能的基本思想是研究系统在无限长时间内的能量均值等价于其系综平均值, 因而可以通过有限时间的取样来估计系综平均值, 进而通过上述公式求解得到相关自由能数据. 由于在实际应用中, 无法得到λ连续变化的相关数据, 只能为λ设定一些差异较小的离散值. 通过对研究体系在不同λ对应的状态下的多次取样, 估计出研究体系的能量以后, 由下面的自由能计算公式, 就可得到相关的自由能变化数据:

其中 N 为计算过程中涉及到的研究体系状态的个数, 也即λ的个数; nλ为对研究体系的某个状态λ的取样次数; Vλ则是研究体系在状态λ下的势能; 坐标 ri, k 则与研究体系在状态λ下的玻尔兹曼概率相关, 其中下标i 和 k 分别指对研究体系每个状态的取样数以及涉及到的研究体系状态的总个数. 尽管 Thermodynamic Integration 计算自由能的相关理论出现得很早, 其准确、通用的计算方法却是在最近才提出的. Thermodynamic Integration 方法与 Free Energy Perturbation 原理以及操作流程上比较类似, Thermodynamic Integration 方法在原理上比较严谨, 计算所得结果较为精确, 计算相关研究体系自由能变化的时候需要构建热力学循环, 通过两次自由能数据的计算, 可以求得自由能变化的数据. 同样, Thermodynamic Integration 也需要对研究体系的各个状态进行大量的样本采集, 需要消耗大量的机时. 因而对于较为简单的研究体系来说, 采用 Thermodynamic Integration 方法计算自由能数据的可以得到很好的结果, 而对生物大分子体系来说, 对计算资源的需求较大.

Jarzynski Equality

Jarzynski Equality 是由 Jarzynski 于 1997 年提出的. 与 Umbrella Sampling 类似的是, 在使用Jarzynski Equality 计算自由能的时候, 也要用到额外的谐振势引导研究体系访问高自由能区域, 以采集一般分子动力学模拟中研究体系无法访问的高自由能区域的样本. 而与 Umbrella Sampling 方法不同的是, 外加的谐振势不是静止的, 而是随时间做恒速运动的. 该谐振势的具体表述如下:

其中, μ0 为谐振势的初始位置, v 为谐振势的运动速度. 这种谐振势做恒速运动的模拟称为操控式分子动力学模拟(Steered Molecular Dynamics Simulation). 操控式分子动力学模拟是受原子力显微镜中使用的单分子操纵的启发而开发的, 是单分子操纵对应的理论模拟. 在操控式分子动力学模拟中, 谐振势引导研究体系从自由能低的区域访问自由能高的区域时, 会对研究体系做功. 由于自由能差是反应过程的热力学上的可逆功, 因此操控式分子动力学中谐振势做的不可逆功 W 就是该反应过程自由能差 ∆G 的上限, 即:

Jarzynski 通过大量的理论模拟, 证实了在有足够多的不可逆功Wi的取样时, 反应过程的自由能差∆G与不可逆功 Wi 有如下关系, 即 Jarzynski Equality:

其中 β 为温度因子.这个公式的成立很大程度上依赖于一些不可逆功接近自由能差时的取样, 而这类至关重要的取样是随着操控式分子动力学的速度增长呈指数形式递减的. 因此从理论上讲, 在使用 Jarzynski Equality 计算自由能数据的时候, 不可逆功数据的采样越多, 得到的结果就越准确. 在实际应用中, 尽管理论上自由能差与操控式分子动力学的速度无关, Jarzynski Equality 的可靠性仍仅限于不可逆功的标准偏差在kBT 的数量级上的慢速过程. 而且在计算自由能数据的时候, 需要一系列的不可逆功的数据, 不可逆功的数目要足够多才能得到可靠的结果.

 Jarzynski Equality 提出十几年以来, 大量的研究工作说明 Jarzynski Equality 在求解简单体系自由能数据的时候, 能够得到与 Umbrella Sampling 准确度相当的结果. 并且, Jarzynski Equality 在操作的便捷性与计算效率上面都较 Umbrella Sampling有更好的表现. 然而在体系稍微复杂一些的时候, 功的耗散以及外加谐振势对生物大分子构型的干扰都会对所得自由能数据的准确度产生较大的影响。

Metadynamics

Metadynamics 是 Parrinello 等于 2002 年提出的一个计算多维自由能势能面的方法. 这个方法与Umbrella Sampling 以及 Jarzynski Equality 一样, 都是利用外加额外的势能, 能够更快地对研究体系的高自由能区域取样 . 与前两种方法不同的是 , Metadynamics 在计算自由能的时候, 对研究体系所加的额外的势能为高斯型的排斥势能. Metadynamics方法类似于“添坑”的过程, 通过对反应路径上的低自由能区域持续地添加排斥势能, 阻止研究体系再次访问反应坐标上已经采集过数据的区域, 不断地促使研究体系向高自由能区域运动, 最后再根据所添加的一系列排斥势能来计算相关的自由能数据. 这个方法最初应用在粗粒度模拟中研究相关体系在变量空间中的动力学行为. 之后, Metadynamics 的排斥势函数持续地应用在分子动力学模拟当中. 这里要讨论的是 Metadynamics 最新发展的计算自由能的方法. 假设 Metadynamics 要对 N 维的反应路径求解自由能势能面, 使反应坐标 S 为 N 维的向量:

那么体系的哈密顿量就为

其中 VG(S,t)为额外添加的排斥势能, 具体表达式为:

其中 ω 分别为 σi为第 i 个变量排斥势能的高斯分布的高度和宽度. 最终根据额外添加的一系列排斥势能, 就可以计算相关自由能数据

其中 U(R)为势能函数.Metadynamics 在求解自由能的时候拥有很多优势, 首先, 它通过一系列的外加排斥势能, 一点点地把研究体系从低自由能区域推向高自由能区域, 加速了对研究体系高自由能事件的取样; 其次, 此方法可以设定高纬度的反应路径, 进而求解高纬度的自由能数据, 例如二维的自由能势能面(图 3); 再次, 由于研究体系是从低自由能区域一点点地推高到高自由能区域, 此方法在寻找最优反应路径方面拥有特有的优势; 最后, 与上一点相关的是, 此方法不需要事先对研究体系的反应路径有预判, 而 Umbrella Sampling 和 Jarzynski Equality 则在计算自由能数据之前, 就要设置好反应路径. 鉴于上述优点, Metadynamics 迅速成为一种流行的自由能计算方法, 有很多工作使用此方法计算相关的一维以及二维的自由能数据]. 一开始, Metadynamics 计算自由能数据是由独立于常用的分子动力学模拟软件之外的插件PLUMED执行的. 现在此方法已经集成到了多个计算模拟软件中, 例如: NAMD, ORAC, CP2K, 以及 CPMD等.

各个自由能计算方法的原理以及自身优缺点

Free Energy Perturbation 和 Thermodynamic Integration计算自由能数据所需的热力学循环示意图

(a) Ytreberg 改进 Jarzynski Equality 所用到的限制项Ua(ϕ, θ); (b) 本课题组提出的混合算法得到的自由能数据与 Umbrella Sampling 结果比较吻合

Metadynamics 方法得到的配体-蛋白相互作用体系的二维自由能数据

各个自由能计算方法的原理以及自身优缺点

源文链接：《自由能计算方法及其在生物大分子体系中的适用性问题》