手机版

钙钛矿分子动力学模拟 — 离子迁移与相变的微观机制

发布时间：2026-07-16 来源：科研学术网

小中大

字号：

钙钛矿分子动力学模拟在近五年从一个相对小众的方向变成了热门——太阳能电池效率一路飙升到26%以上，稳定性问题却一直卡脖子，而离子迁移被认为是导致器件退化的元凶之一。我从2019年开始做这个方向，从MAPbI₃做到FAPbI₃，再到CsPbI₃和二维Ruddlesden-Popper钙钛矿，最大的感受是：钙钛矿的分子动力学模拟比普通无机材料的MD难做太多，因为它的”软”——晶格软、有机阳离子会转、卤素离子会跑，这些动态行为对模拟方法的选择极其敏感。选错方法，你可能算出一个完全错误的离子迁移机制。

一、钙钛矿分子动力学模拟的特殊难点

钙钛矿分子动力学模拟有四个让做传统无机MD的人头疼的特点。

第一，有机阳离子的旋转动力学。MA（甲铵）和FA（甲脒）阳离子在钙钛矿笼里不是静止的，它们在皮秒到纳秒的时间尺度上快速旋转和重新取向。MA的C-N轴在室温下的旋转相关时间约1-5 ps，FA因为分子更大且有平面性约束，旋转慢一些但更复杂（面内旋转+面外翻转）。这个旋转行为直接耦合到PbI₆八面体的倾斜模式，进而影响能带结构和载流子动力学。你的模拟必须能正确描述这种有机-无机耦合，否则算出来的相变温度和离子迁移行为都会偏差。

第二，非谐效应极强。MAPbI₃的晶格热导率只有约0.3-0.5 W/mK，比大多数氧化物低一个数量级以上。这是因为PbI₆八面体的大振幅非谐振动和MA阳离子的旋转产生了极强的声子散射。经典MD力场（比如MYP力场）通常用简谐键伸缩项，对非谐效应的描述不够。AIMD能自然包含非谐性，但计算代价限制了可以跑的体系大小和时间尺度。

第三，离子迁移的能垒低。MAPbI₃中I⁻的迁移能垒只有0.3-0.6 eV，比氧化物中氧空位迁移的能垒（通常在1-3 eV）低得多。低能垒意味着在室温下就有可观的迁移速率——这是钙钛矿光伏器件不稳定的物理根源。但低能垒也给模拟带来挑战：你需要跑足够长的轨迹才能采样到足够多的迁移事件来统计扩散系数。如果轨迹太短，可能一个迁移事件都没观察到，得出”离子不迁移”的错误结论。

第四，相变序列复杂。MAPbI₃从低温到高温经历正交相→四方相→立方相的两次相变，相变温度分别在160 K和330 K左右。模拟相变需要捕捉到PbI₆八面体倾斜模式的序参量变化和MA阳离子的有序-无序转变，这对力场的精度和模拟时间都提出了很高要求。

二、力场参数化：有机阳离子（MA/FA）的电荷分配

钙钛矿分子动力学模拟中，力场参数化的第一个拦路虎就是有机阳离子的电荷分配。MA和FA是有机分子，带+1的净电荷，但这+1分布在C、N、H多个原子上。不同的电荷分配方案会导致完全不同的阳离子-无机框架相互作用强度，进而影响八面体倾斜、相变温度和离子迁移能垒。

经典的MYP力场（Mattoni等人提出的MAPbI₃力场）对MA阳离子使用了来自DFT计算的ESP（静电势拟合）电荷：N约-0.4e，C约+0.2e，H约+0.2e（分布在3个N-H和3个C-H氢原子上）。这个方案的问题是ESP电荷对构象敏感——MA旋转时原子所处的化学环境变了，电荷也应该变，但经典MD力场用的是固定点电荷。所以MYP力场在描述MA旋转和八面体倾斜耦合时有一定偏差。

后来的改进版本（比如Jiang等人的力场）引入了可极化力场的思路——至少对Pb和I使用Drude振子模型来处理电子极化。测试表明可极化力场对离子迁移能垒的描述比固定电荷力场更准——因为I⁻在迁移过程中经过的鞍点位置电子极化环境变了，固定电荷描述不了这种动态极化效应。

对于FA阳离子，情况更复杂。FA是一个共轭体系（N-CH-N骨架），π电子在整个分子上离域。用简单的原子点电荷分配无法准确描述它的静电势分布。我试过三种方案：ESP点电荷、CHELPG电荷、以及多极展开（包含原子偶极）。结果表明多极展开方案对FA旋转势垒的描述最好，但实现复杂度高。如果计算资源有限，至少应该用RESP（Restrained ESP）电荷——它在ESP电荷的基础上加了约束，减少了构象依赖性。

一个容易被忽视的细节：MA和FA的电荷分配必须在钙钛矿环境里做，不能在真空中做。真空中优化的MA阳离子构象和电荷分布，和在PbI₆笼里的完全不同。规范做法是把MA/FA放在一个足够大的钙钛矿团簇里做DFT计算，然后在嵌入环境下拟合电荷。

三、AIMD vs 经典MD：什么时候必须用第一性原理

钙钛矿分子动力学模拟中，AIMD（从头算分子动力学，通常是Born-Oppenheimer MD）和经典MD的选择是一个经典的精度-效率权衡。但钙钛矿体系有一些特殊考虑，让这个选择不那么简单。

先说AIMD的优势。AIMD能自然描述电子极化、电荷转移、化学键的断裂和形成——这些在经典MD里需要特殊的力场形式（比如ReaxFF）才能处理，而且参数化极其复杂。对于钙钛矿体系，AIMD在以下场景是必需的：

离子迁移的过渡态搜索和能垒计算。经典MD力场算出来的迁移能垒严重依赖力场参数，不同力场对同一体系能给出差一倍的能垒。只有AIMD（配合NEB或者constrained MD）能给出可靠的值。
相变过程的模拟。相变涉及PbI₆八面体倾斜模式的协同变化和MA/FA的有序-无序转变，这是一个高度非谐的过程。经典MD力场中的简谐键项和固定电荷模型很难准确描述。
有机阳离子的分解反应。MA和FA在高温或湿度条件下会分解（脱质子化、脱氨），这涉及化学键断裂，必须用AIMD或者反应力场。

但AIMD的局限也非常明显。用VASP跑AIMD，一个300-400原子的MAPbI₃超胞，在普通计算集群上跑100 ps需要几万核时。而研究离子迁移和相变往往需要纳秒级的轨迹——差了至少一个数量级。所以纯AIMD通常只能跑几十ps，只适合研究快过程（比如MA旋转动力学），不适合研究慢过程（比如I⁻的长程扩散）。

经典MD的优势在于可以跑到纳秒甚至微秒级别，能采样到更多的离子迁移事件。但代价是力场精度——特别是对迁移能垒的描述。我的实用策略是：用经典MD做大规模和长时间尺度的探索性模拟，然后用AIMD对关键结论做交叉验证。比如：经典MD发现某个温度下I⁻的扩散系数异常高，然后用AIMD在同样条件下跑一个短轨迹，看迁移能垒是否一致。

还有一个中间方案：机器学习力场（MLFF）。用AIMD数据训练的神经网络势能面可以同时保留AIMD的精度和经典MD的效率。现在VASP和CP2K都支持on-the-fly训练的MLFF，虽然训练过程需要额外的AIMD计算，但一旦训练完成，后续的MD可以跑到纳秒级。

四、离子迁移能垒计算中的采样陷阱

离子迁移能垒是钙钛矿分子动力学模拟中最关键也最容易算错的性质之一。能垒差0.1 eV，室温下的迁移速率就差将近50倍（根据Arrhenius关系exp(-ΔE/kT)）。所以能垒计算的精度要求非常高。

最常用的能垒计算方法是NEB（Nudged Elastic Band），在初始态和终态之间插入一系列image，优化找到最小能量路径。在钙钛矿体系里用NEB有三个容易踩的坑。

第一个坑：终态选择。I⁻从一个八面体顶点跳到相邻八面体的等价顶点，看起来只有一条路径，但实际上钙钛矿的非谐晶格会有多种可能的弛豫模式。如果你的初始态和终态结构没有充分弛豫（不仅仅是I⁻的位置，还有周围PbI₆骨架和MA阳离子的取向），NEB可能收敛到错误的鞍点。我遇到过一个案例，NEB给出的迁移能垒是0.35 eV，但检查中间image发现MA阳离子在鞍点处取向异常。把MA的初始取向调整后重算，能垒变成了0.52 eV。所以做NEB之前，先分别对初始态和终态做充分的离子弛豫。

第二个坑：image数量。NEB的image太少会漏掉能量路径上的细节（比如中间可能出现的小能垒），image太多则计算量爆炸。对于钙钛矿中的I⁻迁移，5-7个image通常够了。但如果迁移路径涉及多个原子的协同运动（这在钙钛矿里很常见），建议增加到9-11个image。

第三个坑：climbing image的使用。CI-NEB在找到鞍点后会让最高能量的image爬升到精确的鞍点位置，给出更准的能垒。但CI-NEB在某些情况下（尤其是势能面很平坦的时候）可能不收敛或者收敛到非物理的鞍点。建议先用普通NEB找到能量路径的大致形状，确认合理后再切换到CI-NEB。

除了NEB，还有一种基于MD轨迹的能垒估算方法：统计I⁻在不同位置的出现概率，构建自由能面（Potential of Mean Force, PMF）。PMF方法的好处是不需要预设迁移路径，能自动找到最低自由能路径。但PMF需要非常好的采样——如果某个高能垒区域采样不足，PMF会低估能垒。我一般用NEB确定迁移机制，用PMF验证热力学一致性。

五、从模拟轨迹提取扩散系数的正确姿势

钙钛矿分子动力学模拟做完后，从轨迹中提取I⁻的扩散系数是一个看似简单但细节很多的分析步骤。最常见的方法是算均方位移（MSD），然后通过Einstein关系D = MSD/(6t)（三维扩散）得到扩散系数。

第一个关键问题是MSD的时间区间选择。MSD-时间曲线通常有三个区域：弹道区（很短时间，MSD∝t²）、扩散区（中间时间，MSD∝t）、以及统计噪声区（长时间，数据点太少导致曲线剧烈波动）。只有扩散区的斜率能用来算扩散系数。判断扩散区的方法是画log(MSD) vs log(t)，扩散区斜率应该等于1。

对钙钛矿体系，扩散区通常在10-100 ps范围内（取决于温度）。如果轨迹只有200 ps，你可能只有不到100 ps的有效扩散区数据。这意味着每个温度点需要至少500 ps的轨迹才能得到统计可靠的扩散系数。如果你跑的是AIMD，这个时间尺度很难达到——这就是为什么很多文献里的AIMD扩散系数误差棒很大。

第二个关键问题是轨迹的分段平均。不要把整个轨迹算一个MSD然后线性拟合，这样会严重低估误差。正确的做法是把轨迹分成多个等长的时间窗口，每个窗口算一个MSD和对应的D，然后取平均和标准偏差。分段数越多统计越好，但每段要足够长以覆盖扩散区。一个经验法则是：每段至少是扩散区特征时间的3-5倍。

第三个问题：I⁻的扩散在钙钛矿里不是各向同性的——四方相和正交相中，面内和面外方向的扩散速率不同。所以算MSD时要分解到x、y、z三个方向分别算，得到各向异性的扩散张量。很多人只算总MSD，丢失了方向信息。

最后，算完扩散系数后要做Arrhenius拟合（ln D vs 1/T）得到迁移活化能Ea。Arrhenius拟合至少需要3-4个温度点的数据。如果某个温度点偏离线性趋势太远，不要硬拟合——可能那个温度下发生了相变或者扩散机制变化。钙钛矿的相变温度附近，扩散系数的Arrhenius行为通常会偏离线性。