机器学习原子间势与连续介质模型在柔性InSe扭转双层原子重构研究中的应用

1. 项目概述当柔性二维材料遇上扭转角在二维材料的世界里一个简单的“扭转”操作往往能打开一扇通往新奇物理现象的大门。从魔角石墨烯中发现的超导和关联绝缘态到过渡金属硫族化合物TMDs中的莫尔激子扭转双层体系因其能产生长程周期性的莫尔超晶格并由此调制电子、光学和力学性质成为了凝聚态物理和材料科学的前沿热点。然而要真正理解这些“扭转”出来的新物理第一步也是最基础的一步是搞清楚原子们在这个新的、扭曲的舞台上究竟是如何排兵布阵的。这就是“原子重构”研究的核心。想象一下当你把两层原子晶体相对旋转一个小角度它们并不会像两块僵硬的玻璃板一样简单地叠加在一起。相反为了降低系统总能量原子们会“聪明地”发生微小的位移试图让两层原子在尽可能大的区域内形成能量上更有利的堆叠方式。这个过程就是原子重构。它直接决定了莫尔超晶格内的局域原子排列、应变分布和层间距离而这些微观结构细节正是后续一切奇特电子态如平带、范霍夫奇点和物理性质如铁电性、超导的舞台背景。那么如何精准地预测和描述这种原子重构呢对于像石墨烯这样相对“硬”的材料连续介质弹性理论模型已经能给出不错的近似。但对于一些“软”材料比如我们今天要讨论的主角——硒化铟InSe情况就复杂多了。InSe单层具有极高的柔性杨氏模量低至~50 GPa远低于石墨烯的~1 TPa这意味着它更容易发生形变原子重构的效应会在更大的扭转角范围内显现且重构形成的结构特征如畴壁、节点尺度更小。此时传统的连续介质模型可能会因为忽略了某些原子尺度的细节特别是面外褶皱的能量代价而“力不从心”。因此我们需要一把更精细的“尺子”来丈量这个原子世界。密度泛函理论DFT固然精确但其计算成本随原子数增长极快面对包含成千上万个原子的莫尔超晶格扭转角越小超胞越大DFT直接计算几乎是“不可能完成的任务”。近年来兴起的机器学习原子间势MLIP技术为我们提供了绝佳的解决方案。它通过学习DFT计算产生的高精度数据构建一个既快又准的“代理模型”使我们能够以接近DFT的精度对包含数万个原子的大体系进行完全原子尺度的弛豫计算。本文就将带你深入一项具体的研究工作如何运用机器学习原子间势MLIP与连续介质模型这两种“尺子”去探究高度柔性的InSe扭转双层中原子是如何进行大规模重构的。我们将拆解从势函数训练、大尺度弛豫到与连续模型对比的全过程并分享其中关键的实现细节、遇到的“坑”以及背后的物理考量。无论你是计算材料学的新手还是对二维材料扭转电子学感兴趣的探索者相信这篇详实的“实验笔记”都能给你带来启发。2. 方法论核心如何为柔性InSe打造一把高精度“原子尺”要研究InSe扭转双层的原子重构我们首先需要一把可靠的“尺子”来度量原子间的相互作用。这把尺子必须足够精确能捕捉InSe层内共价键/离子键和层间范德华力复杂的能量 landscape同时又要足够高效能处理数万原子的大体系。我们的策略是双管齐下构建一个基于第一性原理的机器学习原子间势MLIP作为金标准同时发展一个参数源自同一套DFT数据的连续介质模型作为快速近似和对比基准。2.1 基石第一性原理计算设置所有的起点都是高精度的密度泛函理论DFT计算它为MLIP提供训练数据也为连续模型提供参数。软件与泛函选择我们使用Quantum ESPRESSO软件包进行所有DFT计算。对于InSe这类层状材料范德华vdW相互作用至关重要。我们选择了OPTB88-vdW泛函它在描述层间结合方面表现出良好的平衡。赝势方面采用了SSSP库中经过验证的PBE赝势。注意泛函和赝势的选择需要谨慎测试。我们对比过几种vdW修正方法如DFT-D2, D3, vdW-DF发现OPTB88在结合能曲线和晶格常数上与实验值符合较好特别是对于层间距离和堆叠能差的预测。计算参数收敛性为确保数据质量我们对截断能波函数80 Ry电荷密度720 Ry和k点网格进行了严格的收敛性测试。对于不同大小的体系原胞、超胞、扭转超胞我们采用了等效k点密度的策略确保不同尺寸结构的k点采样一致性避免因采样差异引入的系统误差。2.2 锻造“原子尺”机器学习原子间势的训练实战训练一个优秀的MLIP其核心在于构建一个既能全面覆盖构型空间又不会过于庞大、导致训练困难的数据集。2.2.1 训练数据集的设计哲学我们的训练集不是随机采样的而是有明确物理目标的设计覆盖层内变形包含单层InSe原胞在不同单轴、双轴、剪切应变高达±4%下的能量和应力。更重要的是我们生成了6×6的单层超胞并通过分子动力学MD模拟在不同温度50K-600K和应变下采样获取了原子偏离平衡位置的各种构型。这确保了MLIP能学习到所有可能在扭转双层弛豫中出现的层内原子相对位移和面内应变。覆盖层间相互作用这是关键。我们计算了对齐双层原胞在不同面内平移矢量r0和层间距离d0.82-0.94 nm下的结合能。涵盖了平行P和反平行AP两种堆叠方式。这直接教会MLIP如何根据两层原子的相对位置和距离计算结合能。包含扭转体系本身我们将两个能用DFT直接计算的大扭转角21.79°和13.17°体系包括其刚性和弛豫后的构型也加入训练集。这起到了“锚定”作用让MLIP在目标应用场景附近有更准确的表现。2.2.2 模型架构与训练细节我们选择了MACE等变消息传递神经网络作为MLIP的架构。相比一些传统的高斯过程或简单神经网络MACE这类基于图神经网络的模型能天然满足物理系统的平移、旋转和镜面对称性等变性对于学习复杂的原子环境信息更有优势。几个关键的超参数和训练技巧截断半径设为8 Å。这是一个比默认值通常5 Å更大的选择。为什么因为InSe的层间相互作用范德华力作用范围相对较远且堆叠能差对原子环境敏感。更大的截断半径确保了模型能“看到”足够远的邻居从而精确捕捉层间相互作用的细节。损失函数与优化损失函数是能量、力和应力的加权平方和。训练分两阶段前200轮我们给力较大的权重能量:力:应力 1:10:1让模型快速学习力的变化之后启用随机权重平均SWA并调整权重为100:10:1专注于能量的精确收敛。优化器使用Adam的AMSGrad变体。验证与误差最终模型在测试集上的能量均方根误差RMSE低至0.1 meV/atom力的RMSE约为7 meV/Å。力的误差主要来自高温MD采样的构型这些构型中原子位移较大属于数据集的“困难样本”但这个误差水平对于结构弛豫来说已经足够精确。2.2.3 势函数验证不只是看误差曲线训练完成后不能只看测试集误差就万事大吉。我们进行了更物理的验证结合能曲线将MLIP预测的不同堆叠构型如XX‘, MX’, MM‘等的结合能随层间距d的变化曲线与DFT原始数据点直接对比。如图3所示两者吻合得非常好特别是能量最小值的位罝和深度这是决定重构驱动力堆叠能差的关键。弹性常数从MLIP计算的单层InSe应力-应变关系中提取出的杨氏模量和剪切模量与DFT直接计算的结果一致。这验证了模型描述材料力学性质的能力而力学性质直接决定了重构的“阻力”。经过这一系列步骤我们得到了一把为InSe量身定制的、高精度的“原子尺”——MLIP。它继承了DFT的精度又拥有了处理大体系的速度接下来就可以用它来探索扭转双层的广阔世界了。2.3 快速近似“尺”连续介质模型的构建与适配连续介质模型将原子层视为弹性薄膜其变形能用弹性理论描述而层间相互作用用一个依赖于局部面内偏移r0的粘附能密度W(r0)来近似。我们的模型主要基于Enaldiev等人为TMDs发展的框架但针对InSe的特性做了关键调整。模型核心公式系统的总能量是弹性能与粘附能之和E ∫ d²r [U(r) W(d(r), r0(r))]其中U(r)是弹性能密度取决于层的应变张量W是粘附能密度是局部层间距d和面内偏移r0的函数r0(r) ≈ θ ẑ × r 2u(r)包含了扭转和原子位移u(r)的贡献。InSe带来的挑战与修改粘附能拟合我们发现InSe的粘附能W(d, r0)随d和r0的变化行为与TMDs不同。特别是W中与r0反对称部分对应不同堆叠的能量差随d的衰减无法用单一指数函数很好地拟合。因此我们在拟合公式中引入了两个指数衰减项从而更准确地描述了这一依赖关系。参数获取弹性参数拉梅常数λ和剪切模量μ直接从我们用于MLIP训练的DFT单层应变计算中提取。粘附能面W(d, r0)的参数则通过对高密度网格24×24×1 in r0空间上的DFT计算数据进行傅里叶分析得到。关键近似该模型做了一个重要简化假设在任何位置层间距d都自动取该处局部r0所对应的最优值d_min(r0)并且忽略了面外弯曲的弹性能代价。这个近似对于特征尺度较大的重构可能成立但对于InSe这种柔性材料形成的、尺度很小的畴壁和节点就可能带来偏差。至此我们拥有了两把“尺子”一把是原子尺度的精密“游标卡尺”MLIP另一把是宏观快速的“卷尺”连续模型。接下来就用它们去测量扭转InSe的世界。3. 原子重构图景MLIP揭示的柔性InSe扭转双层结构演化有了训练好的MLIP我们便可以放手对一系列不同扭转角从近30°到近0°的P型和AP型InSe扭转双层超大超胞进行完全原子弛豫。这些超胞的原子数最多达到近4万个i40对应θ 1°。计算本身使用了ASE和LAMMPS等工具关键在于收敛标准我们要求所有原子上的力分量最大值小于0.25 meV/Å以确保结构真正弛豫到了极小点。3.1 重构的物理驱动力一场能量博弈原子重构的本质是两种能量之间的竞争收益通过让两层原子在更大面积上形成低能量堆叠如P堆叠下的MX‘/XM’或AP堆叠下的MM‘来降低系统的粘附能。成本使单层材料发生弹性形变拉伸/剪切所需要付出的弹性能。对于InSe其极低的杨氏模量~50 GPa是决定性的。这意味着使其形变的“成本”很低。因此即使堆叠能差“收益”与TMDs相似InSe也愿意在更大的扭转角下就开始“支付”形变成本以换取粘附能收益。我们预测重构现象会在比TMDs和石墨烯更大的角度范围内出现。3.2 P型堆叠θ 近 0°的重构景观弛豫结果完全证实了我们的预测。图4上排清晰地展示了这一演化过程大角度如13.17°虽然仍是连续的莫尔条纹但已经能观察到明显的面外褶皱。层间距不再均匀在高能量堆叠区域如XX‘层间距较大在低能量堆叠区域如MX’层间距较小。中等角度~6°三角形畴的雏形开始出现。低能量的MX‘/XM’堆叠区域开始扩大并形成清晰的三角域被高能量的XX‘堆叠线畴壁所分隔。此时重构的网格已经肉眼可见。小角度 3°形成完全发育的、规则的三角形畴网络。MX‘/XM’三角域被狭窄的XX‘堆叠畴壁隔开并在三个畴壁的交汇处形成XX‘堆叠的节点。这个结构在角度进一步减小时基本保持自相似只是摩尔周期变大。一个关键细节我们通过计算两层中In原子亚层平均面的垂直距离来可视化层间距的分布。图5中的线扫描显示从畴中心MX‘到畴壁XX‘层间距的变化幅度可达0.05 nm以上。这种显著的面外起伏是柔性材料和强层间相互作用共同作用的结果。3.3 AP型堆叠θ 近 60°的重构对称性破缺的故事AP型堆叠的重构更为有趣它展现出两个不同的阶段图4下排第一阶段角度远离60°重构图案与P型堆叠非常相似只是低能量堆叠区域从MX‘/XM’换成了MM‘和2H。此时主导重构的动力仍然是最小化高能量XX‘堆叠区域的面积MM‘和2H的能量差异很小因此它们形成的三角域大小和形状几乎对称。第二阶段角度接近60°当扭转角非常接近60°时堆叠能量的细微不对称性开始起主导作用。DFT计算显示在AP堆叠中MM‘堆叠的能量略低于2H堆叠这与TMDs中常见的2H最稳定不同。于是在重构中能量更低的MM‘域会逐渐“侵蚀”2H域导致两者形状不再对称MM‘域扩大2H域收缩。在最小的角度下如59.18°2H域甚至收缩成一个节点。实操心得在分析AP型堆叠时要特别注意初始堆叠的设定。我们的“反平行”是指将一层旋转180°后再相对扭转这会导致莫尔图案和重构对称性与P型堆叠不同。在构建初始原子结构时务必检查清楚两种堆叠方式的定义否则后续分析可能完全错误。4. 双尺对比MLIP与连续介质模型的较量与启示现在让我们拿出另一把“尺子”——连续介质模型看看它量出的世界与MLIP的原子尺度“金标准”有何异同。图7和图8展示了直接的对比。4.1 连续模型的成功与局限成功之处对于P型堆叠畴的形状与畴壁厚度对于P型堆叠连续模型非常好地再现了MLIP预测的三角形畴形状以及畴壁的宽度。这说明在描述面内原子位移和应变场方面基于弹性理论的连续模型抓住了主要物理即系统通过面内形变来优化堆叠。暴露的局限层间距预测的偏差这是最显著的差异。如图8所示连续模型预测的层间距变化曲线在低能量堆叠区MX‘的最小值、高能量节点XX‘的最大值都与MLIP结果有出入。这主要源于两个近似粘附能拟合误差连续模型对W(d, r0)的拟合在能量最小值附近较为平缓导致预测的最优层间距d_min(r0)本身就有微小偏差。忽略了面外弯曲能这是最关键的物理缺失。连续模型假设每个局部区域都独立地采取其最优层间距d_min(r0)完全忽略了将原子层弯曲成起伏形状所需要的弹性能。对于InSe这种柔性材料在节点和畴壁这些特征尺度很小~1 nm的区域面外曲率很大弯曲能不可忽略。MLIP的原子计算自然包含了这部分能量因此它“认为”将节点处的层间距抬升到像连续模型预测的那么高0.923 nm是“不划算”的最终平衡态下的节点层间距0.908 nm会更低。AP型堆叠畴形状的失效对于接近60°的AP型堆叠连续模型未能成功再现MLIP所揭示的MM‘域与2H域的不对称性。在连续模型的预测中即使角度很接近60°两个域仍然几乎对称图7下排。这说明当两种堆叠的能量差非常小~meV量级时连续模型对能量景观的细微差别不够敏感或者其简化处理如傅里叶展开的截断抹平了这种微妙的不对称性。而全原子MLIP计算则能精确捕捉这一细节。4.2 给我们的启示方法的选择取决于问题这次对比给了我们一个清晰的路线图何时用连续模型当你需要快速扫描大量扭转角、进行趋势性分析或者关注的主要是面内重构图案、应变场和摩尔势时连续模型是一个极其高效且物理图像清晰的工具。它对计算资源的需求极低可以轻松处理任意大的摩尔周期。何时必须用原子尺度方法如MLIP当研究的材料像InSe一样柔性或者当你关心层间距的精确分布、原子级粗糙的界面、以及由极其微小的能量差驱动的对称性破缺时全原子方法就不可或缺。特别是对于后续计算电子结构如能带、态密度局域层间距和原子位置的微小变化都可能对结果产生重大影响。避坑指南在采用文献中已有的连续模型代码时切忌直接套用。务必检查其粘附能拟合公式是否适用于你的材料。我们的经验是至少要用你的DFT数据重新拟合关键参数并像本文一样与少量关键角度的原子尺度计算结果进行比对验证以评估该连续模型在你的体系中的可靠性边界。5. 技术实现深潜从数据准备到大规模弛豫的实操要点理论很美但实现起来细节决定成败。这部分分享一些在具体操作中积累的经验。5.1 构建高质量MLIP训练集的技巧训练集的“质量”比“数量”更重要。我们的数据集只有666个构型但针对性极强。层内采样用MD采样热扰动和应变下的构型比随机扰动更物理、更高效。关键是温度范围和应变范围要覆盖所有可能出现的原子位移。我们的温度上限600K低于InSe熔点应变范围±2%则宽于重构中预期的内部非均匀应变通常1%。层间采样对齐双层的采样网格要足够密。我们采用了6×6的面内平移网格和0.01 nm的层间距步长。对于P堆叠利用其对称性可以减少计算量。“锚点”构型一定要包含目标体系本身哪怕只有一两个小角度的扭转双层。这能极大地提升模型在目标区域的预测精度。我们加入了21.79°和13.17°的扭转体系。使用预训练模型加速在运行MD生成热扰动构型时我们使用了通用的MACE-mp0基础模型来驱动MD。虽然它最终不能精确复现InSe的堆叠能和振动性质但用于快速生成合理的原子运动轨迹是足够的这节省了大量DFT计算时间。5.2 大尺度扭转超胞弛豫的实用策略弛豫一个包含近4万个原子的体系需要一些计算工程上的考量。初始结构构建采用“公度超胞”方法通过公式生成具有精确扭转角θ和摩尔周期L_M的超胞。确保初始两层是刚性扭转的没有原子重叠。弛豫算法选择对于原子数少于4000的体系扭转角2.5°使用ASE内置的优化器如FIRE或BFGS很方便。对于更大的体系我们将训练好的MACE势函数编译成LAMMPS可用的格式利用LAMMPS强大的并行能力和内存管理进行弛豫。LAMMPS的min_style cg或min_style hftn是不错的选择。收敛标准与监控力收敛标准如0.25 meV/Å需要足够严格但也要平衡计算成本。建议同时监控系统总能量和最大位移的变化。可以分阶段弛豫先用较松的标准如1.0 meV/Å快速接近极小点再用严格标准精细优化。可视化与分析弛豫完成后需要将原子位置数据转化为物理图像。我们编写了脚本将每个原胞位置的面内偏移r0(r)和局域层间距d(r)计算并插值成二维彩图如图4。分析畴壁宽度、节点尺寸、面外起伏幅度等关键参数随扭转角的变化趋势。5.3 连续模型求解的数值细节求解连续模型意味着在傅里叶空间中最小化总能量泛函公式4。傅里叶展开截断位移场u(r)的傅里叶级数不能无限展开。我们需要截断到某一阶摩尔倒格矢g_n。通常取到|g|小于某个截断值就足够了。需要测试收敛性不断增加截断阶数直到重构图案和能量不再发生显著变化。能量最小化这是一个关于傅里叶系数u_{g_n}的无约束优化问题。可以使用共轭梯度法等标准算法。由于能量泛函是u_{g_n}的二次型弹性能加上一个周期函数粘附能优化通常能较快收敛。从u(r)到物理量得到u(r)后通过r0(r) ≈ θ ẑ × r 2u(r)计算局域偏移再通过查表或插值得到d_min(r0)和W最终可以画出层间距分布图和能量分布图。6. 总结与展望柔性扭转范德华材料的研究范式通过这项结合机器学习原子间势和连续介质模型的研究我们清晰地揭示了高度柔性的InSe扭转双层中原子重构的独特行为由于其极低的形变能成本显著的重构包括面外褶皱和面内畴形成在很宽的扭转角范围大到6°甚至更高内就已出现。这为在相对“宽松”的转角制备条件下观测和调控InSe莫尔超晶格提供了可能。从方法论角度看这项工作展示了MLIP在二维材料扭转体系研究中的强大能力。它成功跨越了从DFT的精确但昂贵到连续模型的快速但近似的鸿沟为我们提供了一种既能处理数万原子大体系、又能保持量子力学精度的新工具。尤其对于InSe这类柔性、且电子结构对堆叠和层间距极度敏感的材料原子尺度的细节至关重要MLIP几乎是目前唯一可行的全原子弛豫方案。未来的工作可以沿着几个方向展开电子性质计算将MLIP弛豫得到的精确原子结构作为输入进行大尺度电子结构计算例如通过紧束缚模型或机器学习力场结合电子紧束缚方法预测InSe莫尔平带、范霍夫奇点的位置以及可能的关联电子态。探索更复杂的异质结将方法推广到InSe与其他二维材料如石墨烯、hBN、不同TMDs组成的扭转异质双层。此时层间相互作用更加复杂对训练集设计和MLIP泛化能力提出更高要求。动力学与有限温度效应利用训练好的MLIP进行分子动力学模拟研究重构结构在不同温度下的稳定性、畴壁的动态运动等。连续模型的改进是否可以发展一个包含面外弯曲能项的“高阶”连续模型或者将MLIP预测的力场作为“教师”来训练一个更准确的、基于神经网络的“粗粒化”势函数以兼顾效率与精度最后分享一个深刻的体会在计算材料学中没有一种方法是万能的。MLIP、连续模型、第一性原理、紧束缚模型……它们是一套相辅相成的“组合工具”。研究的艺术在于根据具体问题的尺度、精度要求和计算资源灵活选择和搭配这些工具。本研究正是一个范例用DFT奠定基础用MLIP进行高精度探索和标定用连续模型进行快速扫描和物理理解。这种多尺度、多方法联动的思路将是未来探索复杂材料体系尤其是像扭转范德华材料这样充满惊喜的领域不可或缺的研究范式。