机器学习势函数在氧化镓多晶型相变模拟中的应用与验证

1. 项目概述与核心挑战

氧化镓（Ga2O3）作为下一代宽禁带半导体的明星材料，这几年在功率电子和深紫外光电器件领域的热度一直居高不下。它的优势很明显：超宽的禁带宽度（4.8-5.3 eV）、极高的临界击穿电场（约8 MV/cm），这些特性让它非常适合用来做高压、高温、高效率的功率器件。但搞材料研究的人都知道，Ga2O3有个让人又爱又恨的特点：它是个“多面手”，存在至少五种已知的稳定多晶型结构，包括最常见的单斜β相，以及κ、α、δ、γ等亚稳相。更复杂的是，这些相之间能量差很小，在生长、退火或受到辐照、应力时很容易发生相变。你想精确控制薄膜的相结构来获得特定的电学或光学性能？对不起，传统实验手段像是在黑箱里摸索，成本高、周期长，而且很多微观动力学过程根本“看”不到。

这就是计算材料学大显身手的地方了。我们想用分子动力学（MD）模拟来“观看”原子在皮秒到纳秒时间尺度上的舞蹈，理解相变是怎么一步步发生的。但MD模拟的基石——描述原子间相互作用力的势函数（也叫力场）——在这里成了最大的瓶颈。传统的经验势（比如经典的Buckingham势、Tersoff势）对于Ga2O3这种离子性、配位复杂（Ga有4配位和6配位）的多元氧化物体系，精度往往不够，预测的晶格常数、弹性性质、相稳定性顺序可能和第一性原理计算（如密度泛函理论，DFT）结果相差甚远。而直接用DFT做MD（即从头算分子动力学，AIMD）呢？精度是够了，但算力成本是天文数字。一个包含几百个原子的体系跑几十皮秒，可能就需要耗费超级计算机数周甚至数月的时间。对于研究涉及成千上万个原子、时间尺度在纳秒以上的相变、缺陷演化或辐照损伤过程，DFT-MD基本是不可行的。

所以，问题的核心就变成了：我们能否找到一种方法，既能拥有接近DFT的精度，又能获得像经典力场那样的计算速度，从而对Ga2O3这种复杂体系进行大规模、长时间的原子尺度模拟？机器学习势函数（MLIP）的出现，让这个“鱼与熊掌兼得”的梦想照进了现实。它本质上是一个用机器学习模型（如神经网络、高斯过程）拟合从DFT计算中得到的大量（原子构型，能量/力）数据对的“超级插值器”。一旦训练完成，这个模型在预测新构型的原子间作用力时，速度可以比DFT快几个数量级，同时精度损失很小。

我们这次的工作，就是针对Ga2O3体系，开发并系统验证了两套基于高斯近似势（GAP）框架的机器学习势函数：高精度的soapGAP和追求极致速度的tabGAP。我们的目标不仅仅是复现某个单一相的性质，而是要打造一个“通用型”势函数，能够同时精确描述Ga2O3的所有五种主要多晶型（β, κ, α, δ, γ），并且能合理地外推到非晶、液态甚至非化学计量比等无序结构。最终，我们利用这个工具，首次在原子尺度上揭示了Ga2O3液固相变过程中一个此前未被实验观测到的、氧和镓亚晶格迁移率显著不同的三阶段动力学机制。

2. 势函数构建：从数据到模型的双轨策略

构建一个靠谱的机器学习势函数，七分靠数据，三分靠训练。数据是地基，模型是建筑，地基不牢，再漂亮的模型也会塌。对于Ga2O3这种多晶型体系，构建训练数据库的挑战在于，它必须足够“广”且足够“精”。

2.1 训练数据库的精心设计与构建

我们的数据库包含了1630个原子构型，总计108,411个原子环境。这些构型不是随意生成的，而是经过了精心设计和分类，旨在全面覆盖Ga2O3可能存在的相空间。我把它们分为四大类，这就像为机器学习模型准备了一份从“标准答案”到“开放性问题”的完整题库：

1. 核心晶体相（高对称性构型）：这是数据库的“压舱石”。我们不仅包含了所有五种实验已确认的多晶型（β, κ, α, δ, γ）的完美晶体结构，还额外加入了几个理论上预测的亚稳相结构。对于这些晶体，我们通过施加均匀的压缩和拉伸应变（例如从-3%到+3%），并弛豫内部原子位置，生成了它们在应变下的局部能量极小点构型。这一步至关重要，因为它让模型学到了晶体的弹性响应，而不仅仅是零压下的平衡结构。

2. 非化学计量比相与纯镓相：材料的性质不仅取决于完美晶体，缺陷（如氧空位、间隙原子）往往扮演关键角色。为了确保势函数在描述缺陷、表面或非理想化学配比区域时物理正确，我们引入了非化学计量比的GaOx（x≠1.5）构型，甚至包括纯金属Ga的相。这部分数据保证了模型对整个Ga-O二元相图的化学势有正确的理解，避免在模拟缺陷时出现能量严重失准的情况。

3. 无序体相结构：这包括通过熔融-淬火方法产生的非晶态Ga2O3，以及高温液态Ga2O3。我们生成了不同密度（从3.8到5.4 g/cm³）的非晶结构，以覆盖实验上可能的密度范围。液态构型则通过高温AIMD采样获得。这些低对称性、高能量的构型是训练模型“泛化”能力的关键，让模型学会描述远离晶体平衡态的原子环境。

4. 广泛随机结构搜索与分散构型：这部分可以理解为“撒大网捞鱼”。我们使用主动学习或随机结构搜索方法，在广阔的构型空间中采样，包括一些能量很高、原子排列非常随机的结构，甚至气态分散的原子团簇。这些数据虽然“质量”不一定高（能量高），但“多样性”极好，能极大地扩展模型适用的构型空间，特别是对于模拟高能粒子辐照等极端过程非常有帮助。

注意：数据库的“一致性”是生命线。所有这些构型的能量、原子力和维里应力，都必须由同一套经过严格收敛测试的DFT计算参数（我们用的是GGA-PBE泛函）产生。如果数据来源的精度不一致，比如有些用高精度的杂化泛函算，有些用普通的GGA算，那么训练出的模型就会精神分裂，预测结果不可靠。

2.2 两种GAP模型：精度与速度的权衡

有了高质量的数据，接下来就是选择机器学习模型。我们采用了高斯近似势（GAP）框架，它基于高斯过程回归，在材料科学领域已经非常成熟。但GAP模型本身也有不同的“配方”，核心差异在于如何描述一个原子的局部环境（即“描述符”）。

• 高精度之选：soapGAPsoaGAP的核心是使用了高维度的“平滑原子位置重叠”（SOAP）描述符。你可以把它想象成一个超级灵敏的“原子指纹识别器”。它不仅考虑一个原子周围邻居原子的种类和距离，还通过球谐函数等数学工具，精确刻画邻居原子的三维空间角度分布。这种高维描述符对原子环境的微小变化极其敏感，因此能实现非常高的拟合精度。我们的验证表明，soapGAP预测的能量和力与DFT参考值吻合得几乎完美，对于所有重要晶体相，能量误差的标准偏差仅在毫电子伏每原子量级。然而，高灵敏度带来的代价是计算复杂度。在模拟中，计算每个原子的SOAP描述符并进行高斯过程回归评估，是相当耗时的。它的计算速度大约在每秒每原子800个MD步长，这使得它适合用于几千个原子、时间尺度在纳秒量级的“中型”精细模拟。

• 速度之选：tabGAP当我们需要模拟百万原子、数十纳秒的大尺度过程时（比如研究薄膜生长初期的成核、或大尺寸缺陷簇的演化），soapGAP的速度就捉襟见肘了。为此，我们并行开发了tabGAP。它的思路是做“减法”和“查表”。tabGAP使用的描述符维度很低，主要是两体（2b）、三体（3b）和类似嵌入原子法（EAM）的密度描述符。这些低维描述符的区分能力不如SOAP，因此模型的绝对精度会有所下降（从我们的误差分布看，tabGAP的能量和力误差标准差大约是soapGAP的2-3倍）。但妙招在于，由于描述符维度低，我们可以预先将机器学习模型预测的能量和力“表格化”——也就是在低维网格上预先计算好所有可能的输入对应的输出值。在实际MD模拟中，只需要通过快速的三次样条插值来查表，而无需进行复杂的高斯过程计算。这一下子将计算速度提升了近5个数量级，达到每秒每原子约32万个MD步长，比soapGAP快了约400倍。

实操心得：如何选择？在实际研究中，我的策略通常是“混合使用”。先用tabGAP进行大范围的初步探索和筛选，比如寻找可能的相变路径、观察大致的动力学趋势。一旦锁定关键区域或时间节点，再换用soapGAP进行高精度的“特写镜头”分析，获取更可靠的定量数据（如精确的能量差、键长键角分布）。这种“先粗后精”的工作流，能极大提高研究效率。

2.3 训练技巧与短程作用处理

训练一个通用的ML势函数，绝不是把数据丢给算法就完事了，里面有很多“手艺活”。

• 正则化参数（σ）的智慧设置：机器学习模型容易“过拟合”，即对训练数据背得滚瓜烂熟，但对新数据预测很差。为了防止过拟合，我们引入了正则化参数σ，它本质上是告诉模型：“允许你对这类数据的预测有这么大的误差”。这里的关键在于不能一刀切。对于数据库中最核心、我们最关心的完美晶体相（部分(i)），我们设置了非常小的σ值（能量0.002 eV/atom，力0.01 eV/Å），强迫模型必须非常精确地拟合它们。而对于非化学计量比相和无序结构（部分(ii)(iii)(iv)），我们则放宽了要求，设置了更大的σ值（能量0.0035-0.01 eV/atom，力0.035-0.1 eV/Å）。这相当于告诉模型：“这些数据本身可能噪声大、多样性高，你不用追求完美复现每一个点，抓住主要趋势就行”。这种基于物理理解的、分门别类的正则化，是保证模型既在核心区域精准、又在广阔区域平滑的关键。

• 短程强排斥势的显式拟合：我们的势函数还有一个重要的应用场景：模拟高能粒子辐照引起的碰撞级联。在这个过程中，原子会被撞飞，两个原子可能靠得非常近，产生极强的排斥力。DFT计算这种超短程（<1 Å）相互作用非常昂贵且困难，而标准的GAP模型在训练数据稀少的超短程区域外推能力很差。我们的解决方法是“分而治之”。我们先用一个解析的、物理上合理的ZBL型排斥势对来描述超短程的强排斥作用。然后，在训练GAP模型时，我们不是让它直接学习DFT的总能量，而是学习DFT总能量与这个预设排斥势之差。这样一来，GAP只需要学习从短程到中长程的、相对平滑的相互作用部分，大大降低了学习难度，并保证了在整个距离范围内势能面的物理合理性。

3. 性能验证：当通用势函数遇上复杂多晶型

模型训练好了，是骡子是马得拉出来溜溜。验证一个势函数，不能只看它在训练集上的误差（那叫“考试作弊”），更要看它在各种真实的物理场景下，能否做出正确的预测。我们对soapGAP和tabGAP进行了一系列系统性的“期末考试”。

3.1 晶体结构性质的精确复现

首先是最基本的：五种多晶型的平衡晶格常数。下表展示了我们的GAP预测与DFT计算以及实验值的对比：

注：γ相是一种有本征无序的缺陷尖晶石结构，其晶格常数是多个随机胞的平均值。

可以看到，无论是soapGAP还是tabGAP，对所有五种晶体结构的晶格常数预测都与DFT结果高度一致，最大误差仅0.35%（出现在tabGAP预测的κ相a轴上）。这证明了我们的势函数具备了准确描述多种晶体结构的基础能力。

更严格的考验是相稳定性顺序。我们计算了五种相在均匀应变下的能量-体积曲线（即状态方程）。结果显示，两种GAP势都完美复现了DFT给出的相稳定性顺序：β相最稳定，其次是κ相，然后是α、δ，γ相能量最高。特别重要的是，它们准确地捕捉到了在高压下（对应体积收缩区域）β、κ、α、δ相能量非常接近（在±8 meV/atom以内）的“简并区”。这个区域对应着约24-44 GPa的压力，是实验上观察到α→β/κ→β固态相变的热力学平衡点。我们的势函数能复现这个微妙的能量平衡，意味着它有能力用来模拟高压下的相变过程。

3.2 动态与热学性质的预测能力

静态结构过关了，动态性质呢？我们通过分子动力学模拟来检验。

• 径向分布函数（RDF）：我们在900 K和0 bar条件下对每个相进行了NPT-MD模拟，并计算了RDF。无论是用soapGAP还是tabGAP，得到的RDF与基于DFT的AIMD结果几乎无法区分。RDF的第一峰（~2 Å）对应Ga-O平均键长，所有相都一致；而大于4 Å的长程峰形则反映了不同相中氧堆垛和Ga占位的对称性差异，我们的势函数也准确地再现了这些差异。

• 声子谱：声子谱是晶格动力学的全面反映，也是检验势函数“动态”精度的试金石。我们计算了β, κ, α, δ相的声子色散曲线。soapGAP的结果与DFT计算吻合得很好，仅在高频光学支有微小偏差。tabGAP的预测偏差稍大，特别是在低频声学支的某些点（如β相的Y点），但整体趋势和声子带范围仍然是正确的。对于研究热输运、热膨胀等与声子密切相关的性质，soapGAP是更可靠的选择；而对于关注结构演化、相变路径等过程，tabGAP提供的动力学图像在多数情况下也是可信的。

• 热膨胀系数：我们模拟了零压下不同温度时各相的晶胞体积变化。soapGAP预测的热膨胀率顺序为：β < α ≈ κ ≈ γ < δ，这与已有的第一性原理计算结果一致。tabGAP的预测在α和γ相上略有不同，它略微高估了α相的热膨胀，而低估了γ相的。这反映了低维描述符在区分某些复杂相时的局限性。需要特别注意的是γ相，由于其缺陷尖晶石结构的内在无序性，其在高温下的热膨胀行为偏离了简单的抛物线形式，我们的模拟也捕捉到了这一特点。

3.3 应对无序：液态与非晶态的挑战

一个“通用”势函数的终极试炼，是看它能否描述完全无序的体系。我们构建了密度为4.84 g/cm³（与2123 K实验液态密度一致）的液态Ga2O3模型，以及通过熔体淬火得到的非晶模型。

对于液态，soapGAP和tabGAP模拟得到的RDF、偏径向分布函数（PRDF）以及Ga-O键角分布，都与小体系AIMD的结果高度吻合。RDF中第一峰（Ga-O键）与第二峰（O-O和Ga-Ga最近邻合并）的强度比、峰位都一致，表明两种势函数都能准确描述液态中的短程有序（主要是离子性的Ga-O配位）。

对于非晶态，我们以不同的淬火速率（AIMD: 170 K/ps, soapGAP: 34 K/ps, tabGAP: 3.4 K/ps）获得了玻璃态结构。三种方法得到的非晶RDF也表现出极好的一致性：第一峰更强更窄，第一、二峰之间的谷接近零，并在~4.5 Å处出现很浅的第三峰（对应Ga-O次近邻壳层）。键角分布的主要特征，如60-70°区间的信号缺失、150°处的低肩峰，也被准确复现。尽管由于淬火速率差异，键角分布峰值有微小偏移，但关键趋势一致。这表明我们的势函数能够可靠地用于研究Ga2O3的非晶形成、晶化等过程。

踩过的坑：淬火速率的影响模拟非晶态时，淬火速率是一个不可忽视的参数。我们的模拟中，AIMD因为计算量限制，只能用极高的淬火速率（170 K/ps），这可能导致得到的非晶结构“冻结”了更多高温液态的特征，与实验上缓慢冷却得到的结构有差异。而用tabGAP则可以实现慢得多的淬火（3.4 K/ps），更接近真实条件。在比较不同方法得到的结构细节（如精确的键角、环状统计）时，必须考虑淬火速率的差异。我们的策略是，关注那些对淬火速率相对不敏感的普适性结构特征（如平均配位数、主要键长范围），这些才是势函数精度的真正体现。

4. 实战应用：揭示液固相变的三阶段微观机制

有了经过严格验证的通用势函数，我们就可以去探索那些用传统方法难以触及的物理过程了。我们选择了一个具有重要实际背景且机理复杂的问题：Ga2O3从液态到最稳定的β固相的相变过程。这对理解熔体法晶体生长（如导模法）的初始凝固机制至关重要。

我们设计了一个包含两个液固界面的模拟体系：中间是一块完美的β相（100）面 slabs，两端是经过预平衡的液态Ga2O3。整个体系包含11520个原子，采用周期性边界条件。我们使用计算高效的tabGAP，在1500 K和1 bar下进行了长达10纳秒的NPT-MD模拟。

模拟结果揭示了一个清晰的三阶段相变动力学过程，其能量演化与结构演变如下图所示（概念示意图）：

[液态区域] --- [界面区] --- [β相固体] --- [界面区] --- [液态区域] (阶段I: 慢速过渡) (阶段II: 快速转变) (阶段III: 仅Ga迁移)

阶段I：慢速过渡（0 - ~3000 ps）在相变初期，有序化过程从液固界面开始，缓慢地向液态内部推进，速度约为7 Å/ns。这个阶段最有趣的现象是氧（O）亚晶格和镓（Ga）亚晶格表现出了截然不同的行为。通过对原子位移分布的分析，我们发现一部分O原子在界面附近迅速排列成面心立方（fcc）堆垛序（这是β相中氧亚晶格的堆垛方式），并且一旦形成，这些O原子的运动就受到强烈约束，其位移分布集中在2 Å附近的一个尖峰内。相反，Ga原子的位移分布则很宽，表明它们在整个慢速过渡阶段都保持着较高的迁移率。这意味着，界面处的初始有序化是由氧亚晶格的快速“定位”主导的，而Ga原子则相对自由地在这些初步形成的氧骨架中随机占据四面体和八面体间隙位。

阶段II：快速转变（~3000 - ~3100 ps）当剩余的液态层厚度减小到约14 Å的临界尺寸时，相变过程突然加速。在大约100 ps的极短时间内，发生了一个明显的一级相变，剩余的液态区域迅速完成有序化。此时，整个体系的氧亚晶格已经完全转变为长程有序的fcc堆垛，但Ga原子的排列仍然存在大量缺陷，并未完全占据β相的正确位置。体系从一个“长程有序氧骨架+无序Ga”的状态，转变为“完全有序氧骨架+缺陷Ga”的状态。

阶段III：仅Ga迁移（3100 ps - 结束）在氧亚晶格完全锁定后，相变的后续过程就完全由Ga原子的迁移和重排来驱动了。从模拟快照可以看到，那些错位的Ga原子像缺陷簇一样，在已固化的氧骨架中集体缓慢扩散，逐渐移动到它们能量最低的平衡位置，最终形成完美的β相晶体。这个过程比前两个阶段要慢得多，因为Ga原子在固态中的扩散需要克服更高的能垒。

机制解读与意义这个三阶段机制的核心在于氧和镓亚晶格迁移率的巨大差异。在界面区域，氧离子由于带电性强（O²⁻），与周围离子的库仑相互作用强，一旦找到合适的fcc堆垛位置，就被牢牢“钉住”，迁移率很低。而Ga离子（Ga³⁺）虽然也带电，但在液态或界面这种松散环境中，其迁移受到的约束相对较小。因此，相变前沿的推进呈现出“氧先搭好骨架，Ga随后填空并调整”的异步协作模式。这个微观机制此前未被实验直接观测到，我们的模拟为其提供了原子尺度的理论解释。

实操心得：大尺度模拟的采样与分析

1. 体系尺寸：研究液固相变，界面面积和液态区厚度要足够大，以减小周期性边界条件的影响，并让成核、生长过程自然发生。我们的体系在界面法线方向长达100 Å，是合理的。
2. 温度与压力：1500 K高于Ga2O3的熔点（约2073 K），但我们在1 bar下模拟，过冷度足够驱动凝固。压力设置为0或1 bar，取决于你想模拟常压还是特定压力下的生长。
3. 分析方法：除了看总能量、温度、压力等宏观量随时间变化，原子位移分布函数是揭示不同物种动力学差异的利器。局域结构分析（如基于氧原子位置的常见邻位分析CNA）可以定量标定每个原子的“固态序”。可视化至关重要，定期输出原子构型并用VMD、OVITO等软件渲染，能直观地看到相变前沿的推进和缺陷的演化。

5. 常见问题与模型使用指南

在实际使用我们开发的soapGAP和tabGAP势函数进行研究时，可能会遇到一些典型问题。这里我结合自己的使用经验，做一个集中解答和提示。

5.1 模型选择与性能权衡

Q1：我的研究问题该用soapGAP还是tabGAP？

这是一个最常被问到的问题。选择取决于你的核心需求：

• 首选soapGAP的情况：

  • 需要最高精度：例如，计算相变能垒、缺陷形成能、表面能等对能量极其敏感的性质。
  • 研究声子相关性质：如热导率、热膨胀的精确计算。
  • 体系规模不大：原子数在几千以内，模拟时间在几纳秒以内。
  • 研究非化学计量比或极端缺陷：soapGAP对化学环境的敏感性更高，在成分偏离化学计量比时更可靠。

• 首选tabGAP的情况：

  • 大尺度模拟：原子数超过1万，甚至到百万原子级。
  • 长时间模拟：需要模拟数十纳秒以上的动力学过程，如薄膜生长、粗化、辐照损伤的长期演化。
  • 高通量筛选或初步探索：需要快速测试多种初始条件或参数。
  • 计算资源有限：tabGAP的速度优势可以让你在普通工作站或小型集群上完成可发表的研究。

Q2：tabGAP的精度损失到底有多大？会影响结论吗？

从验证数据看，tabGAP对于五种主要晶体相的静态结构（晶格常数、相对能量）预测非常准确，误差在可接受范围内。它的主要误差来源于对高能无序构型和非化学计量比区域的描述能力稍弱。因此：

• 如果你的研究聚焦于完美晶体或接近化学计量比的有序相变（如β/κ/α相之间的转变），tabGAP的精度足以支持定性甚至定量的结论。
• 如果你的研究涉及大量缺陷、非晶、或严重偏离化学计量比的情况，则需要用soapGAP的结果进行交叉验证，或至少对关键结论用soapGAP做检查。
• 对于液态、非晶等无序体系，我们的测试表明tabGAP在描述短程有序（RDF、配位数）方面表现依然良好，可用于研究结构演化趋势。

5.2 模拟设置与参数优化

Q3：使用GAP势进行MD模拟，有哪些关键的参数设置？

1. 截断半径（Cutoff）：这是最重要的参数之一。我们的GAP势函数内置了描述符的截断半径。绝对不要在MD模拟软件（如LAMMPS）中设置一个小于此值的截断半径，否则会丢失关键的相互作用信息。通常需要设置一个比GAP截断半径稍大（如0.5-1.0 Å）的邻居列表皮肤距离（skin），并确保及时更新邻居列表。
2. 时间步长（Timestep）：对于Ga2O3这样的氧化物，由于包含质量较轻的氧原子，建议使用较小的时间步长以确保能量守恒。通常0.5 fs到1.0 fs是安全的选择。开始可以先试用1.0 fs，监控系统总能量和温度的漂移情况，如果漂移过大，则减小到0.5 fs。
3. 热浴和压浴（Thermostat/Barostat）：对于NPT或NVT模拟，推荐使用较“温和”的控温控压方法，如Nosé-Hoover链（NVT）或Melchionna改进的Nosé-Hoover（NPT）。相关弛豫时间参数（如Tdamp,Pdamp）应设置为时间步长的数百到数千倍（例如，1 fs步长下，Tdamp设为100-500 fs）。
4. 长程静电相互作用：Ga2O3是离子性较强的材料。我们的GAP模型在训练时，其底层DFT数据已经包含了电子云分布带来的静电效应。因此，在MD模拟中不需要也不应该再额外添加Ewald或PPPM等长程静电求和。GAP势已经是一个有效的“短程”势，包含了所有必要的相互作用。

Q4：如何将GAP势函数集成到LAMMPS中运行？

我们的soapGAP和tabGAP都以quip或GAP库兼容的格式提供。以LAMMPS为例，大致步骤如下：

1. 编译LAMMPS时，需要激活ML-GAP或ML-QUIP包（具体名称取决于版本）。
2. 在LAMMPS输入文件中，通过pair_style命令调用GAP势。对于tabGAP，由于是表格形式，可能需要指定对应的参数文件。
3. 确保势函数文件（.xml或.gap文件）和参数文件在正确的路径下。 一个简化的LAMMPS输入脚本片段可能如下所示（具体命令需根据实际安装调整）：

# 定义单位、原子类型等 units metal atom_style atomic # 读取数据文件 read_data your_system.data # 设置配对势为GAP pair_style quip pair_coeff * * Ga2O3_soapGAP.xml "Potential xml_label=GAP_2023" 14 8 # 或对于tabGAP # pair_style table gap 10000 # pair_coeff * * Ga2O3_tabGAP.table GAP # 设置邻居列表 neighbor 2.0 bin neigh_modify every 1 delay 0 check yes # 定义计算 compute pe all pe/atom compute stress all stress/atom NULL virial # 运行能量最小化 minimize 1.0e-6 1.0e-8 1000 10000 # 运行NPT MD fix 1 all npt temp 300 300 0.1 iso 0 0 1.0 timestep 0.001 # 1 fs thermo 1000 run 1000000 # 1 ns

注意：这只是一个示例，实际使用时务必查阅相关软件包的具体文档，并正确设置原子类型映射、单位转换等。

5.3 结果分析与可靠性判断

Q5：如何判断我的MD模拟结果是否物理可靠？

1. 能量守恒检查（NVE系综）：对于微正则系综（NVE）模拟，系统总能量（动能+势能）应该是一个常数（除了极小的数值积分误差）。总能量的漂移应远小于系统动能或势能的涨落。如果发现明显漂移，首先检查时间步长是否过大。
2. 结构合理性：定期可视化原子构型。观察是否有不合理的键长（如两个原子异常接近）、原子飞离体系（“飞原子”现象）或晶体结构明显畸变。对于Ga2O3，可以快速计算一下平均Ga-O键长（应在1.8-2.0 Å左右）和配位数。
3. 与已知实验或DFT数据对比：对于你模拟的特定性质（如晶格常数、密度、RDF），与已有的可靠实验数据或高精度DFT计算结果进行对比。这是验证势函数在你所研究具体问题上是否可靠的最直接方法。
4. 系综平均值的收敛性：确保你的模拟时间足够长，使得所观测的物理量（如温度、压力、体积、RDF）的系综平均值已经达到稳定，并且统计误差足够小。
5. 敏感性测试：如果结论非常关键，可以尝试改变一些模拟参数，如初始构型、温度/压力控温器的参数、时间步长等，看主要结论是否依然成立。

Q6：在模拟相变时，如何确定相变是否发生以及相变路径？

1. 序参量（Order Parameter）：定义能够区分不同相的物理量。对于Ga2O3多晶型，常用的序参量包括：
   • 氧亚晶格的堆垛序：通过计算氧原子的位置并做傅里叶变换或直接分析其近邻层排列，可以识别fcc（β相）、hcp（α相）或4H堆垛（κ相）。
   • Ga原子的配位多面体：分析每个Ga原子周围氧原子的配位数（4或6）和具体排列方式。
   • 全局对称性指标：如计算体系的X射线衍射（XRD）图谱或电子衍射图谱，与已知相的图谱对比。
2. 局部结构分析：使用OVITO等软件中的Wigner-Seitz细胞分析、配位数分析、公共邻位分析（CNA）或键对分析等工具，对每个原子进行“染色”，从而在原子尺度上直观看到不同相的区域及其界面。
3. 热力学量的突变：在一级相变点，通常可以观察到体积、能量、焓等热力学量随温度或压力发生不连续的跳跃。
4. 可视化与动画：将模拟轨迹制成动画，是理解相变微观机制最直观的方式。你可以清晰地看到晶核如何形成、生长，界面如何移动，缺陷如何产生和湮灭。

开发并验证一套适用于复杂氧化物体系的通用机器学习势函数是一项艰巨但回报丰厚的工作。它就像为材料科学家打造了一台功能强大的“原子显微镜”，让我们能够以接近DFT的精度，去观察和操控百万原子、纳秒尺度的动力学过程。对于Ga2O3而言，soapGAP和tabGAP的发布，无疑将加速其在功率电子、深紫外探测、高温柔感等领域的材料设计与工艺优化进程。而我们所揭示的液固相变中氧、镓亚晶格异步演化的微观机制，不仅深化了对这一材料本征行为的理解，也为通过外场（如应力、电场）调控其相变和微观结构提供了新的理论线索。未来，这套方法可以进一步扩展到Ga2O3的掺杂体系、异质结界面、以及与其他宽禁带半导体（如Al2O3, In2O3）组成的多元复杂体系，在更广阔的材料设计空间中发挥威力。