1. 项目概述当机器学习遇见二维材料原子重构在二维材料的研究中一个核心的挑战是如何精确地描述和预测其原子结构在复杂环境下的行为。比如当我们把两层二维材料像拧麻花一样相对扭转一个角度就会形成一种被称为“莫尔超晶格”的长程周期性结构。这种结构会催生出许多单层或对齐双层材料中没有的新奇物理现象比如魔角石墨烯中的超导和关联绝缘态。理解这些现象的第一步就是搞清楚原子们在这个扭曲的“舞台”上究竟是如何排兵布阵的——这个过程就是原子结构弛豫或重构。传统上计算材料学的“金标准”是密度泛函理论DFT。它能从量子力学第一性原理出发高精度地计算原子间的相互作用和能量。但它的“阿喀琉斯之踵”是计算成本体系原子数一旦超过几百个计算量就会呈指数级增长变得几乎不可行。而一个典型的、具有小扭转角比如1度的莫尔超晶格其超胞可能包含数万个原子这远远超出了DFT的直接处理能力。这就引出了我们项目的核心如何突破计算尺度的限制去窥探大尺度莫尔体系中原子的精细重构我们选择了两种互补的“武器”机器学习势函数MLIP和连续介质模型。前者像是一个从DFT数据中“学成归来”的超级学徒能以接近师傅DFT的精度但快几个数量级的速度去处理成千上万个原子的弛豫问题。后者则像一位宏观力学家它不关心每一个原子的具体位置而是将材料视为连续的弹性薄片用一套基于弹性理论的方程来描述其整体变形计算成本极低。我们的“战场”是硒化铟InSe这种高度柔性的二维材料。它的杨氏模量远低于石墨烯或过渡金属硫族化合物TMDs这意味着它的原子层更容易被弯曲和拉伸。我们预感到这种高柔性会让InSe的扭曲双层在更宽的扭转角范围内发生剧烈的原子重构从而可能孕育出独特的电子性质。本文将详细拆解我们如何训练一个专用于InSe的MLIP如何构建并调整连续介质模型以及这两种方法在揭示InSe扭曲双层强原子重构现象上的“对决”与互补。无论你是计算材料学的新手还是对机器学习在科学计算中应用感兴趣的开发者这篇文章都将带你深入这个交叉领域的前沿。2. 核心方法论解析从第一性原理到跨尺度模拟要研究InSe扭曲双层的原子重构我们需要在精度和效率之间找到平衡。DFT提供了基准精度但尺度受限经典经验势函数效率高但往往难以准确描述二维材料层间复杂的范德华相互作用和化学键合变化。我们的策略是构建一个从原子尺度到介观尺度的多层次研究框架。2.1 密度泛函理论DFT计算一切的基础所有工作的起点是高质量的DFT计算它为后续的MLIP训练和连续介质模型参数化提供了“真理数据”。计算设置与参数选择我们使用Quantum ESPRESSO软件包进行所有DFT计算。软件包的选择并非随意Quantum ESPRESSO在平面波赝势框架下具有优秀的并行效率和广泛的泛函支持特别适合处理二维材料体系。交换关联泛函与范德华修正这是关键之一。对于InSe这类层状材料层间相互作用主要由范德华力主导。我们采用了OPTB88泛函它包含了非局域的范德华修正vdW-DF2类型。这个选择是基于大量测试的标准的广义梯度近似GGA泛函如PBE会严重低估层间结合能而单纯的vdW修正又可能过高估计。OPTB88在描述结合能曲线方面取得了较好的平衡这对于后续准确拟合层间相互作用势至关重要。赝势我们使用了SSSPStandard Solid-State Pseudopotentials库中的赝势。具体是In的PBE赝势来自PSlibrary和Se的PBE赝势来自GBRV数据集。使用经过验证的高效赝势能在保证精度的前提下显著降低平面波截断能。收敛参数平面波截断能动能截断设为80 Ry电荷密度截断设为720 Ry。这个设置远高于一般体系的收敛要求目的是为了确保作用在原子上的力和应力计算足够精确。因为MLIP训练不仅需要能量还需要准确的原子力和应力作为监督信号。k点网格我们为不同大小的体系采用了等效的k点密度。例如对于原胞单层和双层使用12×12×1的网格对于6×6的单层超胞使用2×2×1的网格。确保k点采样一致性能避免因采样差异引入的系统误差让MLIP学习到纯粹的原子构型与能量/力之间的关系而不是k点网格的 artifacts。注意在准备MLIP训练集时DFT计算的收敛标准必须格外严格。力的收敛阈值通常要设到每原子几个meV/Å的量级远高于通常结构弛豫的标准~25 meV/Å。否则训练数据中的噪声会严重影响MLIP的预测准确性和稳定性。2.2 机器学习势函数MLIP的构建让原子“学会”相互作用MLIP的核心思想是用一个灵活的机器学习模型通常是神经网络来拟合从DFT数据中学到的势能面。我们选择了MACEMessage Passing with Atomic Cluster Expansion框架它是一种等变神经网络能自然地满足物理系统的旋转、平移和置换对称性这对于保证势函数的物理正确性至关重要。2.2.1 训练数据集的设计覆盖所有可能场景训练一个稳健的MLIP数据集的质量和覆盖面比数量更重要。我们的训练集精心设计了多个模块旨在让模型学会InSe体系所有重要的原子相互作用模式层内相互作用原胞单层包含完全弛豫的结构以及施加了单轴、双轴和剪切应变高达±4%的25种构型。这教会模型InSe单层在均匀应变下的弹性响应。6×6单层超胞这是关键的一步。我们通过分子动力学MD模拟在50K到600K的多个温度下对原始单层以及施加了-2%到2%应变的单层进行采样。MD模拟使用了一个通用预训练模型MACE-mp0来驱动。这一步的目的是采样原子偏离平衡位置的所有可能热涨落构型。对于将要发生的原子重构层内原子必然经历复杂的非均匀位移这个数据集确保了MLIP能够处理这些“离平衡态不远”的原子环境。层间相互作用对齐的原胞双层将两个弛豫后的单层以平行P和反平行AP两种方式堆叠。然后系统地在层间平面内位移矢量r0的不可约部分进行网格采样6×6网格同时改变层间距d从0.82 nm到0.94 nm步长0.01 nm。这产生了572个构型完整地映射了结合能随堆叠方式和层间距变化的势能面是描述扭曲双层中可变堆叠区域能量的基础。连接微观与介观大角度扭曲双层我们加入了两个DFT尚可处理的扭曲角度21.79°和13.17°的刚性结构及其弛豫后的结构。这为数万个原子的大体系弛豫提供了一个“桥梁”让MLIP在训练阶段就接触到莫尔超晶格这种特殊的周期性边界条件。整个训练集共666个构型虽然数量不算庞大但每个构型都经过精心设计覆盖了从晶格振动、均匀应变到复杂堆叠的相空间。图1的示意图直观地展示了这些构型类别。2.2.2 模型训练与验证精度与泛化能力的考验我们使用MACE框架设置了两层神经网络并采用了一个相对较大的截断半径8 Å。通常的默认值是5 Å但对于InSe双层层间距约0.87 nm为了准确捕捉层间相互作用对堆叠的依赖必须保证截断半径大于层间距使模型能“看到”相邻层的原子。训练采用分阶段的策略前200个epoch能量、力和应力的损失函数权重比为1:10:1重点优化力的预测之后开启随机权重平均SWA并将能量权重提高到100以精细收敛总能量。最终模型在训练集和验证集上的能量均方根误差RMSE低至0.1 meV/原子力的RMSE约为7 meV/Å。力的误差主要来自高温MD采样的构型这是可以接受的因为这些构型本身偏离平衡较远。验证环节至关重要。图2的散点图显示MLIP预测的单层和双层能量与DFT计算结果高度吻合。图3则更令人信服地展示了MLIP完美复现了DFT计算的关键结果不同高对称堆叠方式如XX‘ MX’/XM‘ MM‘ 2H下结合能随层间距d变化的曲线。特别是它准确捕捉到了AP堆叠中MM‘和2H两种堆叠在结合能和最优层间距上的微小差异以及P堆叠中MX’/XM‘这对孪生堆叠的简并性。这些细节的成功复现是MLIP能够可靠用于研究扭曲双层中复杂堆叠竞争的前提。2.3 连续介质模型从原子细节到宏观变形连续介质模型采取了完全不同的思路。它不再追踪每个原子的坐标而是将每一层InSe视为一个连续的弹性薄膜。模型的能量泛函由两部分组成E_total ∫ [U_elastic(r) W_adhesion(d(r), r0(r))] d²r其中U_elastic(r)是弹性应变能密度。它由单层的弹性常数拉梅常数λ和剪切模量μ和局部应变张量u_ij(r)决定描述了原子层被拉伸或剪切所需要的能量代价。这些弹性常数直接从我们DFT计算单层应变-能量关系中获得。W_adhesion(d(r), r0(r))是层间粘附能密度。它描述了在某个位置r当两层之间的局部面内偏移为r0(r)、层间距为d(r)时单位面积的结合能。这个函数的核心来自对前述DFT层间结合能数据对齐双层在不同r0和d下的能量的数学拟合。模型的关键简化与核心方程面外变形能忽略模型假设弯曲单层所需的能量远小于面内应变能因此忽略了面外变形的能量代价。这是一个常见的近似但对于特征尺寸很小的褶皱如域壁和节点这个近似可能会失效。局域堆叠近似模型假设在扭曲双层的每一个位置r其层间距d(r)会自动调整到该处局部堆叠r0(r)所对应的最优值d_min(r0)。这样粘附能就只变成了局部面内偏移r0(r)的函数˜E(r0) W(d_min(r0), r0)。傅里叶展开与最小化将˜E(r0)和原子面内位移场u(r)都用莫尔超晶格的倒格矢进行傅里叶展开。通过变分法最小化总能量E_total可以求解出位移场u(r)的傅里叶系数从而得到重构后的原子位置近似为r0(r) ≈ θẑ × r 2u(r)。我们将DFT计算的结合能数据拟合为一个包含多个指数项和三角函数的解析表达式公式3参数如表I所示。图6显示这个拟合能很好地描述不同堆叠下结合能随层间距的变化趋势尽管在最小值附近的变化平坦区域其精度不如MLIP。两种方法的本质区别MLIP是自下而上的它从原子尺度相互作用出发通过数值弛豫得到每个原子的精确位置。连续介质模型是自上而下的它基于宏观弹性理论和平均化的粘附势通过求解偏微分方程或其傅里叶空间形式得到连续的位移场。前者精确但计算量大后者高效但做了物理近似。我们的研究正是要探究在InSe这种高柔性材料中这些近似在多大程度上仍然成立。3. InSe扭曲双层的强原子重构现象拥有了MLIP这把“利器”我们就可以对包含数万个原子的大尺度莫尔超胞进行全原子弛豫了。我们系统计算了从接近30°到接近0°P堆叠或60°AP堆叠的一系列扭转角。结果清晰地揭示由于InSe单层的高柔性原子重构现象在非常宽的扭转角范围内都极其显著。3.1 重构的物理图像能量博弈原子重构的本质是两种能量之间的竞争收益通过变形使双层中更多区域处于低能量强结合的原子堆叠构型如P堆叠中的MX‘/XM’ AP堆叠中的MM‘。成本使单层发生面内应变和面外弯曲所需要的弹性能。InSe的杨氏模量约为50 GPa远低于典型的TMDs150-350 GPa和石墨烯~1 TPa。这意味着让InSe层变形的“成本”很低。因此即使扭转角较大莫尔周期较短变形的“收益”也容易超过“成本”从而导致重构在比TMDs大得多的扭转角下就开始出现。3.2 P堆叠小扭转角的重构特征对于P堆叠扭转角接近0°弛豫后的结构呈现出清晰的图案图4上排面外褶皱层间距不再是均匀的。在能量有利的MX‘/XM’堆叠区域两层原子靠得更近结合更强而在能量不利的XX‘堆叠区域两层原子被推开层间距增大。这种褶皱在所有扭转角下都存在并随着扭转角减小莫尔周期增大而加剧。面内畴结构形成当扭转角小于约5°时清晰的三角形畴结构开始出现。这些三角形区域是MX‘/XM’堆叠的“能量洼地”它们被一个由域壁和节点构成的网络分隔开。域壁是连接两个能量简并的MX‘和XM’堆叠区域的过渡带其堆叠构型近似为“DW”Domain Wall。节点则是三个畴相遇的点通常是能量最高的XX‘堆叠。收敛行为当扭转角减小到约3°以下时畴、域壁和节点的形状和尺寸基本不再随角度变化达到了一个“收敛”的图案。这意味着重构的微观特征有一个固有的长度尺度~1-2 nm一旦莫尔周期远大于此图案就会自我重复。3.3 AP堆叠扭转角近60°的重构演变AP堆叠的行为更为有趣呈现出两个不同的阶段图4下排大角度阶段如59.18°重构图案看起来与P堆叠类似只是低能量区域变成了MM‘和2H堆叠取代了MX‘/XM’高能量节点仍是XX‘堆叠。此时主导重构的动力仍然是最小化高能量XX‘区域的面积。小角度阶段如60.24°, 60.08°当扭转角非常接近60°时图3中揭示的AP堆叠结合能曲线的不对称性开始显现。MM‘堆叠的能量略低于2H堆叠。因此在能量驱动下MM‘畴会逐渐“吞噬”2H畴的面积导致两者形状不再对称。在原子弛豫结果中2H畴会收缩成一个点状的节点而MM‘畴扩张。这种细微的能量差异能被MLIP精确捕捉并驱动了不对称的重构。3.4 层间距变化的定量分析图5展示了沿莫尔超胞一条路径的层间距变化。可以清晰地看到在能量有利的畴中心MX‘/XM‘或MM’层间距达到最小值~0.86-0.87 nm。在能量不利的节点XX‘中心层间距达到最大值。随着扭转角减小层间距的波动幅度褶皱幅度显著增大。例如在θ1.08°时层间距在0.86 nm到0.91 nm之间变化波动高达0.05 nm这对于层间电子耦合会产生重要影响。这些原子尺度的重构细节特别是域壁和节点的精确几何形状、尺寸以及层间距的局域变化是理解后续莫尔超晶格中电子能带平坦化、局域态形成等物理现象的结构基础。MLIP为我们提供了这幅前所未有的高分辨率“原子重构地图”。4. 机器学习势函数与连续介质模型的正面较量既然我们拥有了全原子尺度的“黄金标准”MLIP和高效的连续介质模型一个自然而然的问题就是连续介质模型在预测InSe扭曲双层重构时到底有多准它在什么情况下会失效本节我们将两种方法的结果进行逐项对比。4.1 面内畴结构的预测对比图7展示了连续介质模型预测的重构图案层间距分布图。与图4的MLIP结果对比我们可以发现对于P堆叠小扭转角连续介质模型的表现相当出色。它准确地再现了三角形MX‘/XM’畴的形状、域壁的走向和厚度以及XX‘节点的位置。这说明当重构的特征尺寸畴的大小远大于原子尺度时连续介质模型将材料视为弹性连续体的近似是有效的弹性能与粘附能竞争的核心物理被模型很好地抓住了。对于AP堆叠近60°连续介质模型出现了明显偏差。特别是在扭转角非常接近60°时如60.08°MLIP结果显示2H畴已经收缩为一个小节点而MM‘畴占据主导。然而连续介质模型预测的2H和MM‘畴仍然近乎对称仅表现出轻微的弯曲。模型未能充分捕捉到MM‘与2H堆叠之间那微小的能量差异所驱动的不对称重构。这是因为连续介质模型的粘附势˜E(r0)是从拟合的DFT数据中通过傅里叶变换得到的在拟合和变换过程中这些非常细微的能量差异可能被平滑掉或未能被低阶傅里叶分量充分描述。4.2 层间距预测的深度剖析图8将两种方法计算出的、沿着超胞相同路径的层间距剖面进行了直接叠加对比。差异更加显著和有趣最优层间距的偏差即使在能量最低的畴中心如MX‘连续介质模型预测的最优层间距也与MLIP结果有微小差别。这源于图6所示的根本原因连续介质模型对结合能曲线W(d, r0)的拟合在最小值附近曲线较平坦处的精度有限。MLIP则能更精确地复现DFT的细节。不过对于判断畴是否形成这个差异影响不大因为关键因素是不同堆叠之间的能量差而非绝对层间距。节点处层间距的严重高估这是最关键的发现。在XX‘堆叠的节点处连续介质模型预测的层间距~0.923 nm显著高于MLIP的结果~0.908 nm。为什么连续介质模型的假设它采用“局域堆叠近似”即认为每个点的层间距d(r)都自动等于该处局部堆叠r0(r)所对应的最优值d_min(r0)。对于XX‘堆叠d_min本来就较大。MLIP揭示的物理在实际的原子弛豫中节点尺寸很小半径约1 nm。要让原子层在如此短的距离内从畴中心的低间距弯曲到节点处的高间距会产生不可忽略的面外弯曲能。连续介质模型完全忽略了这部分能量代价。能量权衡全原子弛豫MLIP会做一个权衡与其付出巨大的弯曲能将节点处的层间距顶到理论最优值不如接受一个略低的层间距从而节省弯曲能。因此实际的层间距在节点处被“压制”了。MLIP由于直接计算原子间作用力自然包含了这种弯曲效应。大角度下的失效对于扭转角较大的情况如图8中θ21.79°连续介质模型预测的层间距变化曲线虚线与MLIP结果实线相差甚远。这是因为在大角度下莫尔周期短原子重构很弱层间距变化主要反映的是刚性莫尔图案中不同堆叠的快速交替。此时“局域堆叠近似”和基于弹性理论的位移场展开都不再适用。连续介质模型本质上是为了描述弛豫主导的体系对于刚性或弱弛豫的莫尔图案其预测能力很有限。4.3 方法论启示与适用范围界定通过这场“较量”我们可以得出以下清晰的结论连续介质模型的优势与适用场景高效快速计算成本极低可以轻松扫描大量扭转角。定性可靠对于柔性材料在小扭转角下的重构它能正确预测畴结构的基本形态如P堆叠为理解物理现象提供直观的图像。适用于特征尺度较大的变形当畴、壁等结构的尺寸远大于单层厚度即弯曲能可忽略时其预测较为准确。机器学习势函数的优势与不可替代性原子级精度能提供每个原子的精确坐标包括面外位移。捕捉复杂物理能自动包含所有能量项如面外弯曲、键角扭曲等被连续介质模型忽略的效应。揭示细微差异对能量势能面的细微特征如AP堆叠中MM‘与2H的微小能量差极其敏感能驱动出正确的对称性破缺重构。跨尺度能力从几个原子到数万个原子保持一致的精度是连接DFT与宏观尺度模拟的桥梁。实践建议在研究一种新的二维材料扭曲体系时可以先用连续介质模型进行快速扫描初步判断重构发生的临界扭转角以及畴结构的大致形态。当需要定量研究电子结构如计算莫尔能带、分析局域应变、或者研究特征尺寸与层厚相当的精细结构如节点、锐利域壁时必须使用MLIP或同等精度的全原子方法进行弛豫。连续介质模型在节点处对层间距的高估可能会显著影响该处局域电子态的预测。我们的工作表明对于InSe这类高柔性材料其重构特征尺寸节点半径~1 nm已经小到使弯曲能变得重要因此全原子方法的价值更加凸显。MLIP在此展现了其独特优势它以接近DFT的精度解决了传统方法无法处理的大尺度原子弛豫难题。5. 实操、挑战与未来展望5.1 构建与应用MLIP的实操要点基于本项目的经验如果你想为自己的研究体系训练一个MLIP以下是一些关键的实操建议和避坑指南训练集构建是成败关键广度优于数量与其用MD生成成千上万个相似构型不如精心设计几十个覆盖不同物理场景的构型。务必包含① 零应变及各种应变下的原胞② 有限温度下大超胞的MD快照采样热涨落③ 覆盖所有感兴趣堆叠方式和层间距的双层构型④ 如果可能包含一两个中等尺寸的目标体系如我们的扭曲双层的DFT弛豫结果。力的收敛至关重要用于训练MLIP的DFT计算必须将力的收敛阈值设得非常严格例如 1e-4 Ry/bohr。不准确的力数据是MLIP训练不收敛或预测不稳定的主要元凶。能量锚点在训练集中加入孤立原子的能量可以帮助模型确定能量的绝对参考点改善其泛化能力。模型选择与超参数调试等变架构对于材料体系强烈推荐使用像MACE、NequIP、Allegro这类等变神经网络架构。它们内置了物理对称性比普通神经网络需要更少的数据就能达到更好的精度和泛化能力。截断半径对于层状材料或涉及表面吸附的体系务必确保截断半径大于相互作用的特征长度。对于范德华层状材料截断半径必须大于层间距。验证与测试一定要严格区分训练集、验证集和测试集。验证集用于调参测试集用于最终评估模型在“未见过的”数据上的表现。测试集应包含训练集未覆盖的堆叠方式或应变模式。大体系弛豫技巧软硬件使用LAMMPS等支持MLIP的成熟分子动力学软件进行弛豫。对于数万原子的体系需要足够的GPU或CPU资源。弛豫策略可以先使用较弱的收敛标准如力 0.05 eV/Å进行快速预弛豫再使用严格标准如 0.001 eV/Å进行精细弛豫。有时结合FIRE算法和共轭梯度法能获得更好的效果。可视化检查弛豫结束后务必用VESTA、OVITO等工具可视化原子结构。检查是否有原子飞离、层间是否发生不合理穿透、重构图案是否对称美观。这是发现潜在问题如势函数不稳定、弛豫未收敛的最直接方法。5.2 当前工作的局限性与未来方向本研究为理解柔性二维材料扭曲双层重构提供了有力工具和新的见解但仍有可拓展的空间MLIP的泛化能力边界我们训练的MLIP专注于InSe单层和双层的基态附近构型。如果要模拟高温相变、缺陷迁移或化学反应需要在训练集中加入相应的高能量构型。一个MLIP通常只在其训练数据覆盖的相空间区域内可靠。动态过程与有限温度效应本文只研究了零温下的静态弛豫结构。实际材料处于有限温度下原子会热振动畴壁可能会波动甚至移动。利用训练好的MLIP进行分子动力学模拟可以研究重构结构的温度稳定性、畴壁动力学等丰富物理。与电子结构计算的耦合原子重构的最终目的是理解其对电子性质的影响。未来的工作可以将MLIP弛豫得到的原子结构直接输入DFT软件进行电子结构计算。更进一步的思路是开发同时学习势能和电子哈密顿量的机器学习模型实现从原子结构到电子能带的一站式预测。扩展到异质结与更多材料本文聚焦InSe同质双层。该方法可无缝推广到InSe与其他二维材料如石墨烯、hBN、TMDs形成的异质结研究其界面重构。也可以为其他高柔性二维材料如黑磷、GeSe等构建MLIP探索其扭曲体系的共性规律。5.3 对实验的启示我们的理论预测对实验观测有直接指导意义扫描探针显微镜SPM预测的~0.05 nm层间距起伏和清晰的三角形畴结构应该可以通过高分辨率的原子力显微镜AFM或扫描隧道显微镜STM进行验证。特别是节点处被“压制”的层间距是一个可以检验的关键预测。光学光谱原子重构会显著改变层间耦合和局部应变场从而影响激子电子-空穴对的能量和空间分布。预计在重构形成的畴网络边界域壁处会出现局域的激子态在荧光光谱上可能表现为额外的峰或展宽。电学输运畴壁和节点作为原子堆叠的“缺陷”网络可能会散射载流子影响材料的迁移率。同时重构导致的周期性应变场也可能产生赜磁场影响电子输运。通过将机器学习势函数与连续介质模型相结合我们不仅高效精准地揭示了高柔性InSe扭曲双层中丰富的原子重构物理更展示了一条通往更复杂二维材料体系多尺度模拟的可行路径。这项工作就像为材料科学家提供了一副兼具“显微镜”视野和“广角镜”效率的新眼镜让我们能更清晰地审视这些扭转的原子世界中隐藏的秩序与奥秘。