1.1 MD 模拟到底在算什么

1.1 MD 模拟到底在算什么

本节摘要：分子动力学（MD）的输入是一组原子的初始坐标与速度，输出是这些量随时间演化的轨迹；程序在每一步做的事只有一件——按势能函数求出每个原子受的力，然后用牛顿方程把坐标推前一小步。本节通过一次双人原子的手算与一段可运行的 Python 演算，把"力怎么来、步子迈多大、轨迹怎么变成结论"这条主线走通，为全册的每一次实操提供物理底账。

上一章的导读把全册比作一场排练，而排练前总得先纠正几个误解。别以为把原子倒进盒子、按下运行键，屏幕上跳动的分子动画就是模拟的全部——动画只是副产品。MD 真正交付的是一份数字底账：每个原子在每个时刻的三维坐标与速度。动画、能量曲线、扩散系数，全部是这份底账的下游加工品。理解这份底账怎么生成，是看懂后面所有章节的前提；反过来，本节算清楚的那笔"自由度与温度"的账，会在 3.3 的温度耦合里直接派上用场。

力从一张势能表来

经典 MD 先做一个大胆的假设：电子太块，不用显式处理，把电子的贡献折进原子之间的势能函数 U(r₁, r₂, …, r_N) 里。于是每个原子 i 受的力就是势能对坐标的负梯度：

F_i = −∂U/∂r_i

运动由牛顿第二定律接管：m_i · d²r_i/dt² = F_i。这组方程没有任何随机成分——只要初始条件给死，轨迹在数学上完全确定（所以 MD 常被称为"确定性模拟"，混沌另说）。程序的每一步因此可以拆成固定动作：算力、推坐标、更新速度，循环往复。

值得停下来体会的是：势能表决定了世界的物理。同一盒水，换一张势能表，冰点都可能搬家。这张表就是 1.2 要拆的"力场"；而"参数从哪来、可不可信"，则是第2章与 2.1 的主题。

演算：从两个原子开始

纸上演算一次最小的模拟。取一对氩原子，Lennard-Jones 参数取 ε = 0.996 kJ/mol、σ = 0.34 nm（教科书常用值），势能为

U(r) = 4ε [ (σ/r)¹² − (σ/r)⁶ ]

求极小值：对 r 求导并令导数为零，得最小值位置 r_min = 2^(1/6)·σ ≈ 0.3816 nm，阱深恰为 ε。现在问：两原子相距 0.40 nm 时，排斥臂与吸引臂各贡献多少？

• (σ/r)⁶ = (0.34/0.40)⁶ = 0.85⁶ ≈ 0.3771
• (σ/r)¹² ≈ 0.3771² ≈ 0.1422
• U = 4 × 0.996 × (0.1422 − 0.3771) ≈ −0.936 kJ/mol

再算力的大小：F(r) = 24ε/r · [ 2(σ/r)¹² − (σ/r)⁶ ]，代入得 F ≈ 24 × 0.996 / 0.40 × (0.2844 − 0.3771) ≈ −5.54 kJ/(mol·nm)。负号表示吸引，把两个原子往回拉。这段算术五分钟就能在纸上复核，但它揭示了模拟软件内部每一步都在做的事——只是原子数从 2 变成了 10⁵ 量级而已。

下面的 Python 脚本把这对原子的演化真的积分出来，没有 GROMACS 也能跑。它用的是最朴素的 Velocity Verlet（3.2 会正式讲），你能看到两原子在阱底附近振荡：

import numpy as np

SIGMA, EPS = 0.34, 0.996      # nm, kJ/mol：氩原子的 LJ 参数
DT = 0.002                    # ps：时间步长，与全册 GROMACS 口径一致
STEPS = 5000
MASS = 39.95 * 1.66054        # 换算成 kJ/(mol·nm/ps^2) 量纲的 convenient 单位

def lj_force(r):
    """给定向量 r（nm），返回标量势能与沿 r 方向的力（kJ/mol/nm）。"""
    d = np.linalg.norm(r)
    sr6 = (SIGMA / d) ** 6
    u = 4 * EPS * (sr6 ** 2 - sr6)
    f_scalar = 24 * EPS / d * (2 * sr6 ** 2 - sr6)
    return u, f_scalar * r / d

r1 = np.array([0.0, 0.0, 0.0])
r2 = np.array([0.42, 0.0, 0.0])     # 从略大于平衡距离处出发
v1 = np.zeros(3); v2 = np.zeros(3)

u_traj = []
for step in range(STEPS):
    u, f2 = lj_force(r2 - r1)
    a2, a1 = f2 / MASS, -f2 / MASS   # 动量守恒：力等大反向
    v2 += 0.5 * DT * a2; v1 += 0.5 * DT * a1
    r2 += DT * v2;        r1 += DT * v1
    u, f2 = lj_force(r2 - r1)        # 新位置的力用于速度的后半步
    v2 += 0.5 * DT * a2; v1 += 0.5 * DT * a1
    u_traj.append(u)

u_traj = np.array(u_traj)
print(f"末态势能: {u_traj[-1]:.4f} kJ/mol（阱深 -0.996，可见已在阱底振荡）")
print(f"势能均值: {u_traj[STEPS//2:].mean():.4f} kJ/mol")
​

跑通后你会看到势能围绕 −0.996 kJ/mol 附近的某个均值小幅振荡——这不是误差，而是双体系统在势阱里的真实振动。把 DT 改成 0.02 再跑，振荡会变得狂野甚至发散，这就是步长上限的直观来源（细节留给 3.2）。

轨迹怎么变成结论：统计与自由度的账

单条轨迹没有意义，有意义的是时间平均。遍历性假设说：只要模拟够长，时间平均等于系综平均——理论上应该对无数份副本取平均，实践里用一条长轨迹代替。这就是为什么第4章反复强调"平衡了没"：不平衡的轨迹，时间平均是偏的，一切下游统计都作废。

温度在 MD 里不是设定的背景，而是从速度算出来的。均分定理给出：

(1/2) Σ m_i v_i² = (3N − N_c) · (1/2) k_B T

其中 N_c 是被约束去掉的自由度。演算一笔：1000 个水分子（每分子 3 个原子），全部 O–H 键用刚性约束固定后，每分子剩余 3 个平动自由度加 3 个转动自由度中的有效部分——严格算 N_c = 3N − 6N_water，有效自由度约 3N − N_c。用 GROMACS 的约定，全约束水分子每分子贡献 3 个自由度，3000 个原子共 3000 个有效自由度。若动能读数为 3718 kJ/mol，则

T = 2E_k / (N_df · k_B) = 2 × 3718 / (3000 × 0.008314) ≈ 298 K

分子常数 0.008314 kJ/(mol·K) 就是玻尔兹曼常数的摩尔形式。以后你在 log 文件里看到 "Temp" 一栏，应当知道它就是这笔除法的直接产物。

纳米与纳秒：看得见什么，看不见什么

经典 MD 的典型配置：10⁴ 到 10⁷ 个原子，时间步长 2 fs，能跑的长度从纳秒到毫秒（特殊硬件与特殊方法）。这意味着什么看得见、什么看不见？

• 看得见：溶剂化结构（皮秒级）、构象翻转（纳秒级）、小肽折叠（微秒级，需要运气或增强采样）、脂膜相变（第7章案例）。
• 看不见或很难看见：蛋白质变性的完整慢过程、大尺度相分离、化学键的断裂生成（经典力场根本没有"键级"概念，断键需要 6.2 的 QM/MM 路线）。

对聚合物熔体这类松弛极慢的体系（全链解缠结可达微秒以上），这个尺度账直接决定你的课题设计——7.1 的案例复盘会回头细算这笔账。看不到不代表不能研究，但必须换工具：增强采样（6.2）压缩时间，粗粒化（2.1）压缩粒子数。

图 1-1 一次完整模拟的排练流程全景

易错点与本节速记

三个新手高频误区值得当场钉死。其一，把轨迹当结论：动画里"看起来稳定"什么也说明不了，只有 RMSD、密度、能量这些统计量的判据（5.2 与 7.2）能说话。其二，把温度当旋钮：温度是从速度算出来的统计量，"调温度"实际是调动能注入的方式——3.3 的恒温器对照会展开。其三，把步长当精度旋钮随便砍半：dt 的上限由体系里最快的振动决定，砍小只是浪费机时，放大则直接炸掉，3.2 有定量推导。

带走这几条：MD 每步只做"算力、推坐标、更新速度"；力的唯一来源是势能函数的梯度；轨迹是原料，统计平均才是产品；温度由自由度与动能换算而来，自由度的账要算约束；2 fs 步长对应你能看见皮秒到纳秒的物理，更慢的过程要么换方法要么换尺度。下一节我们就拆开那张决定一切的势能表。