面向固体原子与连续介质同步模拟的快速高精度准原子方法

《Computer Physics Communications》：Fast and accurate quasi-atom method for simultaneous atomistic and continuum simulation of solids

　　

研究人员报告了一种新颖的混合方法，用于在关键区域（如接触表面、裂纹区域等）进行固体的原子尺度同步模拟，同时对其他部分进行连续介质建模。连续介质通过不同尺寸的准原子(quasi-atoms)进行处理，这些准原子构成了复合介质。准原子之间的相互作用势参数经过优化，以使复合介质的弹性性质与原子介质的弹性性质相匹配。该优化方法在概念上与在线机器学习(ML)方法一致，使其在计算上非常高效。这种方法允许直接应用标准的分子动力学(MD)软件包，并辅以基于ML的优化器。新方法被应用于模拟具有简单成对Lennard-Jones(LJ)势以及描述共价键的多体Tersoff势的系统。利用LAMMPS软件，研究人员模拟了不同尺寸粒子的碰撞。将新方法获得的模拟结果与全原子模拟进行比较，证明了其准确性、有效性和在计算速度上的压倒性优势。此外，研究人员将该方法与其他混合方法进行了比较，特别是与最接近的方法——原子到连续介质(AtC)方法进行了比较。证明了该方法在计算速度和实施便利性方面的显著优越性。最后，研究人员讨论了该方法在模拟其他现象方面的可能扩展。

论文解读：《面向固体原子与连续介质同步模拟的快速高精度准原子方法》

研究背景与意义

随着计算能力的指数级增长，分子动力学(MD)已成为研究原子尺度材料不可或缺的工具。然而，对于大规模系统的全原子模拟在计算上仍然难以实现。为了克服这一限制，介观和混合方法应运而生，允许在不同区域以不同分辨率进行处理。许多物理过程（如固体中的裂纹扩展、纳米颗粒聚集、烧结、自组装、枝晶生长或相分离动力学）迫切需要同时具备原子描述和连续介质描述，以捕捉局部化学性质和长程形貌。传统的连续介质或纯粹的粗粒化方法往往在此类情况下失效。尽管现有的多尺度策略（如耗散粒子动力学(DPD)、光滑粒子流体动力学(SPH)、自适应分辨率方案(AdResS)等）通常用于流体或软物质，但可能无法忠实再现晶体固体所需的特性。此外，现有的混合方法（如粗粒化MD、异质多尺度方法(HMM)或准连续(QC)方法）大多针对液体或气体，或者存在求解器开销大、设置复杂以及对动态系统适用性有限等问题。为此，研究人员开发了一种新颖的快速高精度混合MD-介观方法，旨在解决晶体固体（如弹性、断裂模拟等）的动态系统描述难题。该研究由Artem Chuprov、Egor Nuzhin、Alexey Tsukanov及Nikolay Brilliantov共同完成，发表于《Computer Physics Communications》。

关键技术与方法

为实现上述目标，研究人员采用了以下几个主要关键技术方法：首先，提出了“准原子”概念，将其作为连续区域的动态离散节点，而非真实原子的集合，其相互作用势旨在重现固体的宏观弹性响应。其次，构建了包含真实原子和逐级放大尺寸准原子的混合晶体模型，准原子尺寸按因子2递增，相应的晶格常数和质量也随之增加。第三，采用基于径向基函数(RBF)的代理优化(surrogate optimization)方法，通过最小化目标弹性张量与模拟弹性张量之间的平均相对误差来校准准原子间的相互作用势参数（如ε_k和σ_k）。第四，利用Python-LAMMPS桥接技术，将高性能C++ MD引擎与灵活的基于Python的优化器相结合，实现了参数优化的自动化与高效运行。

研究结果

混合模型描述

在多尺度建模中，系统在关键区域采用原子描述，而在远离关键区域的部分，准原子尺寸可逐渐增大至所需值。对于Lennard-Jones(LJ)势，每个原子和准原子类型有两个参数(σ_i, ε_i)，交叉相互作用参数通过Lorentz-Berthelot规则确定。准原子间的附加谐波相互作用也被视为优化参数的一部分，以提高优化精度。

仿真结果

优化得到的有效势参数

优化程序在30-40次迭代后迅速收敛，目标参数的相对误差小于0.94%。这证实了该方法在数值上的高效性。优化程序自动调整准原子势参数以精确匹配目标弹性模量，且特定原子势的选择不影响准原子模型的转移性。

粒子碰撞：原子与多尺度建模对比

通过比较小颗粒（半径75 ?）和大颗粒（半径0.1 μm）的碰撞，验证了多尺度方法的准确性。在纯排斥和吸引两种相互作用机制下，复合粒子与原子的碰撞动力学均表现出极好的一致性。对于微米级大颗粒，模拟结果与赫兹(Hertz)非粘附接触理论和约翰逊-肯德尔-罗伯茨(JKR)粘附接触理论进行了比较，在机械接触区域（重叠量ξ>0）吻合良好。然而，在缺乏机械接触的负重叠区域（ξ<0），观察到与JKR理论的显著偏差，这归因于原子离散性导致的非接触原子间范德华力未被JKR理论考虑。此外，该方法在Tersoff势（用于硅颗粒）下同样表现良好，证明了其在超越成对势场景下的普适性。

多尺度方法的计算效率

随着准原子类型的增加和系统尺寸的增大，计算加速效果急剧提升。例如，对于半径为30 ?的颗粒，复合颗粒模拟比全原子模拟快十倍；对于半径为40 ?的颗粒，加速效果更为显著。在粒子半径35 ?时，使用四种原子/准原子类型的模拟比使用两种类型快约四倍。

本方法与其它方法的比较

通过与原子到连续介质(AtC)方法在平衡铜系统（T=300 K）中的计算时间对比，发现即使在多尺度方法仅使用两种类型（N=2）且分辨率显著高于AtC的情况下，其计算速度仍明显更快。随着N的增加（如N=4），这种优势将更加明显。这表明所提出的混合方案在计算速度、部署简便性以及目标指标的精确度上均优于现有的多尺度方法。

结论总结

研究人员开发了一种新颖的混合多尺度模拟方法，能够准确描述原子尺度的相互作用和连续介质动力学。该方法的核心在于在关键区域使用全原子模拟，而在系统的其余部分使用由尺寸逐级增大的准原子构成的连续介质描述。准原子的相互作用势参数通过代理优化程序进行校准，以确保有效介质的弹性性质与真实原子介质一致。该方法通过Python-LAMMPS桥接技术实现，能够直接利用标准MD软件包。通过粒子碰撞模拟验证了该方法的高精度和计算效率优势。此外，该方法可推广应用于裂纹扩展、压痕等更广泛的固体力学多尺度问题，并有望在未来扩展到热传导、动态断裂扩展等非平衡过程，甚至通过强化学习进一步提升优化效率。