计算实际例子_自由能计算的相关笔记

川长思鸟来 2023-01-07 10:08 163阅读 0赞

最近为了准备以后的课题读了一下Tuckerman的Statistical Mechanics: Theory and Molecular Simulation，以下简要整理一下第八章Free energy calculations的笔记。以下每个部分的编号是我自己加的，并非书中的章节编号。

1.  自由能的计算

自由能的计算可以利用公式

![Image 1][] （书中用记号

![Image 1][] 表示自由能，本文会用

![Image 1][] 表示自由能。下文中

![Image 1][] 可能表示其它含义，例如表示某个状态），自由能可以帮助我们判断反应在热力学上是否能够发生等相关性质，我们还可以从

![Image 1][] 和

![Image 1][] 推算出许多有用的物理量。自由能计算的难点在于涉及配分函数

![Image 1][] ，

![Image 1][] 的计算需要对所有粒子的坐标和动量进行全空间积分。

--------------------

1.1 通过势场定义状态

判断反应热力学能否自发主要根据反应的自由能变化（假设为正则系综）。假设我们考虑的体系有两个状态

![Image 1][] 和

![Image 1][] ，我们需要考虑这两个状态之间转化的反应。首先我们要说明一下如何定义状态

![Image 1][] 和状态

![Image 1][] ，直观上来说，

![Image 1][] 个粒子对应了一个坐标空间

![Image 1][] ，其中某一些点对应于状态

![Image 1][] ，有一些对应于状态

![Image 1][] ，还有一些点可能难以定义；例如，对应于一个配位反应，状态

![Image 1][] 对应配体和受体没有配位的状态，而状态

![Image 1][] 对应配体和受体配位的状态，那么坐标空间中有一些位置我们可能难以定义配体和受体是否算发生了配位（例如配体和受体的距离并不是特别近，但是也没有很远，二者之间有部分的相互作用，但是相互作用的强度可能远小于配位效应的作用）。书中给出的定义是

![Image 1][] 和

![Image 1][] 分别对应两个不同的势能函数

![Image 1][] 和

![Image 1][] ，例如之前提到的配位反应的例子，

![Image 1][] 是没有配位的状态，对应

![Image 1][] 里面不包含配体和受体的相互作用项，而

![Image 1][] 是配位状态，对应

![Image 1][] 里面包含配体和受体的相互作用项。

--------------------

1.2 Free energy perturbation formula

我们关注

![Image 1][] 和

![Image 1][] 的自由能之差，计算方法为：

![Image 1][]

![Image 1][] ，其中

![Image 1][] 或

![Image 1][] 。

实际中我们通过分子动力学模拟采样从而计算自由能，我们可以采取一些操作从而使得

![Image 1][] 的表达式变成某种类似于采样后计算某一物理量均值的形式。例如，我们知道假设我们希望计算状态

![Image 1][] 下某个物理量

![Image 1][] 的均值，则我们可以写出：

![Image 1][]

我们可以用类似的方法处理

![Image 1][] ：

![Image 1][]

**这告诉我们计算**

**![Image 1][] 可以通过首先设定势函数为**

**![Image 1][] 并且进行分子动力学模拟采样，然后对于采样到的构型计算物理量**

**![Image 1][] 的平均值。**

**公式**

**![Image 1][] 被称为free energy perturbation formula，这个公式本身是严格的，但是实际应用的时候由于采样数量有限，其应用会受到限制（见下文）。**

这里我们可能会思考，

![Image 1][] 和

![Image 1][] 应该分别适用于描述坐标空间中的不同的点或者区域，例如还是之前提到的配位的例子，

![Image 1][] 适合于描述配体和受体距离很远而没有发生配位的状态，实际上由于

![Image 1][] 中不包含配位效应的吸引势能项，在

![Image 1][] 下配体和受体反而可能因为位阻而倾向于排斥；但

![Image 1][] 适合于描述配体和受体距离很近且完全配位的状态，

![Image 1][] 中也包含配体和受体的吸引势能项，二者倾向于靠近并且配位。

实际上书中提到，该方法只适用于状态

![Image 1][] 和状态

![Image 1][] 差异较小，或者说二者实际上对应于坐标空间中的构型有较多重叠。我举例说明一下我的理解，还是以配位作为例子：我们假设配体和受体由于有空间位阻，二者之间有排斥效应，同时还有配位效应，并且配位效应相比排斥效应更加显著，这就导致非配位状态

![Image 1][] 和配位状态

![Image 1][] 对应于坐标空间中差异很大的区域（

![Image 1][] 可能对应二者距离很远的构型，而

![Image 1][] 对应于二者距离接近且空间取向合适、发生配位的构型）。如果我们用

![Image 1][] ，即不包含二者配位吸引势能的势场进行分子动力学模拟，那么我们采样到的构型几乎全都是二者相距较远的构型，在有限次数的采样下或许根本采样不到配体和受体配位的构型。但我们要计算的物理量是

![Image 1][] ，由于我们这里假设配位效应相比于位阻排斥效应更强，所以当配体和受体靠近并且呈现配位的构型的时候

![Image 1][] 很小，从而导致

![Image 1][]很大，也就是说配体和受体呈现配位构型对于计算

![Image 1][] 的均值的贡献很大，这些构型对于计算

![Image 1][]非常重要，哪怕我们在

![Image 1][] 势场下进行采样会导致配位构型出现的概率很小。

**总结一下这个例子：在**

**![Image 1][] 势场下进行采样，位阻排斥效应导致配体和受体发生配位的构型出现的几率很小，但是对于**

**![Image 1][]的均值的计算这些小概率的构型却非常重要。**

但是，由于之前提到了，有限次数的采样很可能完全采样不到二者配位的构型，从而导致基于采样计算出的

![Image 1][]的均值相比于真实的均值偏差很大。这是采样数量终究是有限的所导致的问题，如果采样数量趋于无穷，那么这个问题也就趋于被解决了，因为free energy perturbation formula本身是严格的。通过这个例子大家可以看出为什么书中说这个方法适用于状态

![Image 1][] 和

![Image 1][] 对应的构型重叠较多；另外，这个例子也可以帮助大家理解在分子动力学模拟中采样次数不够多可能会导致的一些问题。

为了解决当状态

![Image 1][] 和

![Image 1][] 重叠较少的时候无法直接使用free energy perturbation formula的问题，我们引入一系列中间态，从而使得体系的势能可以缓慢地从

![Image 1][] 变化到

![Image 1][] 。令

![Image 1][] ，引入一系列中间态，对应于势能

![Image 1][] ，从而构成一个势能序列

![Image 1][] ，序列中相邻两项的势能差异很小，并且相邻两项对应的态重叠较多，从而可以对相邻两项使用free energy perturbation formula。利用自由能是状态函数的性质，我们可以得到：

![Image 1][]

上述方法可以看作通过一系列离散的台阶连接状态

![Image 1][] 和

![Image 1][] ，我们要确保相邻两个台阶足够接近。

--------------------

1.3 metapotential & thermodynamic integration formula。

注：以下为了方便起见，用

![Image 1][] 表示

![Image 1][]

接下来介绍一种连续地、绝热地从状态

![Image 1][] 变化到状态

![Image 1][] 的方法，即一条绝热路径，书上对于绝热路径的定义是：

**An adiabatic path is one along which the system is fully relaxed at each point of the path\[1\].** 引入参数

![Image 1][] 用于描述这条连续的绝热路径，从而可以定义metapotential

![Image 1][] 。

![Image 1][]

这里

![Image 1][] ，因此

![Image 1][] ，由此可见 参数

![Image 1][] 就像一个旋钮，通过缓慢打开这个旋钮可以使得势场连续地从

![Image 1][] 变为

![Image 1][] 。例如，一种显然的选择是定义

![Image 1][] 。

如何利用metapotential呢？我们记势场为

![Image 1][] 对应的配分函数为

![Image 1][] ，其表达式为：

![Image 1][]

其中

![Image 1][] 为一个系数。由此可得：

![Image 1][]

**这里在尖括号右下角的**

**![Image 1][] 表示在固定**

**![Image 1][] 的取值的情况下（即固定metapotential的条件下）进行模拟采样后取平均。**

可得：

![Image 1][]

公式

![Image 1][] 被称为thermodynamic integration formula。

当

![Image 1][] 的时候，

![Image 1][]

--------------------

1.4 例子(1)

书中给出了一个例子\[1\]：我们希望求出酶E和抑制剂I结合生成配合物EI的反应的自由能变化，thermodynamic integration formula告诉我们可以设计连续的path使得势场从E和I分离（即不考虑配位效应）转化为E和I配位（即考虑配位效应）。实际的反应是发生在溶液中的，初始E和I分离的时候二者分别被溶剂化并且相距较远，我们可能难以采样到二者配位的构型；一种变通的方法是设计如下的path：首先E和I分别去溶剂化，其次二者在真空中发生配位反应生成EI，最后EI再被溶剂化，由于自由能只是状态函数，因此这个path是合理的。这种方法被称为double decoupling method。

thermodynamic integration formula这一方法的问题在于，为了计算积分

![Image 1][] 我们需要选取许多的

![Image 1][] 的数值并且每一个取值都要进行模拟采样，因此计算量较大。另外，

![Image 1][] 和

![Image 1][] 对应实际中出现的状态，但是

![Image 1][] 中的所有值对应的状态都是我们假想的。

--------------------

1.5 Free energy profile

为了缓解使用thermodynamic integration formula的时候需要采样许多

![Image 1][] 的虚拟态的问题，我们可以通过一些trick：调整

![Image 1][] 和

![Image 1][] ，从而使得我们的采样可以集中在

![Image 1][] 和

![Image 1][] 附近。另外，以下我们将

![Image 1][] 看作一个额外的自由度，即现在相空间中还需要考虑广义坐标

![Image 1][] 以及与之对应的共轭的广义动量

![Image 1][] ，还有一个类似质量的参数

![Image 1][] （实际上这里

![Image 1][] 和

![Image 1][] 的引入只是为了使得哈密顿量的形式显得严格，对于之后的处理没有太大的影响，可以不必纠结

![Image 1][] 和

![Image 1][] ）。这一部分，即书中的8.3节，写得比较迷，个人感觉有些地方有些问题，这里我就按照自己的理解整理了。

仍然采用之前定义的metapotential

![Image 1][] ，并且定义

![Image 1][] 和

![Image 1][] ：

![Image 1][]

我们可以定义一个概率函数

![Image 1][] （关于

![Image 1][] 的定义，书上给出的是

![Image 1][] ，个人认为这里右下角不应该写

![Image 1][] ，即平均值不应该在固定

![Image 1][] 取值的情况下计算）：

![Image 1][]

**![Image 1][] 具有某种类似概率密度函数的性质，有**

**![Image 1][] ，其物理意义就是如果我们将**

**![Image 1][] 也看作一个额外的自由度，从而拓展了相空间（拓展后的相空间加入了**

**![Image 1][] 和**

**![Image 1][]）。**

在拓展后的向空间中取到满足

![Image 1][] 的态的概率密度为

![Image 1][] 。另外由于我们之后关注的是

![Image 1][] 的表达式，其中分子和分母中都包含了对于

![Image 1][] 的积分项，这些项会相互抵消，因此之前说了

![Image 1][] 的表达式中的

![Image 1][] 项并不重要。

**定义free energy profile**

**![Image 1][] ，由此我们将自由能和概率密度联系了起来。**

如前文所说，我们希望在

![Image 1][] 和

![Image 1][] 附近的采样更密集，而在这两点之间区域的采样更少，这就意味着我们希望

![Image 1][] 在

![Image 1][] 和

![Image 1][] 附近的概率密度较大，而在这两点之间的区域概率密度较小（这一点仔细思考似乎倒不是显然的，感觉自己也没有充分理解这里

![Image 1][] 到底意味着什么）；同时我们也发现

![Image 1][] 越大，

![Image 1][] 越小。我们的目标就是选取合适的

![Image 1][] 和

![Image 1][] ，使得

![Image 1][] 在

![Image 1][] 和

![Image 1][] 附近较小（

![Image 1][] 在

![Image 1][] 和

![Image 1][] 附近较大），

![Image 1][] 在

![Image 1][] 和

![Image 1][] 之间的区域较大（

![Image 1][] 在

![Image 1][] 和

![Image 1][] 之间的区域较小），

**也就是说理想的**

**![Image 1][] 函数在**

**![Image 1][] 从**

**![Image 1][] 变化到**

**![Image 1][] 的过程中先上升后下降，即存在一个势垒。**

--------------------

1.6 例子(2)

以下给出一个书上的例子\[1\]，计算过程见文末附录。

考虑两个振动方向正交的谐振子，

![Image 1][] 态表示二者没有耦合，

![Image 1][] 态表示二者完全耦合。

![Image 1][] 和

![Image 1][] 的表达式如下：

![Image 1][]

经过繁琐的计算后（计算过程见附录1）得到：

![Image 1][]

从而可以写出

![Image 1][] 的表达式，这里我们略去了常数项：

![Image 1][]

下图显示了两种选择函数

![Image 1][] 和

![Image 1][] 的方式（函数表达式见图片中的说明）得到的

![Image 1][]和

![Image 1][] 的图像（图片来自于Tuckerman的Statistical Mechanics: Theory and Molecular Simulation）。可以看出相比之下

![Image 1][] 是一种更好的选择，因为其free energy profile在

![Image 1][] 和

![Image 1][] 之间有个明显的势垒。

![c644b918f4b621386ffcdb0fc37f197b.png][]

另外，

![Image 1][] 也可以写作

![Image 1][] ，其中

![Image 1][]

书中8.3节还有一部分讲了adiabatic free energy dynamics (AFED)，暂时还没太搞懂这个方法，以后有空再来填坑吧。

书中8.4节先跳过。

--------------------

1.7 Collective variables (CV)

我们可以通过选取合适的CV描述反应的过程。例如一个简单的双原子反应，我们可以选取两个原子中心的距离作为CV，从而方便地描述这个反应。一般来说CV是体系中所有原子的坐标的函数。

我们可以用Radius of gyration

![Image 1][] 作为长链分子的CV，其定义为：

![Image 1][]

其中

![Image 1][] 是主链上的原子数。 直观上，

![Image 1][] 描述了主链总体的伸展程度。

选定一组CV

![Image 1][] ，其中

![Image 1][] ，我们可以写出这组CV取值为

![Image 1][] 的时候对应的概率密度

![Image 1][] 以及相应的自由能

![Image 1][] ：

![Image 1][]

--------------------

1.8 blue moon ensemble approach

接下来介绍一种计算方法，该方法可以计算沿着某个CV方向有势垒情况下该方向的free energy profile，该方法被称为blue moon ensemble approach。假设只有一个CV：

![Image 1][] ，当

![Image 1][] 的时候对应的概率密度

![Image 1][] 和自由能的计算公式参见本文1.7节。我们想要采样的路径对应的CV是从初态

![Image 1][] 到末态

![Image 1][] ，为此我们引入序列

![Image 1][] ，该序列对应了从初态到末态的一条路径。每次我们取序列中的一个点

![Image 1][] ，固定在

![Image 1][] （即令

![Image 1][] 满足

![Image 1][] 条件）。另外，blue moon ensemble approach并不会直接得到

![Image 1][] （不太清楚为啥不是直接得到

![Image 1][] ），而是得到

![Image 1][] 对

![Image 1][] 的导数：

![Image 1][]

从而我们通过积分进一步求出自由能变化：

![Image 1][]

考虑

![Image 1][] 对

![Image 1][] 的导数：

![Image 1][]

做坐标变换

![Image 1][] ，并且

![Image 1][] ，通过

![Image 1][] 可以得到与之共轭的

![Image 1][] ，并且可以保证

![Image 1][] ，即坐标变换对应的雅克比行列式的绝对值为1（暂时不太清楚如何确保可以找到这样一组坐标变换，使得坐标变换对应的雅克比行列式的绝对值为1，根据书中下文的意思，这里是因为我们的坐标变换中

![Image 1][] 和

![Image 1][] 都参与了，从而保证了坐标变换对应的雅克比行列式绝对值为1；书上提到了，如果只有

![Image 1][] 参与坐标变换， 而

![Image 1][] 不变，那么就无法保证雅克比行列式的绝对值为1。这里具体的数学细节需要进一步深究）。

![Image 1][]

其中，上式的倒数第三行用到了分部积分（分部积分后产生的另一项为何为0我也有点迷...），最后一行尖括号右下角的

![Image 1][] 说明此平均是在

![Image 1][] 的条件下取到的，或者说这个平均对应的权重因子为

![Image 1][] 。

![Image 1][]

上述坐标变换涉及了从

![Image 1][] 变换到

![Image 1][] 和

![Image 1][] ，而另一种坐标变换的方式只涉及从

![Image 1][] 变换到

![Image 1][] ，即首先将动量部分积分后再将

![Image 1][] 变换为

![Image 1][] ：

![Image 1][]

其中

![Image 1][] 是雅克比行列式，

![Image 1][] ，

![Image 1][] ，右下角的角标

![Image 1][] 表示这个平均对应的权重因子为

![Image 1][] （由于这里的坐标变换只涉及

![Image 1][] ，没有涉及

![Image 1][] ，所以我们也可以说这个平均对应的权重因子为

![Image 1][]，因为和动量部分相关的项的积分都被消去了）。

这里有一个小疑问，我记得坐标变换的时候雅克比行列式是需要取绝对值的，这里

![Image 1][] 是正的（否则不能对

![Image 1][] 取对数），那么如何确保

![Image 1][] 是正的呢？我的理解，如果在某个局部

![Image 1][] 小于

![Image 1][] ，可以通过取某一个

![Image 1][] 为

![Image 1][] ；但是会不会某些情况下整个函数的某些局部为正，某些局部为负呢...？欢迎大家讨论、指点（还是数学不够好）。

--------------------

1.9

--------------------

附录1：推导正交的谐振子耦合的free energy profile

![Image 1][]

![Image 1][] 这里积分号里面的项的指数是一个二次型，我们进行如下的换元：

![Image 1][]

二次型

![Image 1][] 对应

![Image 1][] ，其中矩阵

![Image 1][] ，

![Image 1][] ，将

![Image 1][] 进行对角化，可得

![Image 1][] ，其中

![Image 1][] 是正交矩阵，

![Image 1][] 是对角矩阵。经过计算可得：

![Image 1][]

其中

![Image 1][]

坐标变换

![Image 1][]

根据雅克比行列式可以计算出

![Image 1][] （实际上这里的雅克比行列式就是

![Image 1][] ）

所以

![Image 1][]

（吐槽一下，这个计算量也太大了吧...）

--------------------

附录2 含有限制的拉格朗日方程和哈密顿方程

(1) 含有限制的拉格朗日方程

如果对于相空间中自由度的限制有

![Image 1][] 个，则剩余的自由度为

![Image 1][] ，此时体系至少需要

![Image 1][] 个CV进行描述。限制分为两类，一类是holonomic constraint，指限制的方程只包含广义坐标和时间，即

![Image 1][] ；其余种类的限制都属于non-holonomic constraint，例如

![Image 1][] 。

对于有限制的体系，由于坐标之间彼此关联，因此不可以直接变分，需要用Lagrange undetermined multipliers将限制条件引入我们所需要处理的极值问题（作用量极小化，根据最小作用量原理）。假设我们的体系受到的限制满足

![Image 1][] （所有的holonomic constraint写成全微分的形式即可转化为这样的表达式），我们首先写出不受限制的体系对应的作用量

![Image 1][] ，然后用拉格朗日乘子将限制条件包含进去，最后转化为变分求极值。

如果体系不受限制，则：

![Image 1][]

现在我们将限制条件加入，并且引入拉格朗日乘子

![Image 1][] ，可得：

![Image 1][]

![Image 1][] （对应于

![Image 1][] 个方程）

因为我们对于时间不去做变分，所以引入的拉格朗日乘子对应的项只有

![Image 1][] 。

限制条件对应于

![Image 1][] 个方程

![Image 1][] ，因此我们共有

![Image 1][] 个方程和

![Image 1][] 个未知数

![Image 1][] 。

(2) 含有限制的哈密顿方程

接下来推导一下对应的含有限制的哈密顿方程，仿照朗道力学里面推导哈密顿方程的方法。定义广义动量

![Image 1][] ，再定义哈密顿量

![Image 1][] ，可得：

![Image 1][]

由此可得含有限制条件的哈密顿方程：

![Image 1][]

看一下在限制条件

![Image 1][] 下系统的哈密顿量是否守恒：

![Image 1][]

其中最后一步来自于限制条件。由此可见**如果限制条件不显含时间**

**![Image 1][] 的话**

**![Image 1][] ，从而哈密顿量守恒。**

--------------------

附录3 计算含有限制的运动问题的例子

我们以一个简单的例子说明如何具体计算附录2给出的含有限制的运动问题，主要参考了书中1.10节。书中的1.10节标题为Gauss's principle of least constraint，不过仔细读了读，作者只是给出了一些公式，但似乎并没有给出Gauss's principle of least constraint的直观解释。

考虑一个粒子，其坐标为

![Image 1][] ，其受到限制

![Image 1][] 。根据附录2，其运动方程可以写作：

![Image 1][]

在附录2中，我们定义了广义动量

![Image 1][] ，我们还可以定义广义力

![Image 1][] ，当我们将广义坐标

![Image 1][] 定义为笛卡尔坐标

![Image 1][] 的时候广义动量和广义力就会变成我们一般说的动量和外力。

公式

![Image 1][] 告诉我们，粒子的运动既收到来自外场提供的主动的外力

![Image 1][] ，同时由于粒子被限制在平面