分子动力学模拟基础：Lennard-Jones势函数与Verlet数值积分算法

引言

分子动力学（Molecular Dynamics, MD）通过数值积分原子体系的运动方程，追踪每个原子在飞秒至纳秒尺度上的轨迹，进而用统计力学系综平均提取扩散系数、热导率、弹性常数等宏观物性。它是连接微观相互作用与宏观连续介质模型的桥梁，广泛应用于材料设计、纳米摩擦、生物大分子构象分析等领域。一个 MD 程序的三个支柱是：描述原子间相互作用的势函数、时间积分数值算法、以及周期边界与近邻搜索等基础架构。本文逐一剖析这些要素，并给出二维 Lennard-Jones 流体的完整 Python 实现。

理论基础

运动方程与 Lennard-Jones 势

含 $N$ 个原子的体系遵循牛顿运动方程：

$$m_i \frac{\mathrm{d}^2 \mathbf{r}_i}{\mathrm{d} t^2} = \mathbf{F}_i = -\nabla_{\mathbf{r}_i} U(\mathbf{r}_1, \dots, \mathbf{r}_N)$$

最经典的 Lennard-Jones（LJ）12-6 势用排斥项刻画泡利不相容、用吸引项刻画色散力：

$$U(r) = 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]$$

其中 $\varepsilon$ 为势阱深度，$\sigma$ 为零势能距离，最小势能点位于 $r_{\min} = 2^{1/6}\sigma$。对应的力为：

$$\mathbf{F}(r) = \frac{24\varepsilon}{r^2}\left[2\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]\mathbf{r}$$

实际计算在截断半径 $r_c$（通常 $2.5\sigma$–$3.5\sigma$）处截断并做能量平移，配合最小镜像约定与周期边界条件，把力计算复杂度从 $O(N^2)$ 中可交互对数压缩到每个原子仅需处理几十个近邻。MD 常采用约化单位：以 $\sigma$、$\varepsilon$、原子质量 $m$ 为长度、能量、质量单位，时间单位为 $\tau = \sigma\sqrt{m/\varepsilon}$，对氩而言 $\tau \approx 2.16$ ps。

不同物理体系的势函数

势函数的选择决定了模拟能描述的物理——两体势无法给出金属的空位形成能，也无法描述共价键的方向性。

势函数	相互作用形式	适用体系	化学反应	相对成本
Lennard-Jones	两体对势	惰性气体、简单流体	否	1
Morse 势	两体对势	双原子分子振动	否	1
EAM 嵌入原子势	对势 + 电子密度嵌入能	金属与合金	否	3–5
Tersoff 键序势	三体键序	Si、C 等共价材料	有限	10
ReaxFF 反应力场	动态键级	燃烧、催化、氧化	是	50–100

速度 Verlet 积分器

对哈密顿系统，时间积分器必须保辛（symplectic）才能避免长期能量漂移。速度 Verlet 格式为：

$$\mathbf{r}(t+\Delta t) = \mathbf{r}(t) + \mathbf{v}(t),\Delta t + \frac{1}{2}\mathbf{a}(t),\Delta t^2$$

$$\mathbf{v}(t+\Delta t) = \mathbf{v}(t) + \frac{1}{2}\left[\mathbf{a}(t) + \mathbf{a}(t+\Delta t)\right]\Delta t$$

该格式时间可逆、辛、单步局部误差 $O(\Delta t^4)$（位置），全局二阶，且能量误差有界不随时间累积——这正是 MD 普遍弃用同阶 Runge-Kutta 的原因。步长由体系最快振动模决定：含氢体系需 $\Delta t \approx 1$ fs，全重原子体系可取 2–5 fs。温度由均分定理给出 $\langle K \rangle = \frac{f}{2} k_B T$，其中 $f$ 为自由度数；NVT 系综需借助 Berendsen 或 Nosé-Hoover 恒温器。

数值方法与实现

以下代码实现二维 LJ 流体的 NVE 模拟：FCC 晶格初始化、最小镜像约定、$O(N^2)$ 暴力对力计算（教学规模 $N=64$ 足够快），速度 Verlet 积分，输出温度与总能量守恒情况：

import numpy as np

def lj_md(N=64, rho=0.6, T0=1.0, dt=0.004, nsteps=3000, seed=42):
    """二维 Lennard-Jones 流体 NVE 分子动力学（约化单位）"""
    rng = np.random.default_rng(seed)
    n = int(np.ceil(np.sqrt(N)))
    L = np.sqrt(N / rho)
    xs = np.linspace(0, L, n, endpoint=False)
    pos = np.array([[x, y] for x in xs for y in xs][:N])
    pos += 0.05 * rng.standard_normal((N, 2))         # 扰动打破对称
    vel = rng.standard_normal((N, 2))
    vel -= vel.mean(axis=0)                            # 去除整体动量
    vel *= np.sqrt(T0 / (0.5 * (vel**2).sum() / N))    # 标定到目标温度
    rc2 = (0.45 * L) ** 2                              # 截断半径 < L/2

    def forces(pos):
        F = np.zeros_like(pos)
        U = 0.0
        for i in range(N - 1):
            dr = pos[i] - pos[i + 1:]
            dr -= L * np.round(dr / L)                 # 最小镜像约定
            r2 = (dr**2).sum(axis=1)
            mask = r2 < rc2
            inv_r2 = 1.0 / r2[mask]
            inv_r6 = inv_r2**3
            fij = 48.0 * inv_r6 * (inv_r6 - 0.5) * inv_r2
            fi = fij[:, None] * dr[mask]
            F[i] += fi.sum(axis=0)
            idx = np.nonzero(mask)[0] + i + 1
            np.add.at(F, idx, -fi)                     # 牛顿第三定律
            U += (4.0 * inv_r6 * (inv_r6 - 1.0)).sum()
        return F, U

    F, U = forces(pos)
    E0 = U + 0.5 * (vel**2).sum()
    for _ in range(nsteps):
        pos += vel * dt + 0.5 * F * dt**2
        pos %= L                                       # 周期边界
        F_new, U = forces(pos)
        vel += 0.5 * (F + F_new) * dt
        F = F_new
    T_inst = 0.5 * (vel**2).sum() / N                  # 2D: <K>/N = kT
    E_tot = U + 0.5 * (vel**2).sum()
    return T_inst, E_tot / N, abs(E_tot - E0) / abs(E0)

T, e_per_atom, drift = lj_md()
print(f"平衡温度 T* = {T:.3f} (非平衡晶格弛豫释放势能，略高于初始 1.0)")
print(f"单原子能量 E/N = {e_per_atom:.3f}, 能量相对漂移 = {drift:.2e}")

运行可见：体系从受扰晶格初态弛豫，非平衡势能释放转化为动能，温度在初期波动后稳定在 $T^* \approx 1.05$；数千步内总能量相对漂移保持在 $10^{-3}$ 以内（本例约 $4\times10^{-4}$），体现辛积分器的长期稳定性。把 $N$ 提高到数千时，应将暴力对力换成 Verlet 近邻表与元胞链表，复杂度降至 $O(N)$。近邻表的实现要点是：为每个原子记录 $r_c + r_s$（表皮半径，通常取 $0.3\sigma$–$0.5\sigma$）内的邻居，仅当任一原子位移超过 $r_s/2$ 时才重建表，从而在力计算不变的精度下把建表频率降低一到两个数量级。另一个常被忽视的细节是截断势能的一致性处理：若直接在 $r_c$ 处硬性截断而不做能量平移（shifted potential），势能面出现台阶，粒子穿越截断面时受到脉冲力，能量守恒显著恶化；更严格的做法是再加力平移（force-shifted），让势能与力在 $r_c$ 处同时光滑归零。

工程案例分析

问题：液态氩自扩散系数估算

以氩的 LJ 参数（$\varepsilon/k_B = 120$ K，$\sigma = 0.34$ nm）模拟液氩在 100 K、近饱和密度下的自扩散系数。扩散系数由均方位移的 Einstein 关系给出 $D = \lim_{t\to\infty} \langle|\mathbf{r}(t)-\mathbf{r}(0)|^2\rangle/(2d,t)$，$d$ 为空间维数。

状态点	$T^*$	$\rho^*$	$D^*$（典型值）	换算 $D$ (10⁻⁹ m²/s)
过冷液相	0.72	0.85	≈0.012	≈0.9
常温液相	1.00	0.80	≈0.032	≈2.4
高温液相	1.50	0.70	≈0.075	≈5.6

结果与实验值（100 K 下约 $2\times10^{-9}$ m²/s）处于同一量级，验证了 LJ 势对简单流体的定量预测能力。工程上更复杂的体系——水溶液、聚合物、金属熔体——则必须转向更精细的全原子力场（OPLS、AMBER）或 EAM 类势函数。

结论与展望

本文覆盖了 MD 的最小完备知识集：势函数决定物理保真度，速度 Verlet 提供辛且长期稳定的积分，周期边界与近邻表保证计算可行。进阶议题包括：约束算法（SHAKE/LINCS）冻结氢键振动以放大步长、Ewald/PPPM 处理长程静电、多时间步 RESPA 积分、以及恒温恒压系综的严格实现。当前最活跃的变革来自机器学习势函数——以第一性原理数据训练的神经网势头（如 Behler-Parrinello 型、图神经网络型）能以经典力场百倍左右的成本逼近 DFT 精度，正在把 MD 的适用边界从”千万原子纳秒”推向”亿级原子微秒”。