《储能科学与技术》专刊推荐|施思齐等：基于蒙特卡罗模拟的离子导体热力学与动力学特性

作者：刘金平 ¹ 蒲博伟 ¹邹喆乂 ²李铭清 ¹丁昱清 ¹任元 ¹罗亚桥 ¹李杰 ³李亚捷 ¹王达 ¹何冰 ⁴施思齐 ^1,5,6 

第一作者：刘金平（1994—），男，硕士研究生，从事超离子导体离子输运特性的研究，E-mail：[e-mail]；

通讯作者：施思齐，教授，研究方向为电化学储能材料基础科学问题解析、计算方法发展和新材料设计研究，E-mail：[e-mail]。

单位：1. 上海大学材料科学与工程学院，上海 200444；2. 湘潭大学材料科学与工程学院，湖南 湘潭 411105；3. 加州大学欧文分校物理与天文系，美国 加州 92697；4. 上海大学计算机工程与科学 学院；5. 上海大学材料基因组工程研究院，上海 200444；6. 之江实验室，浙江 杭州 311100 

引用本文： 刘金平,蒲博伟,邹喆乂等.基于蒙特卡罗模拟的离子导体热力学与动力学特性[J].储能科学与技术,2022,11(03):878-896. 

LIU Jinping,PU Bowei,ZOU Zheyi,etal.Investigating thermodynamic and kinetic properties of ionic conductors via MonteCarlo simulation[J].Energy Storage Science and Technology,2022,11(03):

878-896. 

DOI：10.19799/j.cnki.2095-4239.2022.0050

摘 要 基于概率统计理论的蒙特卡罗模拟(MC)在20世纪40年代由冯·诺伊曼等提出，其作为一种重要的数值计算方法，已被广泛应用于离子导体的热力学、动力学等性质的研究。然而，MC模拟在相关计算精度、模拟速度以及模拟流程自动化等方面，仍具有较大的提升空间。本文通过分析其在离子导体计算领域中哈密顿量的构建模式(如基于键价和计算或利用团簇展开将通过拟合代表性小尺寸晶胞第一性原理计算总能得到的近邻作用参数给出构型能量)以及结构演变方式(如基于单离子迁移模式假设的构型演变)，提炼出一套用于分析离子导体离子输运特性及相变特性的MC模拟范式，并给出与之相关的半自动化MC模拟程序，预测了石榴石结构离子导体的电导率与迁移离子占据率分别随锂离子浓度的变化趋势。为了更好地拓展MC模拟在离子导体研究中的应用，本文就其在包括正负极材料、电解质及其相关界面的电化学储能材料中的典型热力学与动力学计算案例进行了剖析，包括求解离子导体中的离子扩散问题、模拟迁移离子的分布特征与相关界面层的演变过程等。最后，展望了MC方法目前面临的挑战并给出可能的对策，包括：①精确地获取所有可能发生的事件(如单离子迁移与双离子协同迁移)及其描述(如哈密顿的计算)；②探寻高效算法以精确找到系统的演化轨迹；③精确获得MC模拟中对应的时间步长。关键词 蒙特卡罗模拟；离子导体；离子输运特性；相变特性；晶体生长

明确离子导体的离子输运与相变特性是制备新型全固态电池、电化学转换器件等离子型器件的关键。离子导体的典型结构包括在输运网络中嵌入或迁移实现电荷储存与传输的迁移离子以及构建离子输运网络的仅原点附近振动的骨架离子。高离子电导率的超离子导体的筛选过程通常采用高通量计算方法[如几何分析方法(crystal structure analysis by voronoi decomposition，CAVD)、配位成键理论(bond valence site energy，BVSE)等]分析骨架离子中的离子输运网络。同时，借助基于牛顿力学的分子动力学(molecular dynamics，MD)模拟或基于离子迁移实验的动力学蒙特卡罗(kinetic Monte Carlo，KMC)模拟分析迁移离子的输运特性(例如离子电导率、扩散系数等)。借助群论、巨正则蒙特卡罗(grand canonical Monte Carlo，GCMC)模拟等热力学模拟手段研究离子导体的结构演化。分子动力学(MD)模拟作为典型的动力学模拟手段因其计算精度高、模拟分子运动过程清晰等优势，在离子扩散模式分析及晶体生长模拟等方面具有广泛应用。然而其内禀缺陷限制了MD适用于更广泛的应用场景。例如，第一性原理分子动力学模拟计算速度较慢、难以应用到较大体系的快速模拟；而平均场分子动力学模拟不易获取其所需的高精度势函数。因此，具有类似计算功能的KMC模拟因计算速度快、扩展性强而被广泛应用于评估离子导体的离子输运特性。表1对几种典型的MD模拟与MC模拟进行了对比，尽管两者的应用范围存在众多交叠，但基于计算原理的差异性，MC模拟具有更广泛的时空尺度。此外，在热力学模拟领域(如相图计算、电压平台计算等)，相较于第一性原理相图计算以及空间群分解等模拟手段，MC模拟同样在计算速度与适用性等方面存在巨大优势。

表1   MC模拟与MD模拟的对比

目前，MC模拟还存在计算精度不高、计算过程依赖于较多的输入参数、蒙特卡罗时间与真实时间之间存在差异、无法捕捉细小的离子振动等不足之处。为促进MC模拟的发展与应用，本文介绍了其发展历史，提炼出MC模拟的基本模拟范式与具体计算步骤，以期为科研工作者开展离子导体的相关物化特性研究铺平道路。本研究给出了一套具体的KMC计算程序，并对石榴石型固态电解质的离子输运特性进行了分析，供读者参考。为清晰地了解MC模拟的适用情形，本文就几种典型的离子导体材料(锂离子电池正极材料、负极材料以及固态电解质材料等)介绍了MC模拟的计算结果。最后，本文对当前MC模拟所面临的挑战以及未来可能突破的方向进行了总结，期待该方法的进一步发展。

1 蒙特卡罗模拟方法基础

1.1　蒙特卡罗模拟方法

MC模拟作为一种基于随机抽样统计的模拟手段，由美国数学家冯·诺伊曼等于20世纪40年代在美国的“曼哈顿计划”中提出。到20世纪50年代，进一步提出了MC模拟中可适用于大多数物理平衡系统的著名算法：Metropolis算法。20世纪60年代，KMC模拟开始出现。与此同时，混合蒙特卡罗(hybrid Monte Carlo，HMC)采样算法等也得到了广泛地研究。进入21世纪以来，随着人工智能的快速崛起，蒙特卡罗树搜索算法等更为先进实用的计算方案也得到了快速发展(图1)。根据研究对象的不同，MC模拟可分为两类：一类是针对热力学平衡过程中非时间关联性问题的模拟，如常见的GCMC模拟；另一类是对时间关联的动力学问题的研究，如KMC模拟。这两种MC模拟的最终目标都是尽可能地寻找徘徊在某一低势能面附近的结构组态。从原理上讲，GCMC与KMC之间存在应用场景的重合，即巨正则体系依旧可以进行KMC模拟，反之亦然。多数KMC模拟与热力学MC模拟均可以利用格子气模型建模，故而上述多数模拟均属于格子气蒙特卡罗(lattice-gas Monte Carlo，LGMC)模拟。而格子气模型的MC模拟多属于典型的马尔科夫链(Markov chain)过程，所以这种MC模拟也可以叫作马尔科夫链蒙特卡罗(Markov chain Monte Carlo，MCMC)。这几种模拟手段之间的关系可以通过图3(e)所示的维恩图表示。基于上述几种MC模拟的计算过程大多相似，本文将在格子气MC模拟的介绍中对体系描述、体系演变及新体系接受概率等细节问题进行重点介绍。图1   蒙特卡罗模拟的发展
图2   MC模拟中常见的几种二维格子气模型：(a) 二维正方格子气模型；(b) 二维六边格子气模型；(c) 二维随机格子气模型(紫色的球代表格点上的离子)
图3   MC模拟中的几种能量近似计算：(a) 石墨结构中的锂离子与其他离子之间的相互作用能的计算；(b) 迁移离子的迁移过程中所需要克服的迁移能垒示意图；(c) 迁移离子的不同迁移模式；(d) 不同迁移模式对应的迁移能垒示意图；(e) 几种蒙特卡罗模拟之间的关系

1.2　蒙特卡罗模拟过程

GCMC模拟作为一种适用于巨正则系综限制下经典的热力学MC模拟，被广泛应用于电化学储能材料的合成、多孔吸附等过程的研究。在电化学储能材料计算领域，GCMC常用于模拟石墨等碳基嵌锂材料的结构特征。在模拟过程中，随机选择几种基本的离子运动模式(比如置换、插入、脱出)，并对相应离子运动模式发生的概率进行计算，一般可以通过式(2)所示的计算模式求得。值得指出的是，式中的能量变化不仅可以使用体系的始末态能量表示，也可以使用过渡态能量表示。对于忽略动力学特性的热力学模拟，仅考虑简单的始末态能量是常见的。作为GCMC模拟适用情形的补充，KMC模拟则可以很好地描述离子导体结构演变的动力学过程。KMC模拟的研究最早可以追溯到20世纪60年代对辐射损伤退火问题的模拟。目前，KMC在表面吸附、离子扩散、晶体生长等领域得到了迅速发展和应用。由于模拟机理不同，KMC相对于分子动力学(MD)模拟可以对更大系统的动态演化过程给出更长时间的模拟。离子导体中迁移离子的成功迁移并不频繁，其状态转变会经历长时间的“尝试”，只有极少的“尝试”会成功。KMC模拟不考虑尝试过程只关注可以演变成功的状态，这种状态到状态的转变不受原子振动的限制，不局限于对每个原子(离子)的计算，因此可以实现更长时间尺度的模拟。同时，KMC的离散化思想更适合对微观离散体系进行模拟分析，这也使得离散化模型的搭建在KMC模拟中占据着十分重要的地位。LGMC模拟即为一种典型的离散化处理方式。GCMC模拟与KMC模拟均可以借鉴格子气模型而发展成为一种特殊的LGMC模拟来描述离子在材料体系中的扩散与占位情况。在格子气模型的搭建过程中，时间与空间被分割，按照离子之间的相互作用以及迁移、翻转规则，将模拟体系简化成一系列固定的空间格点与可迁移或翻转的离子。图2展示了几种典型格子气模型的示意图。该模型中，离子位于格点上，沿着边运动，离子与离子之间存在相互作用(比如吸引或者排斥作用)。Ising模型、渗流模型等均是对格子气模型的某种处理与变形。在研究离子导体材料时，MC模拟通常使用体系的能量来确定体系的演变过程。体系能量则通过描述近邻迁移离子间以及迁移离子与骨架离子之间的相互作用得到，这样既可以简化能量的计算过程也可兼顾描述的完整性。常用的能量计算公式如式(1)所示

（1）

式中，m为总的位点数目；为第j个位点被占据所引起的能量变化；为最近邻离子间的相互作用能；和为前置系数(其值取决于结构中离子的排布状态)。如果位点“j”被占据，则，反之。如果位点“j”和位点“i”同时被占据，并且它们互为最近邻位点，则，反之。另外，由Metropolis等提出的系统演变概率计算公式具有精度较高、可有效避免结构在演变过程中进入局域能量极小值等优点，所以在MC模拟中得到广泛应用(利用此算法的MC模拟也可以叫做Metropolis Monte Carlo，MTMC)。依照此算法，组态a到组态b的演变概率(为

（2）

式中，k为玻尔兹曼常数；T为温度，为系统从组态a演变到组态b的能量变化。计算公式如下

（3）

式中，)为组态a(b)的能量，值得指出的是，在动力学模拟中也会考虑到过渡态，如迁移能垒的引入等。仅依托于上述内容来计算离子的扩散系数、离子的迁移率等时间相关量是困难的。为了进一步描述材料结构随时间的一系列演变过程进而计算相邻组态之间的时间关联性，常常通过数学近似的方式获得真实的时间步长

（4）

式中，表示“总迁移概率”；是随机数[∈(0,1)]。由于总概率的计算十分耗时，而在MC模拟过程中一个蒙特卡罗步(MCS)内的离子迁移可认为是同步发生的，所以采用MCS取代真实时间是一种有效的简化计算方案。这种简化已在传统的热力学计算中被广泛采用。对于MC理论在材料学中的应用，其基本思想是将时-空尺度进行分离，因而其主要用于模拟有序系统，很难应用于无定形系统。对于MC模拟中相关名词的概念与定义在表2进行了总结。

表2   MC模拟中的常用概念

2 蒙特卡罗模拟的计算范式归纳与示例

MC模拟已被应用到众多领域，但是如何提高其模拟精度与效率一直是众多科研工作者的重要研究方向。为了清晰了解蒙特卡罗模拟的计算过程，此章节就其计算流程和步骤进行了总结与提炼。一般来说，模型的搭建应尽可能地复现材料的结构特征，比如通过几何分析方式或BVSE来分析材料框架以及离子传输路径，通过格子气模型分析常见的晶态离子导体，通过Ising模型来分析磁性材料等。此外，合适的模型描述也是决定相关模拟精度的关键。通常研究者会通过对模型能量(如内能)的计算，进而以能量的变化来确定演变方向以及对应构型出现的概率。然而，目前还没有一个完美的计算方案可以准确地计算材料构型转变过程中的能量变化。为此，众多科研工作者作了许多努力。比如Arriazu等定义了如图3(a)所示的晶格气体模型，考虑了层内Li离子、层间Li离子以及Li离子与骨架结构的相互作用能，利用GCMC模拟了Li离子在石墨中的插层过程。在模拟中，作者采用一种纯静态的能量计算模式，即只根据系统的总能来确定此构型出现的概率。因其不考虑过渡态能量，所以该模拟方法在动力学模拟中是不够精确的。对于固定的迁移或反转，激活能垒更能准确反映其能量的变化，因此Ven等提出了如图3(b)所示的离子迁移能垒。为避免因离子往返迁移能垒不同而造成的计算困难，其引入了一个动力学解析的激活势垒(∆E_KRA)作为能垒的近似，同时还考虑了目标离子周围离子的排布，采用团簇展开的策略，近似计算了离子的迁移概率。然而该近似不仅忽略了离子迁移的方向性，也忽略了多离子迁移模式。迁移模式的选择(材料结构演变途径)是MC模拟需要解决的第三个关键问题。He等证明了超离子导体中的离子快速扩散并不像典型固体中那样通过孤立的离子迁移发生，而是通过具有低能量势垒的多离子协同迁移进行[图3(c)～(d)]。实际上，迁移模式的选择是激活能计算的前提，而目前大多数报道都是基于单离子迁移模式，这会造成模拟结果的偏差。综上所述，合理搭建构型，确定体系能量描述并选择合适的构型演变模式，是决定模拟精度的关键。本团队通过模块化的思想提出了MC模拟离子导体材料离子输运特性、相变特性的一般模式。这将大大方便后续开发者用来参考并快速修改，进而扩展MC模拟到具体材料体系。相关模拟的框架主要涵盖以下几个部分(图4)：模拟目标的确定、模拟参数与模型的设定、结构演化方式的选择、转换过程中的能量计算、新结构接受的概率计算以及样本数据的获得。在具体问题中可能还需要增加某些新的模块来处理特殊的需求。基于此流程，本文给出了KMC计算石榴石结构离子导体离子输运特性的一般范式，如图5(a)所示，其中模型搭建依托于BVSE与CAVD的分析结果。模拟过程中的输入输出如表3所示。此外，基于对渗流模拟过程(模型搭建→局域连通性分析→渗流团簇分析→迁移网络统计)的分析，提出了通过KMC模拟搭建渗流模拟初始构型，通过KMC给出的离子迁移模式判断局域通堵，进而对离子输运网络分析的一般范式[图5(b)]。KMC模拟程序框图如图6所示，相应的代码与算例可从代码库https://gitee.com/jinpingliu/pKMC-code.git获得。KMC模拟主程序计算步骤如下。图4   MC模拟的一般步骤，图中列举了MC样本数据的生成过程，蓝色一般为动力学模拟需要的模块，橙色一般为热力学模拟需要的模块，绿色为两者通用模块，灰色是其他模拟中可能存在的模块(此处粒子包含带电离子与其他粒子，具体到离子导体中将用离子描述)
图5   KMC模拟流程与计算案例：(a) KMC对Li₇La₃Zr₂O₁₂(LLZO)的建模与模拟过程示意图(模型构建方式或计算内容依托于具体研究的问题，此处以石榴石结构Li₇La₃Zr₂O₁₂为例，采用BVSE与几何分析方法搭建模型)；(b) KMC模拟与渗流模拟的结合示意图(KMC模拟为渗流模拟提供结构模型、迁移离子排布、离子迁移模式等数据，实现对可达离子，有效载流子浓度等的计算)；(c) 通过本KMC模拟程序计算得到LLZO的占据率变化趋势与实验结果的比较；(d) 通过本KMC模拟程序计算得到LLZO的电导率预测趋势与实验结果的比较

表3   MC模拟常用的输入与输出

图6   MC模拟在离子导体中的计算程序详细流程框图（1）读取外部输入信息，如位点连接关系及位点坐标、迁移能垒、反转(迁移)模式等，根据设定的温度、迁移离子浓度等物理参量搭建初始模拟系统。（2）记录此时系统的状态，记为W₀，计算此时系统的能量，为E₀。（3）对系统进行随机演变，得到新的组态，记为W₁，计算此时系统的能量，记为E₁。（4）对比演变前后组态能量的变化，依据概率P接受新的组态，并更新时间。（5）重复(3)～(4)过程，直到体系能量基本维持不变，开始记录样本数据。（6）改变温度，迁移离子浓度等物理参数，重复上述过程，记录多组样本数据，判断是否结束KMC模拟。（7）对样本数据进行计算，得到迁移率、占据率等物理量。需要补充的是，上述程序的开发过程依托于对固态电解质材料的计算，其中相关的输入参数使用到了本团队前期开发的BVSE与CAVD的相关代码。这些输入参数包括骨架结构的离子输运通道，离子在相邻位点之间迁移的平均势垒等。图5(c)～(d)给出了部分计算与实验的结果对比，由图可知，计算得到的Li_xLa₃Zr₂O₁₂的离子占据情况与电导率变化均与实验结果相符。

3 蒙特卡罗模拟在离子导体中的应用

多数MC模拟均可对应到前文对MC模拟的理论基础与计算范式的介绍与梳理。为了描述MC模拟的应用场景，本节将从离子导体的相变过程、电压平台的变化过程、离子输运特性以及晶体生长动力学特性的分析等方面进行详细的介绍与总结。

3.1　相变分析与电压平台的蒙特卡罗模拟

在锂离子电池中，正极材料的离子输运特性直接决定了电池的容量与循环性能。LiFePO₄作为一种被广泛商用的正极材料，在平衡状态下，会自发地分裂出贫锂和富锂两相。如何快速获取相的分离过程十分重要。Xiao等通过DFT计算得到KMC模拟的能垒输入，搭建了如图7(a)所示的模型。模拟结果显示，LiFePO₄和FePO₄之间形成有序的Li_0.5FePO₄相[图7(b)]。Hin等通过KMC模拟了Li_xFePO₄橄榄石纳米晶体中锂离子脱出与嵌入动力学特性。结果显示，锂离子嵌入过程中电池电压的变化包含以下几个阶段：①锂离子插入贫锂相导致电池电势下降；②富锂相Li_1-βFePO₄成核后，电势几乎恒定；③贫锂相完全演变为富锂相以后，电池的电势继续下降[图7(c)～(d)]。除了对LiFePO₄不同Li离子浓度与电压分布关系的计算，Zheng等还对LiNi_xMn_2-xO₄/Li电池系统Ni掺杂下的电压分布进行了计算。作者首先构建了一个5 × 5 × 5的格子气模型，包含1000个可用位点。在模拟过程中，前500个MC模拟步作为弛豫过程，接下来的500个模拟步用于热力学量的计算，最终得到电压分布与x值之间的依赖关系，如图7(e)所示。可以看到随着x的增大，4.7 V电压平台的长度线性增大，4.1 V电压平台的长度线性减小，而总的电池容量则保持不变。当x = 0.5时，只观察到4.7 V的电压平台。有序到无序的转变和4.1、4.7 V两个电压平台的分裂都可以归因于锂离子与骨架结构之间的复杂相互作用。不同于上述对单个半电池的模拟，Kar等通过对两个半电池(LiMn₂O₄正极和碳负极)进行建模，利用GCMC模拟给出了温度等参数变化对Li⁺嵌入与脱出开路电位、电池电流、电池电压、自由能等电池参数的影响。在模拟过程中，当固定温度T、体积V和化学势μ时，每个位置的能量可由式(5)给出

（5）

式中，J是锂离子之间的相互作用能，下标NN和NNN分别表示最近邻和次近邻，如果位点“i”(“j”)被占据，)取值为1，否则为0。J_NN和J_NNN分别为37.5和25.0 meV。位点占据接受概率可通过式(6)进行描述

（6）

其中，e为电子电荷量。计算得到的开路电位结果如图7(f)所示。这组模拟结果有效地解释了测量开路电位剖面实验中所观察到的迟滞现象。不难发现，上述计算式(5)与(6)来自式(1)与(2)的变形，这种计算模式也被大量的MC模拟所采用。图7   MC模拟对正极材料相变与开路电压的分析：(a) MC模拟模型，上图为Li_xFePO₄中Li的分布；下图为对应的能量图。绿色的球代表了锂离子；(b) 在不考虑电荷间作用能的情况下，Li离子沿着a轴的排布。如图所示，LiFePO₄(LFP)和FePO₄(FP)相之间形成了清晰的相边界，Li_0.5FePO₄相收缩到一个有限的区域；(c) Li_xFePO₄橄榄石纳米晶体在室温下恒电流放电过程的MC模拟。灰点代表活性颗粒中的锂原子，表面电流密度为0.5 A/m²。四组图分别为0 s：初始固溶体，0.0001 s：形成两个富锂相并结合在一起，0.00056 s：富锂相生长，10.81 s：贫锂相几乎完全被消耗；(d) 电池电压与活性材料中Li浓度(mol)的关系；(e) MC模拟计算出的Li/LiNi_xMn_2-xO₄电池电压分布，分别为x = 0.1、x = 0.2、x = 0.3、x = 0.5；(f) 全电池(LiMn₂O₄正极和碳负极)开路电位与放电期间正极材料中Ni占据率的关系除了对正极材料相变与开路电压的模拟，MC模拟同样适用于对石墨基、硅基负极材料相特征的模拟。Moon等利用密度泛函理论(DFT)和KMC模拟，研究了锂离子插入c-Si和a-Si体系的热力学和动力学特征。DFT计算的形成能揭示了晶态硅与非晶硅在锂化过程中的相分离机制。晶态硅和非晶硅在锂化作用下的体积膨胀和相变趋势是相似的，而Li的扩散动力学在c-Si和a-Si之间则存在较大差异。在c-Si中，Li迁移势垒为0.6 eV，随着Li浓度的增加，迁移势垒迅速减小到至0.4 eV[图8(a)～(c)]。为了利用KMC模拟非晶硅中锂的扩散，首先要利用体积函数推导出非晶硅中锂迁移势垒。通过KMC模拟发现锂在a-Si中的扩散系数比在c-Si中的扩散系数大一个数量级[图8(d)]。这些研究有助于在原子尺度理解锂化作用的机理，进一步阐明硅化锂的相分离过程。此外，Persson等也通过DFT、团簇扩展和MC模拟预测了Li_xC₆相图，预测结果在高Li浓度(x>0.5)与实验结果相吻合[图8(e)]。图8   Li离子排布对Li离子迁移势垒的影响[(a)～(c)]；(d) a-Si与c-Si中Li离子的扩散系数；(e) 通过MC模拟获得的Li_xC₆的第一性原理相图。图中符号为：G为石墨，Ⅱ为Ⅱ级相，ⅡD为Ⅱ级无序相，Ⅰ为Ⅰ级相，ⅠD为Ⅰ级无序相，Ⅱ'为Ⅱ级相且Li离子组合为2×2为区分KMC与GCMC模拟在电池负极材料模拟中的差异性，Arriazu等分别通过KMC与GCMC对锂离子嵌入石墨的过程进行了分析。图9(a)所示为Dumas-Hérold锂/石墨插层化合物模型，其中蓝色代表锂离子，灰色代表石墨层，锂离子的嵌入过程呈现出多阶段(台阶)现象。离子的运动模式近似为二维运动[图9(b)]，通过KMC模拟及GCMC模拟离子嵌入脱出随时间的变化[图9(c)]，他们发现Daumas-Hérold模型[图9(e)]中提出的结构是由于动力学限制造成的。如果系统的平衡是通过人为设置一个高交换电流密度(GCMC)的方式来实现，结构就变成了一个具有更低能量值[图9(f)]的经典Rüdorff-Hoffmann模型[图9(d)]。因此，可以通过改变溶剂的性质来增加实验系统中的交换电流密度进而抑制DH结构的产生。此外，使用超小的石墨薄片也有利于离子与电解质的交换，进一步抑制DH结构的产生。图9   MC模拟石墨负极中锂离子的排布：(a) 锂/石墨插层化合物模型；(b) 离子在石墨中的传导模型；(c) Li离子插入、脱出过程中的Li离子含量与时间的关系；(d) Rüdorff- Hoffmann锂/石墨插层化合物模型；(e) Dumas-Hérold锂/石墨插层化合物模型；(f) 两种模型对应的能量值

3.2　离子输运特性的蒙特卡罗模拟

除了对开路电压与相变特性等的模拟，MC也可用于对电极材料的离子扩散特性的分析。Ouyang等利用第一性原理计算得到锂离子和Cr离子在纯LiFePO₄和掺Cr的LiFePO₄中沿一维扩散路径迁移的迁移能垒[图10(a)]，采用MC模拟研究了势垒对锂离子二次电池正极材料LiFePO₄电化学性能的影响。模拟中，Cr离子被证明不能在晶体中迁移，基于此，作者搭建了如图10(b)所示的一维离子输运模型。模拟结果显示，随着Cr掺杂量的增加，薄膜的容量不断减小[图10(c)]，随着超胞尺寸的增大，所有掺杂量下的容量也都大大降低。这说明减小粉末正极材料的粒径可以大大提高其容量。图10   (a) Cr离子与Li离子的迁移能垒示意图；(b) 在每条一维通道中，锂离子被Cr离子隔开的示意图；(c) 对于不同的模拟超胞尺寸，容量与掺杂量之间的依赖关系。数字200、300、500、1000和2000指的是超胞的大小
对于固态电解质，其离子输运的动力学特性更加受到关注，KMC在其中的应用更为广泛。由于典型的固态电解质存在三大特点：①相邻位点之间的离子迁移容易发生；②可供迁移离子占据的位点的数量大于迁移离子的数量；③这些可用的位点连接成一个连续的扩散通路。这些特点使得格子气KMC可以很好地适用于固态电解质材料的模拟。Morgan等根据石榴石构型固态电解质的结构特征[图11(a)～(c)]，采用晶格气体MC模拟预测了石榴石结构固态电解质Li_xLa₃Zr₂O₁₂的相关系数变化趋势，发现当x_Li = 3(来自位置能量)和x_Li = 6(来自最近邻斥力)时，其存在特别强的相关效应[图11(d)～(e)]。图11   石榴石框架中离子迁移扩散网络示意图与相关系数随迁移离子浓度的变化：(a) 石榴石结构中由ZrO₆八面体和LaO₈十二面体连接的离子输运通道的三维图。迁移离子随机分布在四面体和八面体间隙中。相互连通的四面体隙和八面体间隙形成三维离子输运网络；(b) 石榴石结构中环形结构示意图，其中四面体和八面体空位被部分锂离子占据。应该注意的是，并不是所有的空位都需要被锂离子占据；(c) 石榴石离子输运通道的二维几何连接。一个八面体间隙与两个四面体间隙相连，一个四面体间隙与四个八面体间隙相连(四个八面体空位与四面体的每个面相连，图中仅展示两个八面体做示例)。在(a)、(b)和(c)中用蓝色表示八面体空位，用黄色表示四面体空位；(d) 集体相关系数随着Li离子浓度(x_Li)的变化；(e) 自相关系数随着Li离子浓度(x_Li)的变化，x_Li指的是Li_xLa₃Zr₂O₁₂中Li离子的化学计量数，取值0～9考虑到键价和理论等静态计算方法的局限性，Chen等提出将键价和理论计算的迁移能垒数据与KMC模拟相结合的方法，得到了指定温度下的绝对电导率[图12(b)]，这显然优于仅考虑几何效应计算得到的电导率变化趋势[图12(a)]。由于迁移离子的具体排布对离子迁移能垒具有极大的影响，所以上述的近似也不够精确。为此，Ven等基于第一性原理团簇展开的方式，对KMC模拟的势垒参数进行了计算。他们发现在无序材料中，应用团簇展开可以近似获得任何构型下的迁移势垒。这意味着在KMC模拟中，任何时刻都能近似地保持细致平衡，并能在热力学平衡状态下计算出所需的动力学量。作者进一步给出了这一研究范式在Li_xCoO₂中的应用。结果表明，迁移机制和激活能垒在很大程度上取决于迁移锂离子周围近邻锂离子与锂空位的排布[图12(c)]。通过模拟，作者发现在所有锂浓度下，锂在层状Li_xCoO₂中的扩散依赖于双空位机制。此外，由于激活能垒受离子浓度影响巨大，扩散系数随锂离子浓度x的变化可以达到几个数量级。综上所述，KMC的模拟一般依赖于其他计算方法的输入(包括模型结构，迁移能垒等)，这对计算结果精度的影响是巨大的。图12   几个KMC计算过程的对比：(a) LLZO的通道连接关系与仅考虑几何效应下的离子电导率示意图；(b) BVSE计算得到的离子输运通道与通过KMC模拟得到的电导率与实验值的比较。实心蓝色曲线和填充圆表示KMC所模拟的电导率。实红线表示有效计算的电导率(来自高温KMC模拟)，虚线表示对较低温度下离子电导率的推断。三角形绿色曲线表示实验电导率数据；(c) 几种典型的离子排布状态与基于局域结构团簇展开的KMC计算得到的离子扩散系数考虑到KMC模拟可以有效地预测和分析晶体材料离子电导率，Grieshammer等开发了一款名为MOCASSIN的固态离子导体MC计算软件，利用该软件对富镧橄榄石中的氧间隙迁移以及掺杂全水合锆酸钡中的质子传导进行了模拟研究。模拟结果揭示了缺陷相互作用对离子电导率的影响。这些效应的组合可能会导致固态离子材料中意想不到的离子输运行为，特别是在具有较多可移动离子的模拟体系中。除此之外，Hoffmann等也开发了一款名为Kmos的软件，在软件描述的系统中，所有的反应空间都可以抽象为一系列离散的空间格点(即为前文的格子气模型)。Kmos通过代码生成新的可用事件来不断演化，以实现较高的模拟效率。

3.3　晶体生长及电解质界面膜(SEI膜)生成的蒙特卡罗模拟

除上述对离子导体材料的模拟之外，MC模拟在晶体生长以及SEI膜层生长等领域也有着重要的应用。基于晶体生长理论的计算模拟是深入研究晶体生长机理，进一步指导晶体生长实验的一项重要技术。MC在模拟晶体生长的过程中分别考虑了原子的沉积、表面扩散、吸附等过程。在KMC模拟中，晶体的生长过程可被视为一个随机过程，使用某些特定的生长模型来描述其中随机事件的发生概率。通常具体模拟步骤如下：①建立微观模型；②对模型进行描述；③模型发生演变；④比较演变前后的结果确定演变概率，这与第二部分的相关总结类似。图13(a)～(g)展示了简立方晶格模型的表面在不同配位下的晶体生长过程，可以看出这几种配位对晶体长大的速度与方向均会产生极大的影响。Barai等基于简化的二维格子气模型通过KMC模拟建立了一个多尺度计算模型，用以捕捉晶体初生离子的形核和生长以及它们聚集成二次过渡金属氢氧化物前驱体离子。结果显示随着溶液pH的增加和氨浓度的降低，二次颗粒的尺寸减小。结合模拟给出的相图[图13(h)]，相关研究可以为实验合成提供有效指导。Chen等采用KMC模拟研究了4H和6H型SiC的竞争生长过程。该模型以硅和碳原子为基本迁移离子，构建了格子气KMC模型。该模拟不仅考虑了同种原子之间的相互作用能(Si-Si)，还考虑了不同类型原子之间的相互作用能(Si-C)。基于此模型，作者进一步研究了生长条件对4H-SiC和6H-SiC竞争生长的影响。结果表明，在生长过程中较高的温度和高比例的Si源会促进6H相的生长。同时，随着温度和Si/C比的升高，晶体表面原子聚集生长的现象越来越明显[图13(i)～(j)]。KMC模拟SEI的生成过程一般分为如下几个过程：吸附、解吸、表面扩散以及钝化(非活性)材料的形成[图13(k)]。基于此流程，Hao等利用KMC研究了充电期间石墨电极上SEI的生长过程。结果显示，锂离子的嵌入速率和充电时间均受锂还原速率和通过SEI的锂传输速率的控制。在较低的温度下，Li的还原速率较高，但SEI的电阻也较高，所以会限制扩散过程。另一方面，在较高温度下，较低的Li还原速率会导致总充电时间增加[图13(l)]。图13   (a)～(g)为简单立方晶格表面不同配位下的晶体生长过程；(h) 晶体尺寸与溶液pH值与氨含量依赖关系的相图。图中显示了在粒径分布为高斯分布的情况下，溶液的平均粒径及其对应的一阶标准差。这里显示的所有数据都是从计算模型中获取的。绿色区域表示次级颗粒较小，黄色区域表示次级颗粒较大。浅蓝点的二次粒径分布标准差较小，紫色点的二次粒径分布标准差较大。在pH值和氨含量方面，沉积较大(Dpart约8 μm)和较小(Dpart约4 μm)的二次活性颗粒的最佳操作条件也在图中突出显示(用红色圆圈表示) ；(i) T = 1100 K下的不同Si/C比条件下，两种晶体在沉积速率F = 0.1 ML/s下，6H-SiC在生长中的比例；(j) T = 1400 K下的不同Si/C比条件下，两种晶体在沉积速率F = 0.1 ML/s下，6H-SiC在生长中的比例；(k) SEI膜的生成过程，一般分为吸附、解吸、表面扩散和钝化等步骤；(l)SEI膜厚度、充电速度与温度的依赖关系目前对晶体生长的MC模拟仍存在许多问题，比如Maazi等模拟晶体生长时，在随机矩阵位点上需要重复数百万次晶粒生长过程。这对于较大的晶格系统来说计算时间会很长。为了提高晶粒生长模拟的计算速度，研究人员对经典MC模型进行了两种修正：①所有阵点都以相同的概率进行随机重定向，即每个蒙特卡罗步都不重复；②位置能量除以一个参数，该参数表示离子存在时，作用在晶界上的有效作用。这些举措极大地加速了大颗粒的生长，减弱了邻近小颗粒的生长，有效地提高MC模拟的速度。

4 结论与展望

本文总结了MC模拟的发展历史，归纳出了一套适用于多数MC模拟的范式，包括以下几个部分：模拟目标的确定、模拟参数与模拟模型的设定、结构演化方式的选择、转换过程中判据的提取(如体系能量)、新结构接受概率的计算、样本数据的获得等。根据此范式，列出了本团队前期开发的一套典型MC模拟代码供读者参考使用，并给出了一个典型的计算案例：预测石榴石结构固态电解质的有效载流子浓度、迁移离子分布与迁移离子浓度的依赖关系，模拟结果与实验规律基本一致。除此之外，本文就其在正负极材料、电解质以及相关界面中的典型热力学与动力学计算案例进行了剖析，包括求解离子导体中的离子扩散问题、模拟迁移离子的分布特征与相关界面层的生长过程等。这些计算不仅涵盖了实验现象的微观机理解释，也包含了实验结果的预测进而提供具体的实验方案指导。在与实验结合方面，MC模拟可以基于某些实验数据(如迁移离子的分布)搭建初始模型，提出针对具体研究体系的改进策略。如Huang等基于实验观测，提出了层状金属氧化物中过渡金属在八面体位点与四面体位点的迁移会改变迁移离子的输运特性。鉴于每种模拟方法本身的局限性，亦可以采用多方法联用的方式，扩展MC模拟的输出，比如文中提到的KMC与渗流的结合以获得有效载流子浓度。为了促进MC模拟在离子导体计算中的发展，本文总结了以下几个MC方法的发展方向，包括：①在MC模拟中如何快速精确地获取所有可能发生的事件以及对应的转变概率。在材料计算领域，这些事件及其描述主要是通过实验测试以及第一性原理计算的方式获得，比如通过密度泛函理论(DFT)计算迁移能垒，通过分子动力学预判离子的迁移模式等。然而这些参数的计算需要解决微观演变中所涉及到的各种近似，比如基于第一性原理模拟的团簇展开就是对能量计算的一种简化近似。除此之外，体系转换速率因子也直接决定了模拟体系的演变速度，该参数一般需要近似为某一个常数，但对于一些复杂体系(比如存在多种迁移模式的晶态离子导体)，转换因子(如离子的迁移尝试频率等)是极为复杂的，虽然通过声子谱等手段可以进行较为精细的计算，但对于较大体系或者较长时间的模拟，相关的计算成本是巨大的。②晶格气体MC模拟过程中仅考虑了当前对模拟体系有利的转移路径，这将导致模拟会很快进入某一亚稳态。当前常用的Metropolis算法虽然能够一定程度上改善此问题，但对于存在多个亚稳态的模拟体系，以上算法还是很难精确地寻找到系统的演化轨迹，并且其计算成本也极为高昂。③MC模拟过程中，模拟事件之间存在先后顺序，这不仅制约了其计算效率也影响了其模拟准确性。比如在极短时间内离子的迁移存在一定的先后顺序，如何确定迁移的先后顺序会直接影响离子迁移的成功率。④MC模拟始终在“自己的”时间域内进行，与真实时间存在差异。虽然可以通过数学手段对MC时间(也叫做蒙特卡罗步)与真实时间进行转换，但这种转换仅为统计上的一种近似，依赖于众多初始输入参数(如尝试频率等)。

相关文章：