35、生物材料建模与模拟：分子动力学与蒙特卡罗方法

Tomato

于 2025-09-19 16:06:28 发布

阅读量229

点赞数

CC 4.0 BY-SA版权

分类专栏：生物材料前沿创新文章标签：分子动力学蒙特卡罗方法 Verlet算法

本文链接：https://blog.youkuaiyun.com/tomato/article/details/151887905

生物材料前沿创新专栏收录该内容

36 篇文章 ¥499.90

订阅专栏¥69.90

会员秒杀 ¥9.9 重磅福利

超级会员免费看

生物材料建模与模拟：分子动力学与蒙特卡罗方法

1. 分子动力学算法基础

在由 N 个原子组成的系统中，系统的总势能可视为各原子位置的总函数 U(r1 …… rn)。依据经典牛顿力学原理，系统中每个原子所受的力可表示为：
[
F_i = -\nabla U = -\left(\frac{\partial U}{\partial x_i}\vec{i} + \frac{\partial U}{\partial y_i}\vec{j} + \frac{\partial U}{\partial z_i}\vec{k}\right)
]
其加速度为：
[
a_i = \frac{F_i}{m_i}
]
经过时间 t 后，其速度和位置可表示为：
[
\vec{v} i(t) = \vec{v} {i0} + \vec{a} i t
]
[
\vec{r}_i(t) = \vec{r} {i0} + \vec{v}_{i0} t + \frac{1}{2}\vec{a}_i t^2
]

在分子动力学（MD）模拟中，描述分子轨迹需运用有限差分技术求解牛顿运动方程以获取其速度和位置。求解运动方程的基础是假设系统中粒子的位置、速度和加速度可通过泰勒级数展开：
[
\vec{r}(t + \delta t) = \vec{r}(t) + \delta t\vec{v}(t) + \frac{1}{2}\delta t^2\vec{a}(t) + \frac{1}{6}\delta t^3\vec{b}(t) + \frac{1}{24}\delta t^4\vec{c}(t) + \cdots
]
[
\vec{v}(t + \delta t) = \vec{v}(t) + \delta t\vec{a}(t) + \frac{1}{2}\delta t^2\vec{b}(t) + \frac{1}{6}\delta t^3\vec{c}(t) + \cdots
]
[
\vec{a}(t + \delta t) = \vec{a}(t) + \delta t\vec{b}(t) + \frac{1}{2}\delta t^2\vec{c}(t) + \cdots
]

1.1 常用算法

MD 模拟中最常用的方法有 Verlet 算法、基于基本 Verlet 算法衍生的 Velocity - Verlet 算法和蛙跳法。
- Verlet 算法 ：核心思想是通过 t 时刻和 t - δt 时刻粒子的位置预测 t + δt 时刻粒子的位置。用 t 时刻的位置分别表示 t + δt 时刻和 t - δt 时刻的位置，其泰勒展开式分别为：
[
\vec{r}(t + \delta t) = \vec{r}(t) + \delta t\vec{v}(t) + \frac{1}{2}\delta t^2\vec{a}(t)
]
[
\vec{r}(t - \delta t) = \vec{r}(t) - \delta t\vec{v}(t) + \frac{1}{2}\delta t^2\vec{a}(t)
]
两式相加并忽略高阶项，t + δt 时刻的位置为：
[
\vec{r}(t + \delta t) = 2\vec{r}(t) - \vec{r}(t - \delta t) + \delta t^2\vec{a}(t)
]
该算法每次积分仅计算一次力，简单方便，但缺点是不给出速度。在恒温系统中，可通过 t + δt 时刻和 t - δt 时刻粒子的位置计算粒子的速度：
[
\vec{v}(t) = \frac{\vec{r}(t + \delta t) - \vec{r}(t - \delta t)}{2\delta t}
]
- Velocity - Verlet 算法 ：基于 Verlet 算法，不仅提供原子的位置，还增加了速度和加速度的计算，同时保证计算精度。
[
\vec{r}(t + \delta t) = \vec{r}(t) + \delta t\vec{v}(t) + \frac{1}{2}\delta t^2\vec{a}(t)
]
[
\vec{v}(t + \delta t) = \vec{v}(t) + \frac{1}{2}\delta t[\vec{a}(t) + \vec{a}(t + \delta t)]
]
- 蛙跳法 ：引入某一时刻前后 0.5δt 的速度来求解位置：
[
\vec{r}(t + \delta t) = \vec{r}(t) + \delta t\vec{v}(t + \frac{1}{2}\delta t)
]
[
\vec{v}(t + \frac{1}{2}\delta t) = \vec{v}(t - \frac{1}{2}\delta t) + \delta t\vec{a}(t)
]
t 时刻速度的解可表示为：
[
\vec{v}(t) = \frac{1}{2}\left[\vec{v}(t - \frac{1}{2}\delta t) + \vec{v}(t + \frac{1}{2}\delta t)\right]
]
蛙跳法的优点是提高了计算精度并提供速度和加速度信息，但缺点是速度和位置的求解不同步，导致无法确定某一位置处速度对系统能量的贡献。

1.2 算法对比

算法名称	优点	缺点
Verlet 算法	每次积分仅计算一次力，简单方便	不给出速度
Velocity - Verlet 算法	提供位置、速度和加速度计算，精度高	-
蛙跳法	提高计算精度，提供速度和加速度信息	速度和位置求解不同步

1.3 算法选择流程

graph TD;
    A[开始] --> B{是否需要速度信息};
    B -- 否 --> C[选择 Verlet 算法];
    B -- 是 --> D{是否要求速度和位置同步求解};
    D -- 是 --> E[选择 Velocity - Verlet 算法];
    D -- 否 --> F[选择蛙跳法];
    C --> G[结束];
    E --> G;
    F --> G;

2. MD 模拟的边界条件与系综

2.1 周期性边界条件

在现有计算机资源下，MD 模拟可处理的粒子数量有限，研究的系统可视为宏观系统的缩小版，模拟盒会产生与实际实验不匹配的边界效应。为此采用周期性边界条件，模拟盒被其平移镜像包围。当计算系统中的任何粒子离开盒子时，其镜像中对应的粒子必须从相反方向移入盒子，从而保持系统中粒子数量和密度恒定。

在计算分子间相互作用力时，可使用最近镜像法。在非键长程力场计算中，截断半径法可解决重复计算问题。若分子间距离大于设定的截断半径，相互作用可计为零，最大截断半径应小于盒子宽度的一半。由于截断半径处的势能函数不为零，若完全忽略会导致能量曲线不连续，通常将开关函数乘以原始势能函数，使函数曲线从开关起点到截断半径逐渐调整为零。

2.2 MD 模拟的系综

系综是系统在相同宏观条件下许多可能状态的集合。根据设计需求，MD 需要在不同系综中实现。常用的系综根据不同约束条件可分为微正则系综、正则系综、等温等压系综等。
- 微正则系综 ：是孤立保守系统的统计系综，主要是原子数（N）、体积（V）和能量（E）保持恒定，也称为 NVE 系综。由于能量恒定，系统在相空间沿特定能量轨迹运动，模拟开始时的初始配置和条件难以满足既定能量要求，通常的做法是用较长的平衡时间调整合适的初始结构以实现既定能量条件。
- 正则系综 ：与微正则系统不同，正则系综中温度（T）保持不变，总动量为零，也称为 NVT 系综。在恒温条件下，系统可与外部环境交换能量，难以保持温度恒定，通常的解决方案是连接外部热浴以维持系统热平衡。常用的热浴方法有高斯热浴、Nose - Hoover 热浴和 Berendsen 热浴。
- 等温等压系综 ：是实验中常用的系综，与常规系综不同，系统的压力（P）保持不变，也称为 NPT 系综。温度调节可通过 Berendsen 热浴等方法控制，压力调节通常通过控制盒子体积实现，常用方法有 Anderson 方法、Berendsen 方法和 Parrinello - Rahman 方法。

2.3 系综选择列表

系综名称	特点	适用场景
微正则系综（NVE）	原子数、体积和能量恒定	研究孤立保守系统
正则系综（NVT）	原子数、体积和温度恒定	需要恒温条件的系统
等温等压系综（NPT）	原子数、压力和温度恒定	实验中常用，需控制压力和温度

3. 增强采样方法

MD 模拟虽已成为重要研究方法，但采样不足常限制其应用。这是因为系统的能量景观粗糙，许多局部极小值被高能垒分隔，导致准遍历问题，使采样难以跨越整个配置空间，影响蛋白质动力学表征。在生物系统中，大的构象变化对蛋白质活性至关重要，为此有以下几种增强采样方法。

3.1 并行回火

并行回火，也叫复本交换，能有效提高空间采样的代表性。该算法最初由 Swendsen 和 Wang 在 1986 年为 MC 模拟提出，即并行回火 MC，后 Sugita 和 Okamoto 将其扩展到 MD 模拟，即复本交换分子动力学（REMD）。其基本思想是同时模拟原始系统的 M 个复本，每个复本处于不同温度，模拟过程中各复本通过 Metropolis–Hastings 准则与相邻复本交换构象。高温可广泛采样相空间，低温能精确采样局部区域，REMD 常用于研究聚合物构象、肽/蛋白质折叠机制和晶体结构，以减少达到平衡态的弛豫时间。

3.2 模拟退火

模拟退火算法可看作是用于研究系统状态分布和解决全局优化问题的马尔可夫链 MC 方法的扩展。首先将研究系统加热到高温，使分子或原子能自由移动以找到更合适位置，温度缓慢下降时，分子或原子动能相应降低，运动受限，构象变化变小，最终系统可达到最低全局能量构象。由于初始构象能量高，多次独立模拟应能采样到研究系统所有可能的稳定构象。

目前常用的退火方法有经典模拟退火（CSA）、快速模拟退火（FSA）和广义模拟退火（GSA）。其中，系统较大时，GSA 方法效率优于 CSA 和 FSA。GSA 能更均匀地探索构象空间，在低温下实现构象状态的大幅跳跃，主要用于长程效应的全局优化，如原子参数化、引力系统以及小分子和肽的构象优化，适合表征灵活性较大的系统，常规模拟退火只能模拟一些小蛋白质，而 GSA 可用于模拟大生物大分子复合物。

3.3 伞形采样

伞形采样由 Torrie 和 Valleau 在 1977 年提出，用于解决传统 MD 模拟在高自由能区域的采样问题。当研究的初始和最终状态被高能垒分隔时，传统模拟方法需长时间模拟才能采样到两种状态。为解决高自由能区域的采样问题，可在初始状态和最终状态之间设置反应坐标（ξ），将整个反应路径划分为多个区域（窗口）。

对于每个窗口 i，向研究系统添加仅依赖于反应坐标的偏置势 ωi(ξ)。当添加的偏置势能是与反应坐标相关、强度为 K 的谐波势时，对于给定位置 ξi，有：
[
\omega_i(\xi) = \frac{1}{2}K(\xi - \xi_i)^2
]
将研究系统的位置限制在中心位置 ξi 附近，以增强系统在该窗口的采样。在相应谐波势的限制下，通过 MD 模拟可获得系统在反应路径上的一系列偏置概率分布 Pi(ξ)。相邻窗口的概率分布必须有足够重叠才能得到准确的自由能曲线。最后，使用加权直方图分析方法结合这一系列概率分布，消除附加谐波势的影响，得到研究系统沿反应坐标的准确概率分布，进而求解相关自由能曲线。

伞形采样有严格的理论基础，可获得非常准确的自由能数据，但需要许多平衡轨迹和计算资源进行数据收集，对计算资源和系统大小有较多要求和限制，在生物大分子领域，更多用于测试新提出的自由能计算方法的准确性。还有一种伞形采样的变体，即自适应伞形采样，可及时调整偏置势能，使能量在始末状态间均匀分布。

3.4 元动力学

元动力学基于选定的集体变量，能有效重建研究系统的自由能表面并加速采样。该算法由 Laio 和 Parrinello 在 2002 年提出并应用于 MD 模拟，可形象地描述为“用计算沙填充自由能阱”。算法假设系统可用一些反应坐标描述，模拟中计算系统定期向系统添加一定的高斯排斥势，如下所示：
[
V(S,t) = \sum_{i=1}^{d}\sum_{k:\tau_k < t}W(\tau_k)\exp\left[-\frac{(S_i(q(t)) - S_i(q(\tau_k)))^2}{2\sigma_i^2}\right]
]
其中，q(t) 是系统对应于时间 t 的微观坐标，Si(q) 是基于微观坐标 q 的第 i 个反应坐标，σi 是第 i 个反应坐标的高斯势宽度，d 是反应坐标的数量，τ 是添加高斯势的时间间隔，k 是正整数，W(t) 是时间 t 时高斯势的高度。

随着模拟进行，越来越多的高斯势被添加到轨迹中，直到系统探测到足够的能量景观。当修改后的自由能相对于反应坐标变为常数时，反应坐标开始大幅波动，同时可恢复能量轮廓，它与高斯势之和相反：
[
V(S,t)\to\infty = -F(S) + C
]
通过调整添加高斯势的时间间隔以及高斯势的高度和宽度，可优化准确性与计算成本的比率。

为解决元动力学方法的各种问题和局限性（如收敛性、采样、多反应坐标等），提出了一些新变体，如高斯势逐渐减小的良温元动力学、基于多个复本和相同反应坐标的多步元动力学、基于多个复本和不同反应坐标的偏置交换元动力学、基于不同温度复本和相同反应坐标的并行回火元动力学。元动力学方法主要用于研究蛋白质折叠、化学反应、分子对接和相变，最近也用于模拟肽在表面的吸附以及计算蛋白质 - 配体和蛋白质 - 蛋白质结合的自由能。目前，专门用于元动力学模拟的软件是 PLUMED，该方法已集成到许多常用的分子模拟软件中，如 AMBER、GROMACS、LAMMPS 和 NAMD。

3.5 增强采样方法对比

方法名称	优点	缺点	适用场景
并行回火	有效提高空间采样代表性，减少弛豫时间	-	聚合物构象、肽/蛋白质折叠机制和晶体结构研究
模拟退火	可解决全局优化问题，能采样到稳定构象	计算量大	长程效应全局优化、灵活性大的系统
伞形采样	有严格理论基础，可获准确自由能数据	需大量计算资源，对系统大小有要求	测试新自由能计算方法准确性
元动力学	加速采样，重建自由能表面	存在收敛性等问题	蛋白质折叠、化学反应、分子对接和相变研究

3.6 增强采样方法选择流程

graph TD;
    A[开始] --> B{是否需要解决高自由能区域采样问题};
    B -- 是 --> C{是否关注全局优化};
    C -- 是 --> D[选择模拟退火];
    C -- 否 --> E[选择伞形采样];
    B -- 否 --> F{是否需要提高空间采样代表性};
    F -- 是 --> G[选择并行回火];
    F -- 否 --> H{是否需要重建自由能表面};
    H -- 是 --> I[选择元动力学];
    H -- 否 --> J[选择传统 MD 模拟];
    D --> K[结束];
    E --> K;
    G --> K;
    I --> K;
    J --> K;

4. 蒙特卡罗方法

许多材料系统问题具有概率和统计性质，需要新方法解决。蒙特卡罗（MC）方法基于简单理论准则，采用重复随机采样方法处理复杂系统问题，通过随机采样模拟对象的概率和统计。MC 模拟方法也叫随机模拟方法或统计模拟方法，通过特定“实验”方法可获得事件的频率（或变量的某些特征值），将其视为随机事件的概率（或随机变量的期望值）以解决相关问题。使用 MC 方法模拟过程时，需要生成具有各种概率分布的随机变量，最简单和最重要的随机变量是在 [0,1] 上均匀分布的随机数。

传统 MC 方法简单，但需大量采样才能得到较准确结果，计算量很大。为此，Metropolis 等人提出重要采样方法，即按一定概率选择是否接受新状态。MC 模拟因算法简单、易于实现以及计算机科学的快速发展而广泛应用，可用于金融工程、宏观经济学、计算物理学（如粒子输运问题、量子热力学、空气动力学）、医学、生物学和勘探等领域，如今也用于模拟蛋白质在固体表面的吸附。