分子动力学模拟势函数

我做过一个势函数选型项目，客户的问题是：模拟金属、聚合物和离子晶体能不能用同一种势？答案显然是否定的。我这个项目我把几种常见势函数放在一起对比，帮客户建立了选型思路。

左侧图是三种对势函数：Lennard-Jones（蓝色）、Morse（红色虚线）、Buckingham（绿色点线）。横轴是原子间距 r，纵轴是势能 U。LJ 势是最简单的 van der Waals 模型，U = 4ε[(σ/r)¹² − (σ/r)⁶]，在 r=σ 处 U=0，在 r=2^(1/6)σ 处取最小值 −ε。它的优点是形式简单、计算快，适合稀有气体和弱相互作用体系；缺点是不能描述化学键的方向性和电荷转移。

Morse 势在 LJ 基础上加了一个指数排斥项，U = D_e[1 − e^(−a(r−r_e))]² − D_e。它更适合描述共价键或金属键，因为能反映键的解离行为。Buckingham 势 U = A·e^(−Br) − C/r⁶，排斥项是指数形式，比 LJ 的 r⁻¹² 更接近量子力学结果，常用于离子晶体（如 MgO、SiO₂）。

右侧图是 EAM（Embedded Atom Method）嵌入函数，也就是多体势。横轴是电子密度 ρ，纵轴是嵌入能 F(ρ)。EAM 认为每个原子嵌入在周围原子产生的电子气中，总能量包括嵌入能和两体静电排斥。这种多体效应让 EAM 能描述金属的弹性、缺陷和表面性质，比纯对势更准。我做金属体系时，基本首选 EAM 或 MEAM；做合金则常用 MEAM 或机器学习的 SNAP/ACE 势。

踩坑经验：第一个是参数来源。LJ 参数可以从 OPLS-AA、AMBER、TraPPE 等全原子力场拿；EAM 参数要从文献或 NIST 数据库找对应元素组合，不要乱用。第二个是截断半径：对势通常在 10-12 Å 截断，EAM 的截断要和势函数文件匹配。第三个是混合规则：做异质体系时，不同原子的 LJ 交叉参数常用 Lorentz-Berthelot 混合规则，但金属-有机界面这种组合，混合规则经常不准，需要专门拟合。

工程经验：势函数选择直接决定模拟结果的可靠性。我一般会先问三个问题：体系含什么元素？要研究什么性质？有没有现成的验证数据？如果是金属缺陷/力学，EAM；如果是有机小分子/聚合物，GAFF/OPLS-AA；如果是离子晶体，Buckingham；如果是反应过程，ReaxFF 或机器学习势。没有万能势函数，只有合适势函数。

更多 LAMMPS 势函数设置和参数来源的实战经验，我汇总在 [分子动力学模拟势函数](https://www.keyanxueshu.com/category/md/) 栏目里。[分子动力学模拟势函数](https://www.keyanxueshu.com/)

势函数选型的最后一个检查是迁移性（transferability）。同一个 LJ 参数从纯水迁移到水-蛋白质界面，可能表现很差，因为界面处的电子环境不同。我的做法是在选定势函数后，做一个小的 DFT 对比：取几个代表性构象，分别用 DFT 和力场算能量，看能量面趋势是否一致。如果力场严重低估或高估某些构象能量，就要换参数或换势函数。对于金属体系，EAM 势的好坏可以通过算空位形成能、层错能和表面能来验证。把这些验证数据放在报告里，客户会相信你的模拟不是“Garbage in, garbage out”。