比较分子动力学模拟中Verlet算法、Speed Verlet算法及辛算法的精度

比较分子动力学模拟中Verlet算法、Speed Verlet算法及辛算法的精度

原理

分子动力学模拟是一种研究凝聚态系统的有力方法，对于经典体系，分子动力学假定所有粒子的运动服从牛顿运动方程，只要给出粒子的初始位置、速度（初始条件），边界条件，以及粒子之间相互作用的势函数，就能得到系统在未来任意时刻的状态信息。
由于多原子体系的相互作用非常复杂，难以用解析法求解体系的运动方程，因此用数值方法求解运动方程的近似解。这个过程就不可避免地引入了误差，包括截断误差和舍入误差。前者与不同算法有关，如有限差分法采用$Taylor$展开后截断某些高阶项会引入截断误差，后者来源于计算机本身的精度。两种误差都能通过减小时间步长$\delta t$来减小，当$\delta t$较大时，截断误差起主要作用，但随着$\delta t$的减小会迅速降低，舍入误差随$\delta t$变化不大，并在$\delta t$较小时起主要作用。
进行数值求解的算法有很多种，对于差分法，其共同点是将粒子的位置，速度，加速度展开为$Taylor$级数。下面我们来分析常用的$Verlet$算法和更精确一些的$Speed$ $Verlet$算法的具体步骤以及它们的截断误差。

$Verlet$ $Algorithm$

$Verlet$算法在分子动力学模拟中应用极为广泛，而且原理也很简明。这种方法利用粒子在$t$时刻的位置和加速度，以及$t-\delta t$时刻的位置，计算出$t+\delta t$时刻的位置。
将粒子的位置用$Taylor$公式展开，忽略四阶以上小量：
$$r(t+\delta t)=r(t)+\frac{d}{dt}r(t)\delta t+\frac{1}{2!}\frac{d^2}{d{t}^2}r(t)(\delta t{)}^2+\frac{1}{3!}\frac{d^3}{d{t}^3}r(t)(\delta t{)}^3+o[(\delta t{)}^4]$$
将式中的$\delta t$换为$-\delta t$得：
$$r(t-\delta t)=r(t)-\frac{d}{dt}r(t)\delta t+\frac{1}{2!}\frac{d^2}{d{t}^2}r(t)(\delta t{)}^2-\frac{1}{3!}\frac{d^3}{d{t}^3}r(t)(\delta t{)}^3+o[(\delta t{)}^4]$$
以上两式相加，并忽略$o[(\delta t{)}^4]$小量得：
$$r(t+\delta t)=2r(t)-r(t-\delta t)+\frac{d^2}{d{t}^2}r(t)(\delta t{)}^2$$
两式相减，并忽略$o[(\delta t{)}^2]$小量得：
$$v(t)=\frac{dr}{dt}=\frac{1}{2\delta t}[r(t+\delta t)-r(t-\delta t)]$$
至此，我们得到了$t+\delta t$时刻粒子位置和$t$时刻速度的表达式，并且知道了位置的误差正比于$o[(\delta t{)}^4]$，速度的误差正比于$o[(\delta t{)}^2]$。$Verlet$算法的数值稳定性比简单的欧拉方法高很多，并保持了物理系统中的时间可逆性与相空间体积元体积守恒的性质。然而，这种算法得出的轨迹与速度无关（无法与热浴耦联），并且不能给出同一时刻的位置和速度。

$Speed$ $Verlet$ $Algorithm$

在$Verlet$ $Algorithm$中，我们发现其无法给出同一时刻的速度和位置，$Speed$ $Verlet$（有些地方也叫$Velocity$ $Verlet$）正是用来克服这个缺陷的。
将位置展开到二阶小量：
$$r(t+\delta t)=r(t)+\frac{d}{dt}r(t)\delta t+\frac{1}{2!}\frac{d^2}{d{t}^2}r(t)(\delta t{)}^2$$
引入$t+\frac{\delta t}{2}$时刻的速度：
$$v(t+\frac{\delta t}{2})=v(t)+\frac{1}{2}\frac{d}{dt}v(t)\delta t$$
则$t$