分子模拟/分子动力学

氩的分子动力学模拟

分子动力学模拟是一种模拟多体系统原子运动的工具。分子动力学用于计算经典系统的平衡和传输性质（粘度、热导率、扩散、反应速率、蛋白质折叠时间、结构和表面涂层）。要模拟分子的动力学，可以使用经典的牛顿力学。在实践中，分子动力学模拟是通过以时间的小增量移动来执行的。它是对原子和分子的经典运动方程的求解，以获得系统的演化。由于分子系统通常由大量粒子组成，因此不可能通过解析方法确定这些复杂系统的性质，因此 MD 被应用于这些系统，通过使用数值方法来解决这个问题。

Verlet 算法

Verlet 算法是一种流行的算法，可以预测原子在未来某个时间点的位置。^[1] 该算法使用泰勒级数逼近来近似时间 $t+\Delta t$ （从时间 t 开始的未来一个时间步长）的粒子位置。

泰勒级数基于函数在另一个点 a 的值和导数来近似点 x 处的函数。

$f(x)=f(a)+{\frac {{\textrm {d}}f(a)}{{\textrm {d}}t}}\left(x-a\right)+{\frac {1}{2!}}{\frac {{\textrm {d}}^{2}f(a)}{{\textrm {d}}t^{2}}}\left(x-a\right)^{2}+{\frac {1}{3!}}{\frac {{\textrm {d}}^{3}f(a)}{{\textrm {d}}t^{3}}}\left(x-a\right)^{3}+\ldots$

Verlet 算法的推导

泰勒级数逼近了分子在 $t+\Delta t$ 的位置方程。

$r_{i}(t+\Delta t)\approx r_{i}(t)+{\frac {{\textrm {d}}r_{i}(t)}{dt}}(t+\Delta t-t)+{\frac {1}{2}}{\frac {{\textrm {d}}^{2}r_{i}(t)}{{\textrm {d}}t^{2}}}(t+\Delta t-t)^{2}$

$r_{i}(t+\Delta t)\approx r_{i}(t)+{\frac {{\textrm {d}}r_{i}(t)}{{\textrm {d}}t}}\Delta t+{\frac {1}{2}}{\frac {{\textrm {d}}^{2}r_{i}(t)}{{\textrm {d}}t^{2}}}\Delta t^{2}$

位置的一阶和二阶导数分别对应速度和加速度，加速度通过牛顿第二定律与力相关。

${\frac {{\textrm {d}}r_{i}(t)}{{\textrm {d}}t}}=v_{i}(t)$

${\frac {{\textrm {d}}^{2}r_{i}(t)}{{\textrm {d}}t^{2}}}=a_{i}(t)={\frac {F_{i}(t)}{m_{i}}}$

$r_{i}(t+\Delta t)\approx r_{i}(t)+v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

考虑时间上向后推移一步

$r_{i}(t-\Delta t)\approx r_{i}(t)+{\frac {dr_{i}(t)}{dt}}(t-\Delta t-t)+{\frac {1}{2}}{\frac {d^{2}r_{i}(t)}{dt^{2}}}(t-\Delta t-t)^{2}$

$r_{i}(t-\Delta t)\approx r_{i}(t)-v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

将前一步和下一步的位置相加

$r_{i}(t+\Delta t)+r_{i}(t-\Delta t)=r_{i}(t)+v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}+r_{i}(t)-v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

$r_{i}(t+\Delta t)+r_{i}(t-\Delta t)=2r_{i}(t)+{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

重新排列

Verlet 算法

$r_{i}(t+\Delta t)=2r_{i}(t)-r_{i}(t-\Delta t)+{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

速度 Verlet 算法推导

Verlet 方程的另一种形式，它涉及速度。通过使用当前时间步长的下一个时间步长的泰勒级数

$r_{i}(t+\Delta t)\approx r_{i}(t)+v_{i}(t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

使用下一个时间步长的当前时间步长

$r_{i}(t)\approx r_{i}(t+\Delta t)-v_{i}(t+\Delta t)\Delta t+{\frac {1}{2}}{\frac {F_{i}(t+\Delta t)}{m_{i}}}\Delta t^{2}$

通过替换

速度 Verlet 算法

$v_{i}(t+\Delta t)=v_{i}(t)+{\frac {\Delta t}{2m_{i}}}\left[F_{i}(t)+F_{i}(t+\Delta t)\right]$

另请参阅

维基百科：Verlet 积分 Chemwiki: Lennard-Jones

时间步长

大多数分子动力学模拟的时间步长在飞秒量级，这是化学键振动的量级。分子动力学模拟受最高频率振动的限制，通常时间步长应低于最高频率的十分之一。通常，有机分子分子动力学模拟使用 1-2 fs ( ${10^{-15}}$ s) 的时间步长，因为含有氢的键的振动速度很快。

${\text{Number of time steps}}={\frac {\text{length of simulation}}{\text{time step length}}}$

第一个问题是，一些过程，如蛋白质折叠，需要相对较长的模拟时间才能完整地观察到。为了在 MD 模拟中观察这些过程，我们需要至少微秒级的模拟。

${\text{Number of time steps}}={\frac {10^{-6}}{10^{-15}}}=1{\text{ billion time steps}}$

即使对于计算机来说，这对于时间步长来说也是一个非常大的数字。为了解决这个问题，应用了一些约束条件。许多分子动力学模拟将所有含氢键约束为平衡值。因为氢是一种非常轻的原子，并且具有非常快的键振动。通过将这些含氢键约束为固定在平衡值，快速振动将被忽略。

另一个问题是共价键振动的振动能级在量子能级之间有很大的空间，因此使用经典力学对其描述很差，最好假设键是刚性的。

时间可逆性

一个通过牛顿运动定律传播的孤立系统应该是时间可逆的。也就是说，如果你。在许多情况下，分子动力学算法是时间可逆的。这意味着通过向前运行轨迹，然后向后运行相同数量的步长，可以按相同的步长序列返回起点。因此，两个具有相同初始条件和时间步长的 MD 模拟将遵循完全相同的轨迹。

第二个结果是该系统被称为**确定性**。如果你知道粒子在给定时间点的位移和速度，并且可以计算作用在它们身上的任何位置的力，原则上，你可以积分运动定律以确定系统在任何未来时间点的状态。通过这种方式，我们可以根据我们现在所处的状态来*确定*系统的未来状态（即位置和速度）。

使用分子动力学采样热力学系综

微正则系综 (NVE)

使用近似牛顿动力学的积分器进行的分子动力学模拟（例如 Verlet 算法）将守恒系统的总能量（即动能和势能之和在每个时间步长保持恒定）。因此，这些模拟将对微正则系综 (NVE) 进行采样，因为系统的总能量 (E) 是守恒的。

Verlet 方程以及所有其他积分器都与泰勒级数展开中高阶项的截断相关联的一些误差。当使用更大的时间步长时，此误差会更显着。此外，由于在有限精度计算机（例如 32 或 64 位寄存器）上执行这些计算而导致的舍入误差会产生数值误差。因此，分子动力学模拟的总能量会在模拟过程中偏离其初始值。

采样不同的系综

要采样不同的系综，必须修改运动方程。这通常被描述为使用恒温器或恒压器进行分子动力学模拟。

正则系综 (NVT)

在正则系综中，物质的量 (N)、体积 (V) 和温度 (T) 是恒定的。在这个系综中，吸热和放热过程的能量与恒温器交换。已经开发了几种方法来修改分子动力学模拟的运动方程，以便它们对正则系综进行采样，例如 wikipedia::Andersen 恒温器和 wikipedia::Nosé–Hoover 恒温器。