1. 项目概述与核心问题石英或者说二氧化硅SiO₂是地壳中含量最丰富的矿物之一。从我们日常使用的玻璃、光纤到构成大陆地壳的主体岩石再到记录地球数十亿年撞击历史的“地质硬盘”石英无处不在。然而当这块看似普通的矿物遭遇极端环境——比如一颗小行星以每秒数十公里的速度撞击地球——它的命运会发生怎样的剧变这个问题不仅关乎地球的过去也影响着我们对材料在极端条件下行为的根本理解。长期以来地质学家在撞击坑周围的岩石中观察到石英内部出现一种被称为“平面变形特征”的微观结构这被认为是冲击事件的“指纹”。传统观点认为在高压下石英要么直接变成非晶态的玻璃即“击变玻璃”要么先转变为热力学稳定的高压相——斯石英。但实验观测却充满了矛盾有些实验看到了非晶化有些则声称发现了结晶。这种分歧背后是动态加载条件下材料在皮秒到纳秒时间尺度内经历的、远离平衡态的复杂原子重排过程。实验手段受限于时间与空间分辨率难以捕捉这一瞬息的原子“舞蹈”。这正是计算模拟大显身手的地方。分子动力学模拟可以追踪每一个原子的轨迹但传统模拟的精度受限于描述原子间相互作用的“势函数”。经典势函数往往在远离训练条件的高压、高应变下失准。近年来基于机器学习的原子间势函数横空出世它通过学习海量第一性原理计算数据能以接近量子力学计算的精度实现大规模、长时标的分子动力学模拟为破解这一难题提供了前所未有的工具。我们这项工作的核心就是利用这样一个新开发的、专为二氧化硅体系优化的机器学习势函数对α-石英在高达60吉帕GPa的冲击压力下的相变过程进行了一次“原子级高清直播”。我们想回答几个关键问题在冲击波的瞬间挤压下石英内部的原子究竟是如何运动的非晶化和结晶是先后发生还是同时竞争最终形成的结构是什么它又是如何形成的2. 研究方法与技术栈解析2.1 机器学习势函数从“经验公式”到“数据驱动”传统分子动力学模拟的基石是原子间势函数它本质上是一个描述原子间相互作用力与能量的数学公式。对于二氧化硅这样的复杂体系开发一个既精确又高效的势函数是巨大挑战。经验势如Tersoff、BKS计算快但往往无法准确描述键的断裂与形成以及在极端条件下的相变能垒。我们本次研究采用的是一种基于原子簇展开ACE框架的机器学习势函数。它的工作原理可以通俗地理解为一个极其复杂的“插值器”或“拟合器”数据准备首先利用第一性原理密度泛函理论DFT计算成千上万个不同原子构型包括晶体、非晶、缺陷、不同压力和温度下的结构的总能量、原子受力和应力张量。这些数据构成了“标准答案”数据集。模型训练机器学习模型此处为ACE学习这些数据。它不预设物理公式而是通过调整数百万个参数构建一个复杂的函数使得对于训练集内的任何原子排列模型预测的能量、受力与DFT结果尽可能接近。泛化预测训练好的模型可以用于预测它从未“见过”的新原子构型的性质和受力。关键在于我们验证了该势函数对多种已知和未知的二氧化硅高压相如斯石英、赛石英、rosiaite型等都具有出色的预测精度其能量-体积曲线与DFT结果高度吻合如原文图1a所示。这意味着用它来模拟高压相变我们对其结果有很高的置信度。注意使用机器学习势函数时务必进行严格的“外推”测试。即用训练集未包含的结构来验证其可靠性。我们的势函数在多种高压相上的良好表现是本研究结论可信度的基石。2.2 模拟“大炮”Hugoniostat方法与冲击模拟如何在计算机里模拟一次微型撞击我们采用了Hugoniostat方法。这不是简单地给系统施加一个固定压力而是通过调节系统的体积和温度使模拟体系的宏观状态压力、内能、体积满足冲击波前后的Hugoniot关系——这是描述材料在冲击波中状态变化的物理守恒方程。简单说它能在分子动力学模拟中复现出材料被冲击波扫过后的真实热力学状态。我们的模拟体系包含超过50万个原子冲击方向沿着α-石英的c轴峰值压力设定为56 GPa以匹配对比实验条件。模拟总时长达60纳秒。这个时间尺度对于原子运动来说已经足够长可以观察到明显的结构弛豫和相变过程。2.3 “眼睛”与“大脑”结构识别与分析方法模拟产生了海量的原子轨迹数据如何从中解读出“相”的信息这是另一个挑战。我们采用了两种互补的分析方法传统方法多面体模板匹配与X射线衍射计算多面体模板匹配通过分析每个原子周围近邻的几何排列来判断其局域结构属于哪种理想晶体类型如面心立方FCC、六方最密堆积HCP等。这主要用于分析氧亚晶格的排列。X射线衍射计算通过Debye公式直接计算模拟体系的XRD图谱与实验测得的图谱进行直接对比。这是验证模拟结果与实验是否一致的最直观手段。创新方法动态图卷积神经网络结构识别传统方法对高度无序或复杂有序的结构识别能力有限。为此我们训练了一个基于动态图卷积神经网络的分类器。这个AI模型能够“看懂”原子的三维局部环境并将其分类到25种预定义的结构类别中包括各种晶体相、非晶相甚至熔体。这相当于给模拟过程安装了一个“智能原子分类器”可以实时、定量地统计不同结构相的比例随时间演化其分辨能力远超传统方法。2.4 探索相变路径系统应变与固态NEB计算为了理解实验中观察到的另一种重要相变——向rosiaite型二氧化硅的转变我们进行了另一组“可控变形”模拟。我们不再模拟复杂的冲击过程而是对石英原胞进行系统性的、各向异性的应变沿c轴方向施加0%到40%的压缩应变。在垂直于c轴的a-b平面内施加0%到30%的压缩应变。 对于每一种应变组合我们都对体系进行能量最小化和热退火处理然后利用DG-CNN识别最终结构。这相当于绘制了一张“应变-相图”清晰地揭示了rosiaite相稳定存在的应变条件窗口。此外为了定量计算从石英到rosiaite相变的能垒我们采用了固态微动弹性带方法。SS-NEB可以找到连接初始相石英和最终相rosiaite的最小能量路径并计算出路径上的能量鞍点即能垒。这对于理解相变发生的难易程度至关重要。3. 核心发现冲击压缩下的“先破后立”3.1 动态过程非晶化先行再结晶紧随我们的模拟清晰地揭示了一个两步走的动态过程如原文图2所示快速非晶化 300 ps冲击波加载后在最初的300皮秒内α-石英的长程有序结构迅速瓦解。DG-CNN识别结果显示石英相的原子比例急剧下降至接近零而非晶相的原子比例飙升至近90%。原子瞬间“迷失”了它们在晶格中的位置。缓慢再结晶~1-60 ns随后在数十纳秒的时间尺度上系统并未停留在无序状态而是开始重新组织。一种被称为缺陷型NiAs结构的相开始出现并逐渐增长。与此同时非晶相的比例持续下降。到约60纳秒时样本几乎完全结晶非晶相几乎消失。这个“先非晶化后结晶”的序列调和了以往实验中看似矛盾的现象。在极短的时间尺度上观察你看到的是非晶化如果有足够长的观测时间你会看到结晶。实验室的冲击实验时间尺度纳秒级恰好落在这个转变序列之中。3.2 最终产物有序的“无序”结构——d-NiAs相那么最终结晶出来的是什么不是热力学最稳定的斯石英而是一种亚稳相缺陷型NiAs结构。它的结构特点是氧原子排列成规则的六方最密堆积网格而硅原子则随机地占据其中一半的八面体空隙。这是一种“有序的氧亚晶格”加上“无序的硅亚晶格”的混合态。更精细的分析发现硅原子的分布并非完全随机。在d-NiAs相的“海洋”中存在着许多微小的“岛屿”这些“岛屿”内的硅原子呈现短程或中程有序其排列方式类似于其他已知的高压亚稳相如NaTiF₄型、SnO₂型、P2₁/c型以及seifertite。我们的模拟XRD图谱与Tracy等人的实验图谱高度吻合原文图3a并且都缺失了斯石英和rosiaite相的特征峰却都有一个介于10°到13°2Θ的宽化峰。传统上这个宽峰被归因于部分非晶化但我们的模拟表明在几乎完全结晶的样本中这个宽峰依然存在。我们通过分解XRD图谱的贡献发现这个宽峰主要来源于硅-硅原子间的相互作用。正是硅亚晶格中这种局域、不完美的有序性导致了衍射峰的宽化。实操心得在分析复杂体系的XRD时不能仅凭一两个主峰就断定物相。宽化峰不一定代表非晶也可能是由纳米晶、严重的晶格畸变或亚晶格无序引起。结合局域结构分析如我们的DG-CNN至关重要。3.3 为何是d-NiAs而非斯石英这是一个动力学竞争的问题。在56 GPa下斯石英在热力学上确实更稳定。但从非晶态结晶需要完成两步1) 氧原子重排成HCP格子2) 硅原子在HCP格子的八面体空隙中有序排列形成斯石英的特定图案。 我们的模拟显示第一步氧原子重排在纳秒尺度内就能完成从而形成d-NiAs相。然而第二步硅原子的长程有序化则需要硅原子进行扩散和重排这是一个更慢的动力学过程。在有限的模拟和实验时间内系统“卡”在了d-NiAs这个亚稳态上。相比之下对非晶二氧化硅的冲击模拟发现在更高的冲击温度下两步过程都能在几纳秒内完成直接形成斯石英。这说明温度是加速硅原子扩散、促使向稳定相转变的关键因素。4. 应变状态的关键角色通往Rosiaite相的路径4.1 “应变-相图”的启示第二部分工作我们通过系统性的应变模拟揭示了rosiaite相形成的秘密。我们发现α-石英能否直接、无扩散地转变为rosiaite型二氧化硅强烈依赖于应变状态原文图4b纯c轴压缩如果只沿着c轴压缩即使应变很大25%石英也倾向于直接非晶化而不会形成rosiaite。复合压缩c轴 侧向当在c轴方向施加25%-40%压缩的同时在垂直于c轴的平面内再施加5%-20%的压缩石英就会稳定地转变为rosiaite相。这完美地解释了为何在静态金刚石对顶砧实验中样品受到一定程度的侧向约束观察到了rosiaite相而在纯单轴冲击实验中侧向自由则观察到了d-NiAs相和非晶化。4.2 能垒计算与相变机理SS-NEB计算定量地支持了这一结论原文图4c,d。在存在侧向压缩~10%的情况下当c轴方向压力超过约15 GPa时rosiaite相就变得比石英更稳定且两者之间的能垒很低约150 meV/atom在室温下很容易跨越。这符合实验观测到的转变压力。 反之如果没有侧向压缩rosiaite相要到约40 GPa才稳定且能垒路径上会出现一个能量更低的非晶态中间态使得体系“滑向”非晶化而非转变为rosiaite。机制解读Rosiaite结构可以看作是α-石英结构通过原子协同剪切位移直接转变而来原文图1b示意。这种剪切转变需要特定的应变张量来“催化”。纯c轴压缩无法提供这种剪切驱动力导致结构失稳而非晶化而合适的非静水压应力侧向压缩则能引导原子沿着一条明确的路径滑移完成无扩散的切变相变。5. 综合讨论与地质学意义5.1 时间尺度与应变状态的“双刃剑”效应本研究最核心的启示在于时间和应变状态是决定二氧化硅在高压下命运的两位“主宰”。时间尺度在冲击加载的纳秒-微秒尺度动力学占主导。系统可能来不及找到最稳定的状态而是被“冻结”在像d-NiAs这样的亚稳态或者完成第一步重排氧亚晶格有序化但来不及完成第二步硅原子完全有序化。应变状态非静水压应力即剪切应力可以彻底改变相变路径。它能够打开一些在静水压下关闭的无扩散相变通道如石英→rosiaite同时抑制另一些过程。这解释了自然界和实验室中观察到的二氧化硅高压相变的多样性。天然撞击中石英颗粒随机取向与周围矿物相互作用会产生极其复杂、局域化的应变场。这可能导致在同一块岩石中不同颗粒经历了不同的应变历史从而产生d-NiAs、rosiaite、非晶化乃至斯石英等多种产物共存的现象。5.2 对“平面变形特征”成因的新解释我们的工作为石英的“平面变形特征”这种经典冲击效应提供了更统一的微观机制图景快速加载阶段在冲击波阵面极高的应变率可能导致石英沿着特定晶面发生局域化的剪切失稳这种失稳可能直接诱发向rosiaite型结构的无扩散转变或者先导致非晶化。高压持续阶段在高压平台期如果条件温度、时间、局部应变合适非晶区域或rosiaite相可能作为前驱体进一步重结晶为d-NiAs相或极少量斯石英。卸载阶段在压力释放后这些亚稳的高压相如d-NiAs rosiaite在常压下不稳定会崩塌回非晶态最终在石英晶体中留下永久的非晶质条带即我们观察到的平面变形特征。5.3 机器学习模拟的价值与展望这项工作充分展示了机器学习势函数在极端条件材料科学中的威力。它架起了高精度量子计算与大规模分子动力学模拟之间的桥梁使得在保持第一性原理精度的前提下模拟数百万原子、数十纳秒的复杂动态过程成为可能。这不仅是计算能力的提升更是研究范式的转变——我们可以更自信地用模拟来解读甚至预测实验。未来的方向可以包括更复杂的矿物体系将方法扩展到其他造岩矿物如长石、橄榄石的冲击响应研究。多尺度耦合将原子模拟揭示的本构关系与宏观冲击波模拟耦合构建更真实的地球物理模型。主动学习与势函数发展在模拟中实时发现新的、未预期的原子构型并自动将其加入训练集迭代优化势函数形成闭环的研究流程。6. 常见问题与模拟实践要点6.1 模拟设置中的关键参数与选择体系大小我们使用了超过50万个原子的超大体系。对于冲击模拟体系必须足够大通常10 nm以容纳冲击波阵面、稀疏波以及可能形成的相变区域同时避免有限尺寸效应。一般建议至少在一个维度上大于20 nm。时间步长对于含硅氧的体系由于氧原子质量轻、振动频率高必须使用较小的时间步长如0.5 fs或1 fs以确保能量守恒和数值稳定。势函数验证在使用任何机器学习势函数前必须对其在目标压力-温度-结构区间进行严格测试。我们不仅对比了能量-体积曲线还测试了弹性常数、声子谱确保动力学稳定以及对非晶态的结构描述能力。控温控压方法Hugoniostat中的阻尼参数tdamppdamp需要谨慎选择。过小的阻尼会导致系统振荡过大的阻尼则会使过程偏离真实的等熵加载线。通常需要通过测试观察系统压力、温度是否平滑地达到目标Hugoniot状态。6.2 结构分析中的陷阱与技巧DG-CNN分类器的局限性它只能识别训练过的结构。如果出现了全新的、未训练的结构它会被错误地归类到某个已知类别。因此分类结果需要与传统的序参量如键角分布、配位数、径向分布函数以及XRD计算相互印证。区分“相”与“局域有序”DG-CNN给出的是每个原子的局域结构分类。一个被分类为“seifertite”的原子团可能只是一个具有类似seifertite短程有序的小区域并不代表形成了具有长程周期性的seifertite晶粒。需要结合团簇分析如基于距离的聚类算法来统计真正独立晶粒的数量和大小。XRD计算中的细节计算模拟体系的XRD时必须考虑有限尺寸效应和原子热振动。我们采用Debye公式计算它直接基于原子对距离自然地包含了所有效应。与实验对比时务必使用相同的X射线波长或波长分布和洛伦兹极化因子。6.3 结果的可重复性与扩展性随机性冲击模拟的初始热运动、以及d-NiAs相中硅原子的初始随机分布可能对具体演化路径有细微影响。重要的结论如先非晶化后结晶、最终形成d-NiAs相应在多次独立运行中保持一致。应变模拟的路径依赖第二部分中系统应变的研究其最终结构可能依赖于退火过程的细节升温速率、最高温度、冷却速率。我们采用的协议加热-退火-冷却是为了提供足够的动能帮助体系跨越能垒找到在给定应变下的低能构型。在实际应用中需要根据具体物理过程调整热处理的参数。通过这项研究我们不仅深入理解了α-石英这一特定材料在动态高压下的行为更展示了一套结合前沿机器学习势函数、大规模分子动力学模拟和智能结构分析的完整研究框架。这套方法对于研究其他复杂材料在极端条件下的非平衡相变具有广泛的借鉴意义。从原子模拟中解读出的故事正在帮助我们重新书写那些刻录在岩石中的、关于地球古老撞击历史的篇章。