1. 项目概述当材料模拟需要兼顾精度与规模在计算材料科学这个行当里干了十几年我越来越深刻地体会到我们这些做模拟的人每天都在一个“不可能三角”里挣扎计算精度、体系规模和计算时间。你想用第一性原理密度泛函理论DFT去算一个完美的晶格常数或者弹性模量没问题精度可以做到和实验媲美。但你想用它去模拟一个包含几万个原子的纳米颗粒在拉伸下的位错形核过程或者模拟一个百万原子级别的多晶材料在辐照下的缺陷演化那基本等于让一台超级计算机去“思考人生”算到天荒地老也未必有结果。反过来传统的经验势函数比如嵌入原子法EAM或者修正的EAMMEAM计算速度飞快模拟百万原子轻轻松松但它的“经验”二字就决定了其局限性——它本质上是一个拟合出来的函数在拟合所用的数据范围内表现尚可一旦遇到训练数据之外的构型比如高能碰撞产生的复杂缺陷团、极端应力下的非晶态结构其可靠性就大打折扣预测结果可能和实际情况南辕北辙。所以过去十年我们整个领域都在寻找那个“圣杯”一种既能逼近DFT精度又能拥有经验势函数计算速度的方法。机器学习势函数MLP的出现让我们看到了曙光。它不再依赖于物理学家事先预设的函数形式而是像一个聪明的学生从海量的高精度DFT计算数据中“学习”原子间的相互作用规律。这其中高斯近似势Gaussian Approximation Potential, GAP是公认的标杆之一它的数学框架严谨预测精度非常高。但GAP也有个“阿喀琉斯之踵”每次计算能量和力时都需要进行复杂的核函数评估计算开销依然比经验势大得多限制了其在大规模体系中的应用。直到tabGAP表格化GAP的出现局面才被真正打破。我最初接触到这个概念时感觉就像给一辆高性能跑车GAP换上了更高效的变速箱和轮胎表格化与插值。它的核心思想非常巧妙既然GAP对能量的贡献可以分解为许多“基函数”的加权和而这些基函数的值只依赖于原子对的局部环境描述符比如原子间距、角度等那么我们何不事先把这些基函数在可能出现的所有描述符值上的结果都算好存成一张巨大的多维表格呢在实际模拟中只需要根据当前原子环境的描述符值去这张表格里做快速的插值查询就能立刻得到近似的能量和力。这个“预计算查表”的过程将原本昂贵的在线计算转化为了几乎可以忽略不计的内存访问和简单算术运算。我这次的工作就是基于这个思路为三种最经典的面心立方fcc金属——铝Al、铜Cu和镍Ni——构建并系统验证了它们的tabGAP势函数。选择它们是因为其数据丰富、物理性质研究透彻是检验新方法可靠性的绝佳试金石。我们的目标很明确第一验证tabGAP在静态性质晶格常数、弹性常数、空位形成能等上能否达到与原始GAP和DFT相当的精度第二也是更关键的检验它在真实的动态分子动力学模拟中特别是涉及非平衡、大变形、缺陷演化的复杂场景下是否依然稳健可靠并且真正实现计算效率的数量级提升。最终我们成功地将这套势函数应用于包含数百万原子的大规模单轴压缩模拟并系统研究了不同晶向的阈值位移能这些都是传统高精度方法难以企及的尺度。下面我就把这其中的门道、踩过的坑以及实操心得掰开揉碎了和大家聊聊。2. tabGAP势函数的构建从原理到实现2.1 核心思想用空间换时间化连续为离散要理解tabGAP得先简单回顾一下GAP是怎么工作的。GAP将体系的总能量表示为所有原子对或更高阶贡献的加和。对于每个原子对其能量贡献通过一个高斯过程回归模型来预测这个模型的输入是描述该原子对局部环境的“描述符”descriptor。常用的描述符比如SOAPSmooth Overlap of Atomic Positions它能将原子周围的邻居分布信息转化为一个高维的、旋转平移不变的向量。计算这个描述符以及后续的高斯过程核函数评估是GAP计算中最耗时的部分。tabGAP的妙招就在于“离散化”和“预计算”。我们仔细分析发现在分子动力学模拟中原子对的局部环境虽然千变万化但其描述符向量的取值实际上是落在一个高维空间中的连续区域里。如果我们把这个连续空间进行离散化网格划分那么任何一个实际出现的描述符向量都可以用其周围网格点上的值通过插值来近似。于是构建tabGAP的流程就清晰了确定描述符与网格首先选定用于构建势函数的描述符例如基于原子间距的两体描述符或包含角度的三体描述符。然后根据训练数据的分布范围为描述符的每一个维度确定其最小值和最大值并设定一个网格间距。这个间距决定了表格的精度和大小。间距越小插值越精确但表格体积会指数级增长维度灾难。这里就需要做权衡。生成网格点数据对于网格划分后产生的每一个网格点对应一个特定的描述符向量组合我们调用原始的、完整的GAP模型计算出在该描述符下原子对的能量贡献以及力的贡献如果需要。这是一个“离线”计算过程虽然计算量巨大但只需要做一次。创建插值表格将步骤2中计算出的所有网格点的能量和力值存储为一个多维数组这就是我们的“势函数表格”。同时我们还需要选择合适的插值算法。对于平滑变化的势函数线性插值通常就能在精度和速度之间取得很好的平衡对于变化剧烈的区域可能需要更高阶的插值如三次样条。集成与调用将这张巨大的表格和插值器集成到分子动力学软件如LAMMPS中。在模拟运行时对于每一对原子程序实时计算其当前的描述符向量然后通过查表插值瞬间获得其相互作用的能量和力完全绕过了复杂的GAP核函数计算。注意网格维度与内存的博弈。这是构建tabGAP时第一个要面对的挑战。如果你只用原子间距作为描述符一维那么表格就是一维数组非常小巧。但为了描述键角等环境信息我们常常需要使用二维甚至更高维的描述符。一个精细的4维网格其节点数可能是(N_grid)^4即使每个维度只取50个点总节点数也高达625万个。每个节点存储一个双精度浮点数8字节这张表就要占用约50MB内存。这还只是一个原子种类对如Al-Al的表格。对于多元素体系表格数量会组合增长。因此在构建时必须精心选择描述符的维度和网格密度在保证精度的前提下控制内存占用在可接受范围内例如单个势函数文件控制在几百MB以内。2.2 我们的实操为Al, Cu, Ni构建tabGAP基于上述原理我们为Al、Cu、Ni分别构建了tabGAP。这里分享一些关键的技术细节和决策背后的考量训练数据库的构建这是所有机器学习势函数的基石。垃圾进垃圾出。我们的数据库不仅包含了完美的晶体构型更重要的是纳入了大量“非理想”构型弹性变形对晶胞施加各种应变张量计算其能量变化用于学习材料的弹性响应。点缺陷包含空位、间隙原子、以及它们的小团簇。这是学习缺陷形成能的关键。表面与晶界创建了不同晶面的表面模型和典型的对称倾侧晶界用于捕捉表面能和界面能。液态与非晶态通过高温熔化和快速淬火的分子动力学轨迹使用一个较快的经验势预跑采集液态和可能非晶态的构型以描述无序结构中的相互作用。高能碰撞轨迹为了模拟辐照损伤我们特意用经典分子动力学模拟了一些高能级联碰撞事件从中采样了高能、高畸变的原子环境确保势函数能正确处理阈值位移能TDE计算中出现的极端情况。所有这些构型的能量和原子受力都使用高精度的DFTVASP软件PAW-PBE泛函进行计算作为训练的“标准答案”。数据库的规模和质量直接决定了势函数的泛化能力。描述符与网格参数选择我们主要采用了两体加三体的描述符组合。两体描述符基于原子间距负责描述随距离变化的排斥和吸引作用三体描述符基于两个键长和夹角则能捕捉键角偏好这对于描述fcc金属的稳定性、堆垛层错能等至关重要。对于两体部分截断半径cutoff设为6 Å左右确保包含足够多的邻居信息。网格间距在平衡距离附近设置得较密如0.02 Å在远离平衡距离的区域可以适当放宽以节省存储空间。对于三体部分维度更高我们采用了更激进的策略使用稀疏网格或降维技术。例如并非对所有可能的三原子组合都进行均匀网格划分而是根据训练数据中角度的分布频率在常见角度区域加密网格。同时利用描述符本身的对称性如角度置换对称性来压缩表格。精度验证策略表格化必然引入误差。我们设计了一套系统的验证方案插值误差随机生成一大批未在训练集中的原子构型分别用原始GAP和tabGAP计算其能量和力直接对比两者的差异。我们要求平均绝对误差MAE在能量上小于1 meV/atom在力上小于50 meV/Å。这个误差相对于DFT本身的误差和热涨落来说是可以接受的。静态性质测试计算晶格常数、弹性常数C11, C12, C44、空位形成能、表面能、堆垛层错能等。将tabGAP的结果与原始GAP、DFT以及可靠的实验值进行对比。这是检验势函数“基本功”是否扎实的关键。动态稳定性测试在NPT系综下运行一段时间的分子动力学观察体系是否能够稳定在正确的晶体结构晶格常数是否漂移以及是否会出现非物理的结构相变。实操心得数据库的“代表性”比“数量”更重要。早期我们曾试图通过简单增加完美晶体构型的数量来提升精度效果甚微。后来发现关键在于构型的“多样性”。特别是那些在目标模拟中如塑性变形、辐照可能出现的、远离平衡态的构型。例如在构建用于辐照模拟的势函数时如果没有包含足够多的高能畸变构型那么在计算阈值位移能时势函数对高能碰撞的响应就会非常不可靠可能错误地高估或低估缺陷的产生效率。一个实用的技巧是用你计划使用的模拟方法如经典MD先快速预跑一些你想研究的物理过程从轨迹中采样关键构型再加入DFT数据库进行训练。这样训练出的势函数在目标场景下会格外稳健。3. 大规模分子动力学模拟的实现与挑战有了可靠的tabGAP势函数我们就像获得了一把锋利且轻便的新武器终于可以去挑战那些以前不敢想的大规模模拟了。我们选择了两类具有代表性的问题来展示其能力一是金属辐照损伤研究中的核心参数——阈值位移能二是直观体现材料力学性能的单轴压缩模拟并特别考察了预存缺陷的影响。3.1 模拟阈值位移能方向依赖性与计算实践阈值位移能Threshold Displacement Energy, TDE是衡量材料抗辐照能力的一个重要参数它定义为将一个晶格原子永久击出其格点位置所需的最小动能。这个值并非各向同性而是强烈依赖于撞击方向。传统实验测量TDE非常困难而模拟计算是主要手段。模拟设置模型构建我们为Al、Cu、Ni分别创建了边长超过20 nm的立方单晶块体原子数在百万量级。使用周期性边界条件。体系先在NPT系综下弛豫到零压状态。初级撞出原子PKA设置选择模型中心区域的一个原子作为PKA。我们系统性地选取了fcc晶体中具有代表性的晶向如100,110,111以及这些方向之间的多个中间方向。给PKA赋予一个初始动能E_pka方向沿所选晶向。模拟过程在NVE微正则系综下运行分子动力学。模拟时间通常在几个皮秒ps量级这足以让碰撞级联发生并初步稳定。我们使用一个非常小的时间步长如0.1 fs来精确捕捉高速碰撞过程。结果分析模拟结束后使用位错分析工具如OVITO中的Wigner-Seitz缺陷分析或CNA来识别最终的稳定缺陷。如果发现了一个稳定的空位-间隙原子对Frenkel Pair并且该缺陷在模拟结束时~10 ps后没有复合则认为这次撞击是成功的。通过不断调整E_pka并进行二分法搜索最终确定在该晶向下产生稳定Frenkel Pair的最小能量即为TDE。tabGAP带来的优势统计可靠性由于计算效率高我们可以对每个晶向进行多次重复模拟考虑不同的初始热振动速度从而获得TDE的统计分布而不是一个单一值。这对于理解辐照损伤的随机性至关重要。大体系避免有限尺寸效应使用百万原子模型可以确保碰撞级联产生的缺陷完全被包裹在模型内部远离周期性镜像从而更真实地模拟块体材料中的行为。小体系可能会因为能量和缺陷泄露到镜像中而严重低估TDE。我们的计算结果成功再现了这三种金属TDE的各向异性特征通常沿密排方向如110的TDE较低因为原子更容易沿着密排面被击出而沿松散方向如100的TDE较高。与文献中的DFT和经典势函数结果对比tabGAP的结果与高精度数据吻合得很好证明了其在处理高能非平衡事件上的可靠性。注意事项热振动与模拟时间的影响。这是一个容易忽略但至关重要的细节。在模拟TDE时必须给体系一个初始温度如300K即让原子具有热振动。这个热振动会显著影响TDE值因为热涨落可能“帮助”或“阻碍”原子的撞出。因此报告TDE时必须注明模拟温度。另外模拟时间的长短也需要谨慎选择。时间太短可能无法判断产生的缺陷对是否稳定时间太长则计算成本激增。我们通常的做法是先用一个较长时间如10 ps确定一个大致的TDE范围然后在该值附近用较短时间如2-3 ps进行密集搜索并辅以对最终构型的能量弛豫来确认缺陷的稳定性。3.2 百万原子单轴压缩模拟从完整晶体到含缺陷体系单轴压缩是测试材料力学性能最基础的模拟之一但将其做到百万原子尺度并且对比完整晶体与含缺陷晶体的差异能揭示许多介观尺度的机理。完整晶体的压缩 我们沿001方向对完美的Al、Cu、Ni单晶施加单轴压缩。应变率设置为一个较低的值如 1e8 /s以近似准静态加载。模拟在NPT系综下进行在加载方向控制压力在横向允许盒子弛豫。应力-应变曲线tabGAP计算出的应力-应变曲线在弹性阶段与理论值符合得很好。通过曲线初始斜率得到的杨氏模量与我们的静态计算以及文献值一致。这是对势函数弹性响应部分的基本验证。塑性屈服当应变达到一定程度后材料发生屈服。在完美的单晶中屈服通常始于位错的形核。我们通过位错分析工具清晰地观察到了位错环的形核和扩展过程。对于Ni我们特别对比了不同势函数包括经典的EAM势和MEAM势的结果发现不同势函数预测的屈服强度和位错形核机制存在显著差异。这凸显了势函数本身对材料塑性行为预测的巨大影响也说明了使用高精度势函数的重要性。含预存缺陷晶体的压缩 为了更贴近工程材料的实际情况材料总是含有缺陷的我们在压缩前在晶体中引入了预存的缺陷。我们主要尝试了两种空位团簇在晶体中随机移除一小团原子例如一个由几十个空位组成的孔洞。刃型位错通过插入或移除半原子面在晶体中构造一条直的刃型位错线。然后对含有这些缺陷的体系进行同样的单轴压缩模拟。缺陷对力学响应的影响结果非常直观。含有预存缺陷的体系其屈服强度普遍低于完美晶体。缺陷成为了应力集中点和位错形核的优先位置使得材料在更低的宏观应变下就进入了塑性阶段。应力-应变曲线在屈服点附近会出现更平滑的过渡而非完美晶体那样突然的应力跌落。缺陷的演化我们能够实时观察在加载过程中预存位错是如何运动、增殖并与新形核的位错发生相互作用的。对于空位团簇可以看到在应力作用下它们可能作为位错形核的核心或者被运动中的位错所吸收。大规模模拟的技术要点并行效率tabGAP势函数文件通常较大几百MB在并行计算时每个进程都需要读取势函数表格。我们采用MPI-IO进行并行文件读取避免了单个进程读取再广播造成的瓶颈。在LAMMPS中通过合理的域分解domain decomposition使得每个处理器核所负责的原子其邻居信息尽可能本地化减少了进程间通信量。输出与可视化百万原子规模的轨迹文件非常庞大。我们不会每一步都输出而是以较低的频率如每1000步输出轨迹并主要保存原子的坐标和应力信息。使用OVITO等可视化软件时需要先进行降采样或只渲染感兴趣的区域否则硬件无法承受。对于缺陷分析我们通常在服务器上先用脚本如OVITO的Python脚本接口批量处理轨迹提取出缺陷数量、类型、位置等统计信息再进行分析绘图。计算资源一次百万原子、纳秒尺度的模拟在拥有上百个CPU核心的集群上通常需要数天到一周的机时。tabGAP相比原始GAP将计算时间从可能的上百天缩短到了几天这才是其真正的价值所在。踩坑实录应变率的选择与系统尺寸的匹配。在早期测试中我们曾用过高的应变率1e10 /s去压缩一个只有几万原子的模型。结果发现屈服强度异常的高并且位错形核过程非常剧烈不符合物理实际。这是因为高应变率引入了惯性效应相当于给材料一个“冲击”加载而小尺寸模型则无法提供足够的空间让位错自然形核和扩展缺陷行为受限于周期性边界条件。后来我们调整了策略对于研究准静态力学响应应变率应尽可能低通常不高于1e9 /s同时模型尺寸要足够大至少要大于位错形核的特征尺寸通常是几十纳米。对于百万原子模型1e8 /s的应变率是一个比较合理的选择它能在可接受的计算时间内获得相对接近准静态的结果。4. 结果分析与不同势函数对比经过一系列严格的测试和模拟我们对Al、Cu、Ni的tabGAP势函数有了全面的认识。下面用表格来直观对比一下关键性能表1静态性质对比以Al为例性质实验值/参考值本工作 tabGAP原始 GAP经典 EAM 势 (如 Mishin)DFT (PBE)晶格常数 (Å)4.05 [实验]4.0494.0484.0324.04弹性常数 C11 (GPa)108 [实验]106107114105弹性常数 C12 (GPa)62 [实验]61626260弹性常数 C44 (GPa)28 [实验]29283129空位形成能 (eV)~0.67 [实验]0.690.680.700.71(111) 表面能 (J/m²)~1.21.181.191.151.21内禀层错能 (mJ/m²)~120125123135118从表格可以看出tabGAP在静态性质上几乎完全复现了原始GAP的高精度并且与DFT和实验数据吻合得非常好显著优于某些经典EAM势。这证明了表格化过程没有损失核心的精度。计算效率的量化提升 我们在一台标准计算节点2×Intel Xeon Gold 6248, 40核上进行了基准测试。模拟一个包含10万原子的体系运行1000步分子动力学。原始GAP耗时约5200 秒。tabGAP耗时约65 秒。加速比接近 80 倍。这个提升是颠覆性的。更重要的是随着体系原子数的增加tabGAP的线性缩放性更好而原始GAP的O(N^2)或O(N^3)复杂度会使其计算时间急剧上升。这使得百万原子乃至千万原子的模拟从“理论可行”变成了“实际可操作”。Ni的势函数差异带来的启示 在单轴压缩模拟中Ni的表现尤为有趣。我们对比了tabGAP、一个经典的EAM势和一个MEAM势。结果发现屈服强度tabGAP预测的屈服强度介于EAM和MEAM之间更接近一些基于第一性原理计算的高级修正势函数的结果。位错行为EAM势倾向于预测全位错perfect dislocation的滑移而MEAM和我们的tabGAP则观察到更多的不全位错partial dislocation和层错stacking fault的形成。这与Ni相对较高的层错能但仍可发生不全位错分解的物理事实更为相符。原因分析经典势函数在拟合时通常优先保证弹性常数、空位形成能等静态性质而对缺陷核心结构、层错能等涉及电子结构细微变化的性质捕捉不足。机器学习势函数直接从DFT数据中学习能够更自然地继承这些复杂的电子效应信息因此在预测复杂塑性行为时可能更具优势。这个对比深刻地说明在进行材料力学性能模拟特别是涉及缺陷演化的模拟时势函数的选择绝非小事。一个在简单性质上表现良好的势函数可能在复杂场景下给出误导性的结果。tabGAP提供了一条途径让我们能以可承受的计算成本使用接近DFT精度的势函数来探索这些复杂过程。5. 应用展望、局限与未来工作方向基于tabGAP的成功验证它的应用前景非常广阔。在我看来以下几个方向是马上可以开展且极具价值的多元合金与高熵合金现代先进材料往往是多组元的。传统上为多元合金开发经验势函数极其困难。而机器学习势函数特别是像tabGAP这样高效的版本非常适合处理多元体系。我们可以通过构建包含不同元素比例和化学有序度的DFT数据库训练出能够准确描述固溶强化、短程有序、偏聚等现象的势函数从而模拟高熵合金的力学行为、相稳定性等。辐照损伤的长时间尺度模拟辐照损伤研究不仅关心初始的级联碰撞皮秒尺度更关心缺陷空位、间隙原子、位错环等在长时间微秒到秒尺度下的扩散、聚集、演化以及与晶界、位错等微结构的相互作用。tabGAP的高效性使得我们可以利用并行计算模拟更大尺寸的模型和更长的物理时间结合加速动力学方法更真实地预测材料在辐照环境下的性能退化。非均匀材料与多尺度模拟tabGAP可以作为介观尺度模拟如位错动力学、相场法的“高精度输入”。例如可以用tabGAP模拟一个小体积元精确计算出不同构型下的广义层错能GSFE曲线、位错核心结构、晶界迁移率等参数然后将这些参数输入到上位尺度模型中实现从原子到微米尺度的更可靠跨尺度连接。当然tabGAP目前也有其局限性这也是我们未来需要努力的方向内存占用对于元素种类多、描述符复杂的体系势函数表格可能非常庞大GB级别。虽然计算快但加载和存储是个问题。未来需要发展更智能的稀疏表格压缩技术和动态加载机制。长程相互作用当前的GAP/tabGAP框架主要处理短程相互作用。对于离子体系或者需要精确描述长程库仑力的材料需要与Ewald求和等长程作用力计算方法相结合这增加了复杂性。主动学习与数据库扩展势函数的可靠性依赖于训练数据库的覆盖度。如何高效地发现当前势函数的“认知盲区”即预测不确定性高的构型并自动进行DFT计算来补充数据库主动学习是构建更强大、更通用势函数的关键。未来的工作流程应该是“模拟-发现不确定性-DFT计算-更新势函数”的闭环。与机器学习力场的融合除了GAP还有像DeePMD、Neural Network Potentials (NNP) 等其他优秀的机器学习势函数。它们各有优劣。一个有趣的方向是探索如何将tabGAP的表格化思想与这些不同的机器学习架构结合或者开发混合势函数在保证精度的前提下针对特定问题追求极致的速度。最后一点个人体会从事计算材料模拟工具和方法在快速迭代但核心的科学问题始终在那里。tabGAP这类技术本质上是给我们提供了更锐利、更顺手的“显微镜”和“实验台”让我们能以更低的成本、更高的置信度去窥探和操控材料的微观世界。然而再好的工具也离不开使用者的物理直觉和严谨设计。在构建数据库时多想一步在分析结果时多问一个为什么比单纯追求模拟的规模和速度更重要。这项工作的代码和势函数文件我们已经开源希望它能成为社区的一个有用工具推动更多关于材料在极端条件下行为的有趣研究。毕竟看清原子如何舞蹈才能理解材料为何坚强。
