三、牛顿运动方程与积分算法

1. 牛顿运动方程
对于每个原子i：
$$F_i = m_i a_i = -\frac{\partial U}{\partial r_i}$$
其中F_i为受力，m_i为质量，a_i为加速度，U为势能函数，r_i为位置。
2. 有限差分积分算法
Verlet算法：
$$r(t + \Delta t) = 2r(t) - r(t - \Delta t) + a(t)\Delta t^2$$
- 优点：时间可逆，能量守恒好
- 缺点：速度计算精度低
Velocity Verlet算法：
$$v(t + \frac{\Delta t}{2}) = v(t) + \frac{a(t)\Delta t}{2}$$
$$r(t + \Delta t) = r(t) + v(t + \frac{\Delta t}{2})\Delta t$$
$$v(t + \Delta t) = v(t + \frac{\Delta t}{2}) + \frac{a(t + \Delta t)\Delta t}{2}$$
- 优点：同时更新位置和速度，精度高，稳定性好
- 是目前最常用的MD积分算法
时间步长选择：
- 典型值：2 fs（飞秒，10⁻¹⁵秒）
- 限制因素：体系中最快运动（键长振动，C-H键约10⁻¹⁴秒周期）
- 约束算法（SHAKE/SETTLE/LINCS）：固定键长，允许更大时间步长（2-4 fs）
3. 周期性边界条件（PBC）
- 模拟盒子（如立方体、十二面体、截角八面体）被其镜像复制填满整个空间
- 原子离开盒子一侧时，从另一侧进入（保持原子数恒定）
- 最小镜像约定：每个原子只与最近的镜像相互作用
- 目的：消除边界效应，模拟无限大体系
4. 长程相互作用处理
截断法（Cutoff）：
- 只计算截断距离（如10-12 Å）内的非键相互作用
- 简单高效，但对于带电体系可能引入显著误差
Particle Mesh Ewald（PME）：
- 将长程静电相互作用分解为实空间短程和倒空间长程两部分
- 实空间：截断范围内直接计算
- 倒空间：使用快速傅里叶变换（FFT）在网格上计算
- 优点：精度高（~10⁻⁵），效率合理
- 缺点：需要周期性边界条件