模拟计算分子动力学：从入门到精通的实战经验

分子动力学(MD)模拟是连接微观与宏观的桥梁——从原子运动轨迹出发，计算扩散系数、热导率、力学性质等宏观物理量。无论你用的是LAMMPS、GROMACS还是MS，MD的核心原理是相通的。本文从跨软件的通用视角分享MD模拟的实战经验。

一、MD模拟核心概念

1.1 基本原理

MD模拟的核心是牛顿运动方程：

m_i × d²r_i/dt² = F_i = -∇U

• m_i：原子i的质量
• r_i：原子i的位置
• F_i：原子i受到的力
• U：总势能（由力场计算）

1.2 积分算法

经验：99%的MD模拟用Velocity Verlet就够了。不要纠结积分算法选择。

二、MD模拟标准流程

2.1 七步法

Step 1: 建模 → 初始构型
Step 2: 力场分配 → 势能函数
Step 3: 能量最小化 → 消除重叠
Step 4: 平衡模拟 → NVT/NPT达到目标T/P
Step 5: 生产模拟 → 采样轨迹
Step 6: 后处理分析 → 物理量计算
Step 7: 验证 → 与实验/DFT对比

2.2 能量最小化

必做！ 跳过这一步直接跑MD几乎一定会爆炸。

经验：先用Steepest Descent跑1000步消除大的重叠，再切换Conjugate Gradient精弛豫到力<0.1 kcal/mol/Å。

2.3 平衡阶段判断

平衡的3个标志：

1. 温度稳定：波动<±10K（目标300K时）
2. 能量稳定：总能量波动<5%（NVE）/ 势能波动<5%（NVT）
3. 密度稳定：NPT中密度波动<2%

验证方法：画温度/能量/密度vs时间曲线，观察是否进入平台期。如果还在漂移，继续跑。

2.4 生产模拟

平衡后切换到产数据模式：

三、关键物理量计算经验

3.1 径向分布函数(RDF)

RDF描述原子间距的统计分布：

经验值参考：

• 液态水O-O第一峰：2.8 Å
• 金属Cu-Cu最近邻：2.56 Å
• 高分子C-C键长：1.54 Å

计算经验：

• 需要足够帧数(>1000)
• bin宽度：0.05-0.1 Å
• 截断距离：小于盒子最小边长的一半

3.2 扩散系数

两种方法：

MSD法： D = lim(t→∞) <|r(t)-r(0)|²> / (6t)

Green-Kubo法： D = (1/3) ∫<v(0)·v(t)> dt

经验：

• MSD法更直观，但需要线性区(通常>100ps)
• Green-Kubo法精度高但需要频繁采样速度
• 液体扩散系数典型值：10⁻⁹~10⁻¹⁰ m²/s
• 固体扩散系数：10⁻¹²~10⁻¹⁵ m²/s（需要长时间模拟）

踩坑经验：固体扩散系数极小，常规MD时间(1ns)内可能只看到1-2次跳跃事件。建议用增强采样(NEB/metasynamics)或高温加速。

3.3 热导率

NEMD法（非平衡MD）：

• 在体系两端设热源和热库
• 测量温度梯度
• κ = J / (dT/dx)

Green-Kubo法：

• κ = (1/(3kBV²T)) ∫<J(0)·J(t)> dt

经验：

• NEMD更直观但体系要大(>50×50 Å)
• Green-Kubo需要长轨迹(>1ns)和频繁采样
• 典型热导率：Cu 400 W/mK, Si 150 W/mK, 玻璃 1 W/mK

3.4 力学性质

应力-应变曲线：

1. NPT平衡到目标温度
2. 逐步施加应变(每步0.001)
3. 每步弛豫后记录应力
4. 绘制σ-ε曲线

弹性常数：

• 施加小应变(0.001)在不同方向
• 记录应力响应
• C_ij = σ_j / ε_i

四、不同软件的选择经验

4.1 软件对比

4.2 选择决策

你的体系是什么？
├─ 金属/合金/陶瓷 → LAMMPS
├─ 蛋白质/DNA/药物 → GROMACS
├─ 高分子/聚合物 → LAMMPS或GROMACS
├─ 化学反应 → LAMMPS+ReaxFF
├~ 纳米材料(石墨烯/CNT) → LAMMPS
└~ 教学/快速验证 → MS Forcite

五、常见陷阱与避坑

5.1 时间步长陷阱

5.2 盒子大小陷阱

最小盒子尺寸：大于截断距离×2

如果盒子太小：

• 原子通过周期边界与自身相互作用
• RDF在L/2处截断(L为盒子边长)
• 表面效应未被消除

5.3 统计误差陷阱

MD是统计力学方法，单次模拟结果有统计误差：

经验法则：

• RDF：需要>1000帧
• 扩散系数：MSD线性区>100ps
• 热导率：>1ns轨迹
• 自由能：需要充分采样

降低误差：多次独立模拟取平均，比一次长模拟更好。

5.4 力场验证陷阱

最常见的错误：选了力场不做验证直接跑

正确流程：

1. 验证晶格常数(<2%偏差)
2. 验证密度(<5%偏差)
3. 验证关键物理量(弹性常数/熔点)
4. 验证通过后再跑正式模拟

六、项目报价参考

结语

MD模拟的核心经验是：”先验证力场，再做平衡，最后产数据”。跳过任何一步都可能导致结果不可靠。建议新手从简单体系(如液态Ar或Cu体相)开始练习，掌握流程后再做复杂体系。如有MD计算需求，欢迎联系我们获取定制方案。