机器学习原子间势与连续介质模型在柔性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和层间距离d（0.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万个（i=40，对应θ < 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进行高精度探索和标定，用连续模型进行快速扫描和物理理解。这种多尺度、多方法联动的思路，将是未来探索复杂材料体系，尤其是像扭转范德华材料这样充满惊喜的领域，不可或缺的研究范式。