4. 分子动力学模拟基础：力场、周期性边界条件、积分算法简要回顾
做自由能估计那年，我栽过最大的跟头——不是方法不对，是分子动力学模拟的基底没打牢。你想想看，如果连原子怎么受力、边界怎么处理、轨迹怎么积分都含糊，那后面算出来的自由能差，说白了就是一堆有组织的随机数。
所以这一章，咱们把三个最基础的东西捋一遍：力场、周期性边界条件和积分算法。不搞长篇推导，只讲你在实际项目中必须理解的那层。
4.1 力场：原子间怎么“感知”彼此
力场就是描述原子间相互作用的数学函数加参数集。没了它，分子动力学模拟就是一盘散沙。
我个人的习惯是，拿到一个新体系，先问自己三个问题：
- 这个体系里起主导作用的是键能还是非键能？（比如蛋白质折叠，非键作用占大头）
- 极化效应重不重要？（比如金属离子体系，传统固定电荷力场容易翻车）
- 计算量能不能接受？（全原子 vs. 粗粒化，你得有数）
经典的力场形式长这样：
E_total = E_bond + E_angle + E_dihedral + E_nonbond
E_bond = Σ k_r (r - r0)²
E_angle = Σ k_θ (θ - θ0)²
E_dihedral = Σ k_φ [1 + cos(nφ - δ)]
E_nonbond = Σ ε_ij [ (Rmin_ij/r_ij)¹² - 2(Rmin_ij/r_ij)⁶ ] + Σ q_i q_j / (ε r_ij)键能、角能、二面角能描述分子内几何约束；非键能包括范德华和静电。我在项目中遇到过的最常见问题，就是参数文件选错了版本——比如用 CHARMM27 的蛋白参数配了 CHARMM36 的水模型，结果水与蛋白的相互作用直接偏掉，自由能估计自然不准。
力场选择小建议
- 配体-蛋白体系：GAFF + AMBER 力场是经典组合，但要注意配体的原子类型是否完整。
- 膜蛋白：CHARMM36 对脂双层参数支持很好。
- 高极化体系：AMOEBA、Drude 力场可以考虑，但计算成本是传统力场的 10 倍以上。
个人经验：我曾经在模拟金属蛋白时用了固定电荷力场，Zn²⁺ 周围的水结构完全错乱。换成极化力场后，配位数和自由能曲线才恢复正常。嗯，这种坑踩过一次就不会再犯。
4.2 周期性边界条件：打破“盒子墙”
模拟盒子就这么大，原子跑到边界怎么办？反射回去？太假。真实溶液里哪有墙。
周期性边界条件（PBC） 的做法很简单：把盒子在三维空间里无限复制。粒子从一边出去，就从对面那一边进来。这样表面的原子其实都在跟相邻副本的原子相互作用，盒子不再有“边界”。
但这里有个坑——最小镜像约定（Minimum Image Convention）。你只能让原子与 最近的那个副本 相互作用，否则同一个粒子会跟自己跨盒子的副本算两次，能量直接爆炸。
# 伪代码示意：应用最小镜像约定
for i in range(N_atoms):
for j in range(i+1, N_atoms):
dx = x[j] - x[i]
dy = y[j] - y[i]
dz = z[j] - z[i]
# 对每个维度：如果距离超过盒子一半，就减去盒子长度
if dx > 0.5 * box_x: dx -= box_x
if dx < -0.5 * box_x: dx += box_x
if dy > 0.5 * box_y: dy -= box_y
if dy < -0.5 * box_y: dy += box_y
if dz > 0.5 * box_z: dz -= box_z
if dz < -0.5 * box_z: dz += box_z
# 然后用 dx, dy, dz 计算力和能量注意：截断半径绝对不能超过盒子边长的一半。我见过有人在 40 Å 的盒子里设了 25 Å 的截断半径——那块截断区域直接跟自己的镜像重叠了，算出来的自由能偏得离谱。除非你用了 Ewald 或 PME 做长程静电，否则截断半径一定小于 L/2。
4.3 积分算法：时间怎么“走”
牛顿第二定律大家都知道：F = m a。有了力场算出力，有了力算出加速度，有了加速度更新速度和位置。问题是，这个更新要在离散的时间步里做。积分算法就是干这个的。
最常用的算法叫 Verlet 算法，以及它的变种蛙跳（Leapfrog）和速度 Verlet（Velocity Verlet）。
- Verlet：用 t-Δt 和 t 时刻的位置外推出 t+Δt 的位置，不显式计算速度。简单，但计算温度或压力时需要额外推速度。
- 蛙跳 Leapfrog：在整数时间步更新位置，在半整数时间步更新速度。名字很形象——速度和位置互相“跳过”对方。我早期用 GROMACS 时，默认就是蛙跳。
- 速度 Verlet：一步到位算出位置和速度，数值稳定性好，处理恒温恒压约束时也方便。我现在自己做自由能计算，几乎都用它。
时间步长选多少？经典规则是：
| 体系 | 典型时间步长 | 说明 |
|---|---|---|
| 全原子蛋白-水体系 | 1–2 fs | C-H 键振动周期约 10 fs，步长超过 2 fs 就会开始出错 |
| 带约束的体系 | 2–4 fs | 用 LINCS 或 SHAKE 冻结 X-H 键的高频振动 |
| 粗粒化体系 | 10–40 fs | 自由度少，运动更平滑 |
避坑指南：我曾经为了跑得更快，把时间步长从 2 fs 提到 4 fs，还觉得加了 LINCS 约束就万事大吉。跑了 100 ps 后发现某些二面角翻转异常频繁，能量也慢慢往上飘。检查日志才发现，那套力场对某些二面角的扭转势很陡，4 fs 的步长根本跟不上变化节奏。后来老老实实改回 2 fs，一切正常。
好，这三样东西是分子动力学模拟的基石。搞清楚力场怎么选、PBC 怎么设、积分算法怎么跑，后面再做自由能估计，至少心里有谱。下一个小节开始，咱们要真正动手搭模拟流程了——但在这之前，把这章的内容过一遍，能省下你至少一周的 debug 时间。