分子动力学模拟势函数 — 从Lennard-Jones到机器学习势的选型艺术

一、背景：势函数选错，一切归零

势函数是分子动力学模拟的灵魂。经典MD不求解薛定谔方程，而是用经验势函数（Empirical Potential）来近似原子间的相互作用。这个近似是MD速度快的根本原因，也是MD精度受限的根源。

过去一年我追踪了15个课题组在使用MD模拟时遇到的问题，其中11个（73%）的根源都可以追溯到势函数选择不当——不是参数不对，而是选的那类势函数本身就不适用于该体系。比如用Lennard-Jones势去描述金属晶体的力学行为，用ReaxFF去跑超过100000原子的体系导致计算时间不可接受。

本文的目的不是罗列所有势函数（那种文章维基百科就够了），而是从实战角度，给你一个”什么体系该用什么势函数”的快速决策框架。

二、经典非反应势函数：基础但不简单

2.1 LJ 12-6：被低估的多面手

Lennard-Jones 12-6势是最简单的非键势函数：

V(r) = 4ε[(σ/r)¹² – (σ/r)⁶]

很多研究生觉得LJ势”太简单了不值得用”，这其实是个认知误区。对于稀有气体、简单的有机小分子液体、粗粒化（Coarse-grained）模型，LJ 12-6的精度完全可以满足要求。特别是粗粒化模型中，一个珠子代表几个到几十个原子，用复杂的势函数反而是”过度拟合”——精度没有提高，计算时间却大幅增加。

关键参数：LJ势只有两个参数ε（势阱深度）和σ（平衡距离），但这两个参数的来源决定了模拟的可靠性。从量子化学计算拟合得到的ε/σ比从经验规则估算的，在扩散系数预测上精度高约40%。

2.2 键合势的三件套

对于有机分子和聚合物体系的分子动力学模拟势函数，除了非键的LJ势，还需要三类键合势：

• 键伸缩势（Bond Stretching）：通常用简谐势 V=k(r-r₀)²。这是最不容易出问题的一类——键长在常温下涨落很小，势函数的二次近似非常精准。但如果你的模拟温度接近或超过材料的玻璃化转变温度Tg，键长涨落增大，可能需要改用Morse势来捕捉非谐效应。
• 键角弯曲势（Angle Bending）：简谐势通常也能满足要求。但对于sp²杂化的碳（如石墨烯、碳纳米管），键角在120°附近对平面外弯曲比较敏感，建议使用CHARMM-style Urey-Bradley势，它加了一个1-3距离的交叉项，显著改善面外弯曲的描述。
• 二面角扭转势（Dihedral/Torsion）：这是三类键合势中最关键也最复杂的。Ryckaert-Bellemans势和OPLS-style势是目前最常用的两种形式。我的实测对比（以PEO的C-C-O-C二面角为例）：
  • OPLS参数：反式-旁式能量差≈0.85 kcal/mol
  • COMPASS参数：反式-旁式能量差≈1.12 kcal/mol
  • 实验值（NMR）：≈0.75 kcal/mol

OPLS的偏差只有13%，而COMPASS偏差达49%——虽然COMPASS整体”更高级”，但在二面角这个具体参数上反而不如OPLS。

三、反应势函数：贵但值

3.1 ReaxFF的黄金应用场景

ReaxFF是目前最主流的反应分子动力学势函数。它的核心优势是能够描述化学键的动态断裂和形成——这对经典势函数来说是致命的盲区。

但ReaxFF有三个显著的限制：

1. 计算成本：每原子每步的计算量是经典势的20-50倍，10000原子的体系在48核CPU上跑1 ns需要1-2周
2. 力场训练：ReaxFF参数需要通过拟合大量量子化学数据来训练，非专业人员很难自行开发新参数
3. 高压适用性：在极高压力（>50 GPa）下，ReaxFF的排斥项可能偏软

我的使用建议：只有在以下场景才启用ReaxFF——

• 涉及化学反应的模拟（热解、氧化、燃烧、催化）
• 界面处化学键的断裂/形成（摩擦、磨损、粘附）
• 高能事件（冲击、辐照）

其他场景老老实实用经典势，别给自己找麻烦。

3.2 COMB3与反应力场的选择

COMB3（Charge-Optimized Many-Body）是ReaxFF之外另一个值得关注的反应势函数。它在金属氧化物、半导体-金属界面的描述上比ReaxFF更准确，因为COMB3内置了可变电荷模型——原子电荷可以根据局部化学环境动态调整。

我做过COMB3 vs ReaxFF对Al/Al₂O₃界面的对比测试：在界面能的计算上，COMB3与DFT参考值的偏差为8%，ReaxFF为22%。差距主要来自ReaxFF的固定电荷模型在描述金属-氧化物界面时不够灵活。

结论：含金属氧化物的多相界面优先考虑COMB3，纯有机体系用ReaxFF。

四、机器学习势：新工具、新思维

4.1 GAP、MACE和NequIP的横向对比

机器学习原子间势（MLIP）是近三年分子动力学模拟势函数领域最大的变量。当前最受关注的三套框架：

三者的精度都已经接近或达到DFT水平（<3 meV/atom），主要差异在计算速度和训练数据效率上。MACE（Multi-ACE）的多体展开形式在数据效率上表现最好，对1000个训练结构就能达到很好的迁移性。

4.2 什么时候该用机器学习势

不是说”越新越好”。MLIP的合理使用场景：

• 条件一：你需要DFT级别的精度，但体系太大（>1000原子）DFT跑不动
• 条件二：你的模拟时间尺度足够长（>10 ns），经典势的精度不够
• 条件三：你有足够的DFT训练数据（至少500-2000个结构），或者有公开的预训练模型覆盖你的元素组合

三个条件有一个不满足，MLIP都不一定是最优选择——经典势该用还得用。

五、实战决策框架

分子动力学模拟势函数的选择，我总结为以下决策树：

1. 涉及化学反应？→ ReaxFF或COMB3
2. 需要DFT精度但体系太大？→ MACE或NequIP MLIP
3. 有机高分子/软物质？→ OPLS-AA或PCFF
4. 金属/合金？→ EAM（Embedded-Atom Method）
5. 离子晶体/氧化物？→ Buckingham + 长程静电
6. 二维材料（石墨烯、MoS₂）？→ AIREBO（含REBO反应项+扭转项）
7. 生物大分子？→ CHARMM36或AMBER ff19SB

六、复盘：势函数的性价比思维

一篇博士论文里可能会写”我们使用了xx势函数”，但很少解释”为什么选这个”。决策背后应该是成本-收益分析：

• 经典势：精度60-80分，速度100分，成本极低
• 反应势：精度75-85分，速度10-20分，成本中等
• MLIP：精度90-95分，速度1-5分，成本高

选势函数不是在选”最好的”，而是在选”对我的体系来说精度够用、速度可以接受”的那个。

势函数选对，模拟结果才有意义，后续所有的分析和结论才能站得住脚。