VASP分子动力学模拟：一个高温下MgO熔体结构的AIMD分析

提到分子动力学模拟，很多人首先想到的是GROMACS、LAMMPS这些基于经验势的经典MD软件。但在某些体系中，经验势的局限性会成为精度瓶颈——比如高温熔体中离子的极化效应、化学反应过程中的键断裂与形成、金属表面的氧化过程。这些场景下，VASP的ab initio分子动力学（AIMD）就成了必要的工具。

项目背景：为什么经验势不够用

团队处理过一个MgO高温熔体的局部结构分析项目。MgO是地幔矿物学中的关键组分，客户在做高温高压下硅酸盐熔体的物理性质研究，需要MgO在3000 K熔融状态下的结构数据作为基线参照——包括Mg-O配位数、O-Mg-O键角分布、以及熔体的短程有序程度。

经典MD模拟MgO熔体通常使用Born-Mayer势或Buckingham势，这些势在描述室温下的MgO晶体结构时精度不错。但在3000 K的熔融状态下，离子的极化效应显著增强——经验势中固定的离子电荷（Mg²⁺、O²⁻）无法描述电荷重新分布，导致配位数的预测出现偏差。文献中经验势给出的Mg-O配位数约4.8-5.2，而AIMD的结果约5.6-6.0——这个差异在高精度研究场景下是不可忽略的。

AIMD的计算参数

团队使用VASP进行AIMD模拟，PAW赝势，PBE泛函。模拟超胞为3×3×3（216个原子），截断能500 eV。k点只取Gamma点——AIMD模拟通常只用Gamma点，原因是在高温下电子态的热展宽已经抹平了k点采样带来的能量差异。

时间步长取2 fs——这是AIMD中比较激进的设置，标准做法通常是1 fs。但MgO体系中不存在轻原子（如H），且高温下离子振动频率相对较低，2 fs的能量守恒在测试中表现良好（总能量漂移<1 meV/atom/ps）。

温度控制使用Nose-Hoover链式恒温器，目标温度3100 K（略高于MgO的熔点3098 K，确保体系完全熔融）。先在NVE系综下跑2 ps预热，再切换到NVT系综正式模拟20 ps。

计算成本：AIMD的”贵”是真实的

216个原子的AIMD在64核（2×Intel Xeon Gold 6248）上的速度约为0.3 ps/天。20 ps的总模拟时间需要约67天——这对大多数项目来说是不现实的。

团队最终采用了并行策略：在3个不同的计算节点上同时运行3个独立的模拟（不同的随机初始速度），每个运行10 ps，最后合并统计结果。这种方式虽然不能增加单次模拟的统计时间，但通过多个独立轨迹的叠加可以提高统计可靠性。即便如此，总计算时间也需要约35天。

偏径向分布函数的结构信息

Mg-O的偏径向分布函数（PDF）在1.95 Å处出现第一峰，峰位比室温晶体中的Mg-O键长（2.10 Å）缩短了约0.15 Å——这是高温下热膨胀导致的平均距离变化。第一配位层的积分结果给出Mg-O配位数为5.8±0.4，高于经验势的4.8-5.2，更接近中子衍射实验外推值5.6-6.0。

O-Mg-O键角分布的峰值在82°和170°附近——前者对应于面共享的八面体连接，后者对应于边共享。晶体MgO中O-Mg-O键角的分布高度集中在90°（立方晶系），而熔体中键角分布明显宽化，出现了约30%的键角偏离90°超过20°的情况，说明八面体的畸变程度很高。

客户用AIMD的配位数数据修正了他们硅酸盐熔体热力学模型中的MgO活度系数——这个修正对预测岩浆中橄榄石的结晶温度有直接影响。AIMD的计算成本确实高昂，但对于经验势无法可靠描述的体系，它是获取精确结构信息的唯一途径。

如需了解更多的科研计算案例与项目经验，可以参考 https://www.keyanxueshu.com/ 上的技术文章。