1. 项目概述当机器学习“遇见”钙钛矿微观世界在太阳能电池材料的研究前沿甲脒碘化铅FAPbI3钙钛矿无疑是一颗耀眼的明星。它拥有理想的光学带隙、出色的载流子迁移率和较长的载流子扩散长度这些特性使其光电转换效率在实验室中已突破25%大关。然而与许多明星材料一样FAPbI3也面临着“稳定性”这一阿喀琉斯之踵。尤其是在低温环境下其晶体结构会经历复杂的相变从室温的立方相α相降至约285 K时转变为四方相β相最终在约150 K进入所谓的γ相。这个低温γ相长久以来就像一个蒙着面纱的谜团——实验上观测到其结构存在显著无序性但具体的原子排布、有机阳离子的运动状态以及无机骨架的扭曲模式一直缺乏清晰、统一的原子尺度图像。理解这个低温相至关重要。因为材料的宏观性能无论是光电转换效率还是长期稳定性都根植于其微观的原子结构和动力学行为之中。有机阳离子这里是甲脒离子FA的旋转是否被“冻结”铅碘八面体[PbI6]4-是如何倾斜的这些微观细节直接影响了材料的电子结构、声子谱乃至缺陷的形成与迁移。传统的研究手段如X射线衍射XRD在材料高度无序时可能“力不从心”而第一性原理分子动力学AIMD模拟虽然精确但其巨大的计算成本将模拟的尺度和时间限制在皮秒和数百个原子的量级难以捕捉到相变这类涉及长程有序和慢动力学的过程。这正是机器学习势函数Machine-Learned Potential, MLP大显身手的舞台。简单来说我们可以把它理解为一个极其聪明的“插值器”或“拟合器”。它通过“学习”大量由高精度量子力学计算如密度泛函理论DFT产生的数据构建一个能够以接近DFT精度、但计算速度提升数个数量级的原子间相互作用模型。有了它我们就能进行大规模数万原子、长时间纳秒级的分子动力学模拟以前所未有的清晰度“观看”材料在降温过程中每一个原子的实时运动轨迹从而揭示那些隐藏在实验数据背后的微观真相。本次分享我将带你深入一项前沿研究如何利用基于神经进化势NEP框架训练的机器学习势函数结合大规模分子动力学模拟彻底揭开FAPbI3低温γ相的神秘面纱。我们将一起拆解整个研究的技术路线从势函数的构建与验证到模拟策略的设计再到对八面体倾斜、阳离子取向等关键序参量的深度分析。更重要的是我会结合自己的模拟经验分享在运用这类先进工具时如何设置参数、解读结果以及避开那些容易让人栽跟头的“坑”。2. 技术核心机器学习势函数与模拟策略解析2.1 为什么是机器学习势函数在深入FAPbI3的低温世界之前我们必须先理解手中的“望远镜”和“显微镜”——机器学习势函数。传统上原子尺度模拟主要有两类方法基于量子力学的第一性原理方法如DFT和基于经验公式的经典力场Classical Force Field。第一性原理方法精度高无需预设参数是计算材料学的“金标准”。但其计算量随原子数呈O(N^3)甚至更快的增长进行一次包含百原子、几皮秒的AIMD模拟就可能消耗巨大的计算资源这使得研究相变、缺陷扩散等需要大体系、长时标的过程变得几乎不可能。经典力场通过预设的数学函数如Lennard-Jones势、库仑势描述原子间作用计算飞快可以轻松模拟数百万原子、纳秒甚至微秒的过程。但其精度严重依赖于力场参数的拟合质量。对于FAPbI3这种包含强耦合电子效应、氢键和动态无序的复杂混合体系开发一个普适且高精度的力场异常困难往往顾此失彼。机器学习势函数巧妙地走了“中间道路”。它的工作流程通常如下数据生成用DFT计算成千上万个不同原子构型包括各种键长、键角、晶格畸变、阳离子取向的能量、力和应力。这些构型需要尽可能覆盖材料可能出现的所有物理和化学环境。模型训练选择一个机器学习模型如本研究中使用的神经进化势NEP或流行的深度势能DeepMD将每个构型的原子坐标转化为一种数学描述符Descriptor用以表征每个原子的局部化学环境。然后训练模型使其能够根据这些描述符精准预测出DFT计算出的能量和力。验证与应用用模型预测一系列训练集之外的、已知性质的测试构型如不同相的晶体结构、弹性常数、声子谱等验证其迁移能力和精度。通过验证后即可将其植入分子动力学程序中进行大规模模拟。注意训练数据的质量和广度是MLP成功的生命线。如果训练数据没有充分涵盖你要研究的相变路径、缺陷构型或化学反应那么模型在这些区域的预测将是不可靠的甚至会产生完全错误的结果。这被称为“外推风险”。本研究采用的NEP势其参考数据来源于使用SCAN-rVV10泛函计算的DFT数据涵盖了纯FAPbI3、纯MAPbI3以及它们的混合体系。这种广泛的覆盖确保了该势函数能够准确描述从立方相到各种扭曲相的结构和能量 landscape。2.2 分子动力学模拟的“艺术”系综、参数与尺度有了可靠的势函数下一步就是设计模拟实验。分子动力学模拟并非简单地“按下运行键”其中每一步设置都蕴含着对物理过程的深刻理解。1. 系综选择NPT vs NVT模拟是在一个统计系综中进行的。为了研究温度-压力相图本研究采用了NPT系综恒粒子数、恒压、恒温。这个系综允许系统的体积即晶格常数随温度和压力变化而弛豫是研究结构相变最自然的选择。相比之下NVT恒体积系综会人为约束晶格可能无法正确反映相变时的体积变化。2. 控温与控压方法控温使用了Bussi-Donadio-Parrinello (BDP) 热浴。这是一种随机速度重标定的方法能产生正确的正则分布比简单的Berendsen热浴在理论基础上更严格。控压使用了随机晶胞重标定Stochastic Cell Rescaling, SCR方法。它通过对晶胞矢量引入随机扰动来实现恒压同样能产生正确的等温等压分布。3. 系统尺寸与边界条件研究使用了包含49,152个原子的超胞。这是一个非常关键的选择。相变特别是涉及长程有序的相变可能受到有限尺寸效应的显著影响。如果超胞太小模拟出的相变温度、有序参数可能会严重偏离宏观体系的真实值。通常需要通过测试不同尺寸超胞的结果是否收敛来确定一个“安全”的尺寸。本研究选择的超大超胞基本消除了有限尺寸效应带来的疑虑。4. 升降温速率动力学陷阱的“导演”这是本研究的精髓之一也是理解“亚稳态”为何形成的关键。在模拟中升/降温是通过逐步改变热浴的目标温度来实现的。文中提到的速率是6.34 K/ns。这意味着在1纳秒的模拟时间内系统温度变化6.34开尔文。为什么这很重要在真实实验中冷却速率可能是每分钟几开尔文甚至更慢这给了原子充足的时间通过热涨落翻越能量势垒找到能量最低的基态。但在MD模拟中受限于计算资源我们使用的冷却速率比实验快了数十亿倍。这种“快淬”过程使得系统没有足够的时间去探索整个能量面很容易被“困”在第一个遇到的局部能量极小值亚稳态中而无法到达全局最小值基态。模拟中的策略研究者通常会尝试不同的冷却速率。如果发现一个结构在极慢的冷却速率下计算上极其昂贵才能得到而在较快的速率下得到另一个结构那么前者更可能是热力学稳定的基态后者则是动力学冻结的亚稳态。本研究通过对比不同速率下的结果确认了所观察到的现象不是模拟伪影而是反映了真实的物理动力学过程。3. 低温相的结构揭秘从全局搜索到动态演化3.1 基态结构搜索穷举法与Glazer记号在开始动态模拟之前首先要回答一个基本问题在绝对零度0 K下FAPbI3最稳定的晶体结构是什么这就是基态结构搜索。研究采用了一种“暴力但有效”的穷举采样法以一个立方相的原胞包含12个原子为基础构建2x2x2的超胞。在这个超胞中随机化两种自由度a) 所有FA阳离子的空间取向b) 所有[PbI6]八面体的倾斜模式和幅度。这产生了约100万个不同的初始原子构型。使用机器学习势函数对每一个构型进行能量最小化结构弛豫直到原子受力小于0.1 meV/Å。分析所有弛豫后结构的能量找到能量最低的那个。结果如图1所示能量最低的结构被确定为a⁻b⁻b⁻相采用Glazer记号描述。这里简单解释一下Glazer记号它用上标“”、“-”、“0”来描述相邻八面体沿某个晶轴倾斜的相互关系。“”表示同相倾斜所有八面体朝同一方向倾“-”表示反相倾斜相邻八面体倾斜方向相反“0”表示无倾斜。因此a⁻b⁻b⁻表示在a轴方向是反相倾斜在b和c轴方向也是反相倾斜注意在正交晶系中a, b, c轴不等价。这个a⁻b⁻b⁻结构图1B中所有的FA阳离子都指向同一个方向无机骨架呈现高度有序的倾斜模式。它的能量比紧随其后的竞争者如a⁰b⁻b⁻还要低稳坐基态宝座。3.2 升降温模拟捕捉相变与滞后现象接下来让这个系统“活”起来。从a⁻b⁻b⁻基态开始加热同时从高温立方相开始冷却观察其相变路径。模拟结果图2清晰地展示了一个非对称的、带有滞后现象的相变过程加热路径从基态开始a⁻b⁻b⁻ (GS) → ~190 K → β相 (a⁰a⁰c⁺) → ~315 K → α相 (a⁰a⁰a⁰)。这与实验上观测到的加热过程相符。冷却路径从高温相开始a⁰a⁰a⁰ (α相) → ~285 K → a⁰a⁰c⁺ (β相) → ~120 K →a⁻a⁻c⁺ (γ相)。关键发现来了冷却最终得到的不是能量最低的a⁻b⁻b⁻基态而是一个能量高出约2 meV/atom的a⁻a⁻c⁺结构这就是动力学陷阱的典型表现。系统在快速冷却过程中没有足够的时间完成从β相a⁰a⁰c⁺到基态a⁻b⁻b⁻所需的复杂重构特别是沿c轴从同相倾斜转变为反相倾斜而是选择了一条“捷径”在a和b轴引入了反相倾斜但保留了c轴的同相倾斜从而冻结在了a⁻a⁻c⁺这个亚稳态中。从图2的晶格参数、热容和能量曲线可以明显看到β→γ的相变在冷却过程中是一个连续变化二级相变特征而在加热过程中从亚稳态γ加热则表现出更尖锐的变化一级相变特征进一步印证了滞后现象。3.3 八面体倾斜模式分析可视化结构演化如何定量描述这种倾斜研究计算了每个[PbI6]八面体相对于理想立方位置的欧拉角ψ, θ, φ从而得到其在三个方向上的倾斜角度。图3的温度-倾斜角分布图如同一幅“相变地图”高温区285 K所有倾斜角集中在0度附近对应无倾斜的立方a⁰a⁰a⁰相。中温区285 K ~ 120 Kψ角出现一个约10度的分布而θ和φ角仍为0。这对应a⁰a⁰c⁺相即仅在c轴方向发生同相倾斜。低温区120 Kψ角分布增强至约15度同时θ和φ角也出现了约5度的分布。这表明在a和b轴方向也出现了反相倾斜结构最终演变为a⁻a⁻c⁺。这个分析直观地展示了冷却过程中无机骨架的连续演化先是在c轴“集体歪头”同相倾斜然后在更低的温度下在a和b轴方向开始“错落有致地歪头”反相倾斜但c轴的倾斜模式被“锁定”了未能翻越势垒转变为能量更低的反相模式。4. 有机阳离子的“舞蹈”与“冻结”无机骨架的倾斜只是故事的一半。作为混合钙钛矿的灵魂有机阳离子FA的取向和动力学行为对材料性能有着至关重要的影响。4.1 取向分布从自由旋转到部分冻结为了表征FA的取向研究者选取了两个特征向量连接两个氮原子的矢量r_NN和连接碳与氢的矢量r_CH。通过统计在不同温度下这两个矢量在空间中的角度分布概率P(θ, φ)可以清晰地看到阳离子有序度的变化。图4的结果非常精彩330 K (a⁰a⁰a⁰相)r_NN和r_CH的取向分布几乎是均匀的球面分布图4A, D。这说明在高温立方相中FA阳离子在笼子里近乎自由旋转像一个高速旋转的陀螺。200 K (a⁰a⁰c⁺相)分布开始出现图案图4B, E。r_NN矢量明显倾向于指向晶格的[100]和[010]方向即ab平面内的两个轴。这是因为无机骨架在c轴倾斜后笼子的形状从立方体变成了四方柱限制了FA的某些取向。10 K (a⁻a⁻c⁺相冷却得到)分布图案变得更加尖锐和复杂图4C, F。大部分FA的取向仍“继承”自上一相a⁰a⁰c⁺的偏好方向但同时出现了四个对称的“翅膀”这对应着理想a⁻a⁻c⁺结构中FA的可能取向图1C中的蓝点。最重要的对比图中用红点标出了基态a⁻b⁻b⁻相中FA的取向。在冷却得到的亚稳态结构中这个取向出现的概率极低。通过计算自由能势垒F -k_B T ln(P)发现FA要旋转到基态的取向需要克服超过100 meV/FA的高能垒。这就像FA分子被“冻”在了a⁰a⁰c⁺相遗留下来的取向附近无法跨越鸿沟到达能量更低的基态取向。4.2 旋转动力学时间尺度的测量取向分布是静态快照而动力学则告诉我们这些分子“动”得有多快。通过计算取向自相关函数ACFC(τ)可以量化FA旋转的快慢。C(τ)衰减得越快说明分子重新取向的速度越快。从图5可以看到在高温相C(τ)在几个皮秒内就衰减到零表明快速旋转。随着温度降低衰减速度变慢。在低温a⁻a⁻c⁺相C(τ)衰减得非常缓慢意味着FA的旋转运动几乎被“冻结”。通过将ACF的衰减拟合为指数形式可以提取出旋转特征时间τ_rot。如图6所示τ_rot随温度降低呈指数增长符合阿伦尼乌斯公式。一个关键验证将模拟得到的C-H矢量旋转活化能与实验值来自文献对比发现两者吻合得很好见表I。例如在a⁰a⁰c⁺相模拟得到的C-H旋转活化能为48.5 meV与实验值45 meV非常接近。这不仅验证了机器学习势函数的可靠性也确认了模拟中捕捉到的动力学行为是物理真实的。对基态的洞察研究还特别模拟了真正基态a⁻b⁻b⁻相中的FA动力学。发现即使在120 K下其ACF在长达10纳秒的模拟时间内也几乎不衰减C(τ) ~ 1估算其旋转时间至少超过20微秒。这与实验在低温γ相中观测到的纳秒级旋转动力学截然不同。这提供了最强有力的证据实验上观测到的低温相不是高度有序、阳离子完全冻结的a⁻b⁻b⁻基态而正是我们模拟中得到的、部分无序且阳离子尚存缓慢旋转的a⁻a⁻c⁺亚稳态。5. 与实验的握手NMR与INS光谱验证计算模拟的最终价值在于能否解释和预测实验现象。本研究通过核磁共振NMR和非弹性中子散射INS两种对局部结构极其敏感的谱学技术进行了出色的交叉验证。5.1 核磁共振NMR捕捉局部化学环境的“指纹”研究者对FAPbI3单晶进行了低温魔角旋转固态¹³C和¹⁵N NMR实验95 K。¹³C谱在三次连续的“冷冻-解冻”循环中谱图几乎不变图7A表明碳原子的局部环境相对稳定。¹⁵N谱情况大不相同图7B。谱线由多个重叠的信号峰组成并且这些峰的相对强度在三次循环中发生了微妙变化。这说明氮原子所处的局部化学环境存在分布且这种分布在每次快速冷却淬火时都可能略有不同——这正是结构无序和动力学冻结的特征。为了解读这些谱线研究者用DFT计算了不同结构中的¹⁵N化学位移对于有序的a⁻b⁻b⁻基态所有氮原子等价理论上只应产生一条尖锐的谱线。对于冷却得到的无序a⁻a⁻c⁺结构由于FA取向多样氮原子处于多种不同的局部环境中计算出的化学位移呈现一个宽分布。尽管绝对位移值因未考虑自旋轨道耦合等因素存在系统偏差但计算明确显示只有无序结构才能产生与实验观测相符的宽谱。对实验谱进行分峰拟合发现它可以用8个不同位置的氮信号来很好地描述图7D-F这与模拟中观察到的FA取向集中在有限几个方向的结果定性一致。5.2 非弹性中子散射INS探测氢原子的振动谱INS对氢原子运动极其敏感能提供材料的振动谱声子态密度。研究者计算了三种结构的INS谱图8理想a⁻b⁻b⁻和理想a⁻a⁻c⁺结构由于FA高度有序局部环境单一计算出的谱图由一系列尖锐的峰组成。冷却得到的a⁻a⁻c⁺结构由于FA取向无序氢原子处于多种不同的环境中导致谱峰显著展宽。与实验谱10 K对比实验谱是宽泛的“馒头峰”与冷却得到的无序a⁻a⁻c⁺结构的计算谱高度吻合而与两个有序结构的尖锐谱峰截然不同。这构成了一个完整的证据链模拟预测的无序a⁻a⁻c⁺亚稳态其静态结构特征NMR化学位移分布和动态/振动特征INS谱峰展宽都与实验观测完美匹配。而理论上能量最低的有序a⁻b⁻b⁻基态其预测谱图与实验严重不符。这强有力地证明实验上难以测得的FAPbI3低温γ相其真实面目正是这个被动力学冻结的、部分无序的a⁻a⁻c⁺亚稳态。6. 实操启示与避坑指南基于这项研究以及我个人在计算材料学领域的经验我想分享几个在运用机器学习势函数研究复杂材料相变时的关键点和常见陷阱。1. 势函数的选择与验证是重中之重不要做“黑箱”用户在使用一个现成的MLP之前务必了解其训练数据的来源、范围和局限性。检查它是否在你要研究的相变区间、压力范围、缺陷类型上有充分的训练数据。本研究使用的NEP势其训练数据明确包含了FA/MA混合体系的各种构型这为其在FAPbI3相变研究上的可靠性奠定了基础。必须进行基准测试即使论文声称势函数很好也要用自己的体系做简单的基准测试。例如计算一下常温常压下晶格常数、弹性常数、几个关键声子频率与已知的实验值或高精度DFT结果对比。如果这些基本性质都对不上后续复杂的动力学模拟结果将不可信。2. 分子动力学模拟的参数敏感性时间步长本研究用了0.5 fs对于含氢体系是典型的安全值。步长太大会导致能量不守恒甚至模拟崩溃太小则浪费计算资源。一般原则是小于最轻原子振动周期的十分之一。升降温速率这是影响相变结果最关键的参数之一。永远要意识到你模拟的速率远快于实验。如果条件允许尝试不同的冷却速率。如果某个相变现象在某个临界速率以下出现以上消失那么这个现象很可能与动力学过程相关。本研究通过对比确认了a⁻a⁻c⁺亚稳态的形成是动力学冻结而非模拟伪影。系统尺寸对于涉及长程有序如反相倾斜的相变必须做尺寸收敛性测试。模拟一个2x2x2的超胞和一个4x4x4的超胞看相变温度、有序参数是否变化很大。本研究使用超大超胞等效于16x16x16原胞虽然计算代价大但结果令人信服。3. 数据分析从海量轨迹中提取物理定义清晰的结构序参量对于钙钛矿八面体倾斜角如本研究用的欧拉角和有机阳离子取向如N-N矢量是核心序参量。需要编写可靠的脚本如利用OVITO、MDAnalysis等工具包来自动化分析每一帧轨迹。警惕“平均”的误导在无序相中平均结构可能掩盖重要信息。例如低温γ相的平均倾斜角可能很小但实际每个八面体都有显著的瞬时倾斜只是方向杂乱。因此像本研究一样分析分布函数图3, 图4比只看平均值更有价值。关联函数与时间尺度计算自相关函数时确保模拟时间远长于你要测量的特征时间τ_rot。如果ACF在模拟结束时还未衰减到零你只能给出τ_rot的下限如“ 10 ns”就像本研究中对基态相的分析那样。4. 与实验对比的艺术对比“苹果和苹果”将模拟的“快淬”结构与实验的“快淬”样品对比是合理的。但如果你想对比平衡态性质就需要极慢的模拟速率或增强采样技术。利用多种实验手段单一实验技术可能有局限。XRD对长程有序敏感但对高度无序的相解析困难正如本文提及的γ相XRD难题。NMR和INS对局部环境和轻元素敏感是验证模拟中原子尺度细节的利器。成功的计算研究往往能统合多种实验表征形成一个自洽的图像。承认并量化不确定性MLP有预测误差DFT计算化学位移存在系统误差如忽略自旋轨道耦合模拟时间尺度有限。在讨论中明确指出这些局限性并说明它们如何影响结论的强弱。例如本研究明确指出计算的¹⁵N化学位移绝对值有偏差但分布趋势是可靠的。这项研究为我们提供了一个范本如何将前沿的机器学习势函数、精心设计的大规模分子动力学模拟、深入细致的序参量分析以及多种实验谱学技术紧密结合最终破解了一个困扰领域多年的微观结构谜题。它告诉我们FAPbI3的低温世界并非静止不变的热力学基态而是一个被动力学“冻结”在途中的、充满复杂无序和缓慢弛豫的丰富景观。理解这一点对于设计更稳定的钙钛矿材料例如通过添加剂或应变工程来调控相变路径、避免有害亚稳态的形成具有根本性的指导意义。计算模拟不再是实验的附庸它已经成为在原子尺度上“设计”实验、解读现象、甚至预测新物理的强大引擎。