储能科学与技术, 2022, 11(3): 878-896 doi: 10.19799/j.cnki.2095-4239.2022.0050

储能新材料设计与先进表征专刊

基于蒙特卡罗模拟的离子导体热力学与动力学特性

刘金平,1, 蒲博伟1, 邹喆乂2, 李铭清1, 丁昱清1, 任元1, 罗亚桥1, 李杰3, 李亚捷1, 王达1, 何冰4, 施思齐,1,5,6

1.上海大学材料科学与工程学院，上海 200444

2.湘潭大学材料科学与工程学院，湖南 湘潭 411105

3.加州大学欧文分校物理与天文系，美国 加州 92697

4.上海大学计算机工程与科学 学院

5.上海大学材料基因组工程研究院，上海 200444

6.之江实验室，浙江 杭州 311100

Investigating thermodynamic and kinetic properties of ionic conductors via Monte Carlo simulation

LIU Jinping,1, PU Bowei1, ZOU Zheyi2, LI Mingqing1, DING Yuqing1, REN Yuan1, LUO Yaqiao1, LI Jie3, LI Yajie1, WANG Da1, HE Bing4, SHI Siqi,1,5,6

1.School of Materials Science and Engineering, Shanghai University, Shanghai 200444, China

2.School of Materials Science and Engineering, Xiangtan University, Xiangtan 411105, Hunan, China

3.Department of Physics and Astronomy, University of California, Irvine, California 92697, USA

4.School of Computer Engineering and Science, Shanghai University

5.Materials Genome Institute, Shanghai University, Shanghai 200444, China

6.Zhejiang Laboratory, Hangzhou 311100, Zhejiang, China

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

收稿日期: 2022-01-28   修回日期: 2022-02-12  

基金资助: 国家自然科学基金面上项目.  11874254
国家重点研发计划项目.  2021YFB3802104
上海先进陶瓷结构设计与精密制造专业技术服务平台.  20DZ2294000
之江实验室科研攻关项目.  2021PE0AC02

Received: 2022-01-28   Revised: 2022-02-12  

作者简介 About authors

刘金平（1994—），男，硕士研究生，从事超离子导体离子输运特性的研究，E-mail：506094182@qq.cm； E-mail：506094182@qq.cm。

摘要

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

关键词： 蒙特卡罗模拟 ; 离子导体 ; 离子输运特性 ; 相变特性 ; 晶体生长

Abstract

Monte Carlo (MC) simulation, on the basis of probability and statistics theory, was proposed by Von Neumann et al. in 1940s. As an important numerical method, MC simulation has been used to investigate thermodynamic and kinetic properties of ionic conductors. However, there exists a large improvable space for the MC simulation in the calculation accuracy, simulation efficiency and simulation process automation. In this work, through systematic analysis of the Hamiltonian model in MC simulation (For example, based on bond-valence theory or cluster expansion, the configuration energy is given by fitting the neighbors interaction parameters that can always be obtained by first principles calculations of representative small supercells) and the evolution model (For example, configuration evolution based on the assumption of single-ion jump mode) of material structure, a set of MC simulation paradigms for analyzing the ion transport and phase transition characteristics of ionic conductors are extracted, and the corresponding semi-automatic simulation codes are given which can be used for predicting the respective dependences of the ionic conductivity and the occupancy of the migrated ions in the garnet-structured ionic conductor with the lithium ion concentration. To expand the application of MC simulation in the research of ionic conductor, we further analyze its applications in the typical thermodynamic and kinetic properties calculations of electrochemical energy storage materials which includes anode and cathode materials, electrolytes and the related interfaces, including the ionic diffusion problems, the distribution characteristics of migrated ions and the evolution of related interface. Finally, the current challenges faced by MC methods are prospected and the possible solutions are presented, including: ① Accurately capturing all possible events (such as single-ion hop and two-ion cooperation hop) and their descriptions (such as the Hamiltonian calculations); ② Searching for efficient algorithm to accurately find the evolution trajectory of the system; ③ Accurately obtaining the corresponding actual time in the MC simulation.

Keywords： Monte Carlo simulation ; ionic conductors ; ion transport properties ; phase transition properties ; crystal growth

PDF (11852KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

刘金平, 蒲博伟, 邹喆乂, 李铭清, 丁昱清, 任元, 罗亚桥, 李杰, 李亚捷, 王达, 何冰, 施思齐. 基于蒙特卡罗模拟的离子导体热力学与动力学特性[J]. 储能科学与技术, 2022, 11(3): 878-896

LIU Jinping. Investigating thermodynamic and kinetic properties of ionic conductors via Monte Carlo simulation[J]. Energy Storage Science and Technology, 2022, 11(3): 878-896

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

表1   MC模拟与MD模拟的对比[18-20, 28-32]

Table 1  Comparison between MC simulation and MD simulation[18-20, 28-32]

对比项目MD模拟MC模拟
第一原理MD模拟平均场MD模拟GCMC模拟KMC模拟
目的模拟系统/结构的演变
初衷为了计算统计力学里的系综平均
统计原理时间的平均系综的平均
方法精确地计算出下一时刻离子的分布“试”出每个MC步后相对合理的离子分布，或者计算出每种构型出现的期望
具体操作计算出某个时刻离子的运动状态。选定最小时间单元统计，得到样本数据通过计算体系的转变概率或者离子的迁移概率，更新体系的状态，有时也会根据转换概率值对应到真实时间的变化
模拟时间尺度飞秒～皮秒飞秒～纳秒皮秒～小时
模拟对象尺度微观微观～介观微观～宏观
模拟关注点原子运动轨迹体系组态迁移
应用范围离子迁移，晶体生长，缺陷生成，材料老化等

新窗口打开| 下载CSV


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

1 蒙特卡罗模拟方法基础

1.1　蒙特卡罗模拟方法

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

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   蒙特卡罗模拟的发展

Fig. 1   The development of Monte Carlo simulation


图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   MC模拟中常见的几种二维格子气模型：(a) 二维正方格子气模型；(b) 二维六边格子气模型；(c) 二维随机格子气模型(紫色的球代表格点上的离子)

Fig. 2   Several two-dimensional lattice-gas model which used in MC simulation commonly: (a) two-dimensional square lattice-gas model; (b) two-dimensional hexagonal lattice-gas model; (c) random lattice-gas model (The purple spheres represent the ions at the lattice sites)


图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   MC模拟中的几种能量近似计算：(a) 石墨结构中的锂离子与其他离子之间的相互作用能的计算[45]；(b) 迁移离子的迁移过程中所需要克服的迁移能垒示意图[46]；(c) 迁移离子的不同迁移模式；(d) 不同迁移模式对应的迁移能垒示意图[15]；(e) 几种蒙特卡罗模拟之间的关系

Fig. 3   Approximation of the energy in MC simulation: (a) calculation of the interaction energy between lithium ions and other ions in graphite structure; (b) schematic diagram of the transition barrier during the jump of the migrated ions; (c) different jump modes of migrated ions; (d) different transition barrier of different jump modes; (e) relation diagram between some Monte Carlo simulations


1.2　蒙特卡罗模拟过程

GCMC[42-43]模拟作为一种适用于巨正则系综限制下经典的热力学MC模拟，被广泛应用于电化学储能材料的合成、多孔吸附等过程的研究[44]。在电化学储能材料计算领域，GCMC常用于模拟石墨等碳基嵌锂材料的结构特征。在模拟过程中，随机选择几种基本的离子运动模式(比如置换、插入、脱出)，并对相应离子运动模式发生的概率进行计算[45]，一般可以通过式(2)所示的计算模式求得。值得指出的是，式中的能量变化不仅可以使用体系的始末态能量表示，也可以使用过渡态能量表示[46]。对于忽略动力学特性的热力学模拟，仅考虑简单的始末态能量是常见的。作为GCMC模拟适用情形的补充，KMC模拟则可以很好地描述离子导体结构演变的动力学过程。KMC模拟的研究最早可以追溯到20世纪60年代对辐射损伤退火问题的模拟。目前，KMC在表面吸附、离子扩散、晶体生长等领域得到了迅速发展和应用[28]。由于模拟机理不同，KMC相对于分子动力学(MD)模拟可以对更大系统的动态演化过程给出更长时间的模拟[34, 47]。离子导体中迁移离子的成功迁移并不频繁，其状态转变会经历长时间的“尝试”，只有极少的“尝试”会成功。KMC模拟不考虑尝试过程只关注可以演变成功的状态，这种状态到状态的转变不受原子振动的限制，不局限于对每个原子(离子)的计算，因此可以实现更长时间尺度的模拟[48]。同时，KMC的离散化思想更适合对微观离散体系进行模拟分析，这也使得离散化模型的搭建在KMC模拟中占据着十分重要的地位[49]。LGMC[32, 50]模拟即为一种典型的离散化处理方式。

GCMC模拟与KMC模拟均可以借鉴格子气模型而发展成为一种特殊的LGMC模拟来描述离子在材料体系中的扩散与占位情况。在格子气模型的搭建过程中，时间与空间被分割，按照离子之间的相互作用以及迁移、翻转规则，将模拟体系简化成一系列固定的空间格点与可迁移或翻转的离子。图2展示了几种典型格子气模型的示意图。该模型中，离子位于格点上，沿着边运动，离子与离子之间存在相互作用(比如吸引或者排斥作用)。Ising模型[18]、渗流模型[51]等均是对格子气模型的某种处理与变形。

在研究离子导体材料时，MC模拟通常使用体系的能量来确定体系的演变过程[52-54]。体系能量则通过描述近邻迁移离子间以及迁移离子与骨架离子之间的相互作用得到，这样既可以简化能量的计算过程也可兼顾描述的完整性。常用的能量计算公式如式(1)所示[20, 55]

E=∑j=1mEsite_j×cj+ 12∑j=1m∑i=1i≠jmEnn×cij

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

Pa→b∝ 1,                   if ∆Ea→b<0e(-∆Ea→bkT),     if ∆Ea→b≥0

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

∆Ea→b=Eb-Ea

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

∆t = 1Kij ln μ

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

表2   MC模拟中的常用概念

Table 2  Common concepts in MC simulation

名称释义名称释义
骨架结构离子导体中位于原位振动而不会长程迁移的离子构成的结构相互作用能包括迁移离子之间以及迁移离子与骨架结构之间存在的能量。广义上还应包括外场的作用，包括电场，磁场等对目标离子的作用能
迁移路径离子移动的路径，可以通过BVSE等方法求得。哈密顿量模拟中体系的能量表示形式，依托于结构特性或者模拟对象进行搭建，一般为各相互作用能的和
反转模式一般在磁性材料的模拟中，会选择某些位点磁性的反转，涉及一维、二维、三维等情形蒙特卡罗步也称为一个蒙特卡罗遍历，将所有的操作单元进行一次随机尝试
迁移能垒离子移动过程中需要克服的最大能量值与初态能量值之差弛豫过程体系由高能量状态逐渐过渡到某一低能量状态的过程
反转能垒离子反转过程中需要克服的最大能量值与初态能量值之差接受概率新体系被接受的概率，一般可以通过Metropolis算法等获得
嵌入能垒离子嵌入过程中需要克服的最大能量值与初态能量值之差

新窗口打开| 下载CSV


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

MC模拟已被应用到众多领域，但是如何提高其模拟精度与效率一直是众多科研工作者的重要研究方向[57]。为了清晰了解蒙特卡罗模拟的计算过程，此章节就其计算流程和步骤进行了总结与提炼。一般来说，模型的搭建应尽可能地复现材料的结构特征，比如通过几何分析方式或BVSE来分析材料框架以及离子传输路径，通过格子气模型分析常见的晶态离子导体，通过Ising模型来分析磁性材料等。此外，合适的模型描述也是决定相关模拟精度的关键。通常研究者会通过对模型能量(如内能)的计算，进而以能量的变化来确定演变方向以及对应构型出现的概率[56-59]。然而，目前还没有一个完美的计算方案可以准确地计算材料构型转变过程中的能量变化。为此，众多科研工作者作了许多努力。比如Arriazu等[45]定义了如图3(a)所示的晶格气体模型，考虑了层内Li离子、层间Li离子以及Li离子与骨架结构的相互作用能，利用GCMC模拟了Li离子在石墨中的插层过程。在模拟中，作者采用一种纯静态的能量计算模式，即只根据系统的总能来确定此构型出现的概率。因其不考虑过渡态能量，所以该模拟方法在动力学模拟中是不够精确的。对于固定的迁移或反转，激活能垒更能准确反映其能量的变化，因此Ven等[46]提出了如图3(b)所示的离子迁移能垒。为避免因离子往返迁移能垒不同而造成的计算困难，其引入了一个动力学解析的激活势垒(∆EKRA)作为能垒的近似，同时还考虑了目标离子周围离子的排布，采用团簇展开的策略，近似计算了离子的迁移概率。然而该近似不仅忽略了离子迁移的方向性，也忽略了多离子迁移模式。迁移模式的选择(材料结构演变途径)是MC模拟需要解决的第三个关键问题。He等[15]证明了超离子导体中的离子快速扩散并不像典型固体中那样通过孤立的离子迁移发生，而是通过具有低能量势垒的多离子协同迁移进行[图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

新窗口打开| 下载原图ZIP| 生成PPT

图4   MC模拟的一般步骤，图中列举了MC样本数据的生成过程，蓝色一般为动力学模拟需要的模块，橙色一般为热力学模拟需要的模块，绿色为两者通用模块，灰色是其他模拟中可能存在的模块(此处粒子包含带电离子与其他粒子，具体到离子导体中将用离子描述)

Fig. 4   General steps of MC simulation, in which listed the calculation process of MC sample data. Blue box is the module required by dynamic simulation, orange box is the module required by thermodynamic simulation, green box is the common module of the both, and gray box represent modules that exist in other simulation(Here, particles include charged ions and other particles, which will be described as ions in ionic conductors)


图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   KMC模拟流程与计算案例：(a) KMC对Li7La3Zr2O12(LLZO)的建模与模拟过程示意图(模型构建方式或计算内容依托于具体研究的问题，此处以石榴石结构Li7La3Zr2O12 为例，采用BVSE与几何分析方法搭建模型)；(b) KMC模拟与渗流模拟的结合示意图(KMC模拟为渗流模拟提供结构模型、迁移离子排布、离子迁移模式等数据，实现对可达离子，有效载流子浓度等的计算)；(c) 通过本KMC模拟程序计算得到LLZO的占据率变化趋势与实验结果[61]的比较；(d) 通过本KMC模拟程序计算得到LLZO的电导率预测趋势与实验结果[61]的比较

Fig. 5   KMC simulation process and calculation cases: (a) schematic diagram of modeling and simulation process of Li7La3Zr2O12(LLZO) by KMC (The model construction method or calculation content depends on specific research questions. Here, the garnet structure Li7La3Zr2O12 is taken as an example, and the model is constructed by BVSE and CAVD method); (b) schematic diagram of KMC simulation combined with percolation simulation (KMC simulation provides the structure model, migrated ion arrangement and ion hop mode for percolation simulation, and percolation simulation is used to calculate accessible ions and effective carriers concentration); (c) the comparison of occupancy calculated by this KMC program with the experimental results; (d) the comparison of ionic conductivity calculated by this KMC program with the experimental results of LLZO


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

Table 3  Some inputs and outputs of MC simulation

输入直接输出间接输出
结构文件(包含骨架与迁移离子信息的CIF文件或者POSCAR文件)样本数据，包含：迁移离子分布，迁移离子状态(比如自旋取向)，以及其对应的时间，温度，电场等信息占据率等离子分布数据
离子迁移/反转/嵌入等模式磁矩、内能、比热容等热力学数据
对应的迁移/反转/嵌入能垒扩散系数、相关系数、电导率等动力学数据
弛豫步数开路电压，中间相等其他数据信息
平衡采样步数
温度，电场，磁场，离子浓度等

新窗口打开| 下载CSV


图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   MC模拟在离子导体中的计算程序详细流程框图

Fig. 6   Flow diagram of MC calculation program in ionic conductor


（1）读取外部输入信息，如位点连接关系及位点坐标、迁移能垒、反转(迁移)模式等，根据设定的温度、迁移离子浓度等物理参量搭建初始模拟系统。

（2）记录此时系统的状态，记为W0，计算此时系统的能量，为E0。

（3）对系统进行随机演变，得到新的组态，记为W1，计算此时系统的能量，记为E1。

（4）对比演变前后组态能量的变化，依据概率P接受新的组态，并更新时间。

（5）重复(3)～(4)过程，直到体系能量基本维持不变，开始记录样本数据。

（6）改变温度，迁移离子浓度等物理参数，重复上述过程，记录多组样本数据，判断是否结束KMC模拟。

（7）对样本数据进行计算，得到迁移率、占据率等物理量。

需要补充的是，上述程序的开发过程依托于对固态电解质材料的计算，其中相关的输入参数使用到了本团队前期开发的BVSE[60]与CAVD的相关代码[6]。这些输入参数包括骨架结构的离子输运通道，离子在相邻位点之间迁移的平均势垒等。图5(c)～(d)给出了部分计算与实验的结果对比，由图可知，计算得到的Li x La3Zr2O12的离子占据情况与电导率变化均与实验结果相符。

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

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

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

在锂离子电池中，正极材料的离子输运特性直接决定了电池的容量与循环性能[62]。LiFePO4[63]作为一种被广泛商用的正极材料，在平衡状态下，会自发地分裂出贫锂和富锂两相。如何快速获取相的分离过程十分重要。Xiao等[59]通过DFT计算得到KMC模拟的能垒输入，搭建了如图7(a)所示的模型。模拟结果显示，LiFePO4和FePO4之间形成有序的Li0.5FePO4相[图7(b)]。Hin等[64]通过KMC模拟了Li x FePO4橄榄石纳米晶体中锂离子脱出与嵌入动力学特性。结果显示，锂离子嵌入过程中电池电压的变化包含以下几个阶段：①锂离子插入贫锂相导致电池电势下降；②富锂相Li1-β FePO4成核后，电势几乎恒定；③贫锂相完全演变为富锂相以后，电池的电势继续下降[图7(c)～(d)]。除了对LiFePO4不同Li离子浓度与电压分布关系的计算，Zheng等[65]还对LiNi x Mn2-x O4/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等[19]通过对两个半电池(LiMn2O4正极和碳负极)进行建模，利用GCMC模拟给出了温度等参数变化对Li+嵌入与脱出开路电位、电池电流、电池电压、自由能等电池参数的影响。在模拟过程中，当固定温度T、体积V和化学势μ时，每个位置的能量可由式(5)给出

εiT,V,μ=ciJNN∑NNcj+JNNN∑NNNcj-μ

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

Pi=exp-εf-εie∕kBT

其中，e为电子电荷量。计算得到的开路电位结果如图7(f)所示。这组模拟结果有效地解释了测量开路电位剖面实验中所观察到的迟滞现象。不难发现，上述计算式(5)与(6)来自式(1)与(2)的变形，这种计算模式也被大量的MC模拟所采用。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7   MC模拟对正极材料相变与开路电压的分析：(a)MC模拟模型，上图为Li x FePO4 中Li的分布；下图为对应的能量图。绿色的球代表了锂离子[59]；(b) 在不考虑电荷间作用能的情况下，Li离子沿着 a 轴的排布。如图所示，LiFePO4(LFP)和FePO4(FP)相之间形成了清晰的相边界，Li0.5FePO4 相收缩到一个有限的区域[59]；(c) Li x FePO4 橄榄石纳米晶体在室温下恒电流放电过程的MC模拟。灰点代表活性颗粒中的锂原子，表面电流密度为0.5 A/m2。四组图分别为0 s：初始固溶体，0.0001 s：形成两个富锂相并结合在一起，0.00056 s：富锂相生长，10.81 s：贫锂相几乎完全被消耗；(d) 电池电压与活性材料中Li浓度(mol)的关系[64]；(e) MC模拟计算出的Li/LiNi x Mn2-x O4 电池电压分布，分别为 x = 0.1、 x = 0.2、 x = 0.3、 x = 0.5[65]；(f) 全电池(LiMn2O4 正极和碳负极)开路电位与放电期间正极材料中Ni占据率的关系[19]

Fig. 7   MC simulation for phase transition and open circuit voltage of positive material: (a) MC simulation model (The above figure shows the distribution of Li ions in Li x FePO4, the corresponding energy diagram is shown below, the green balls represent Li ions); (b) the arrangement of Li ions along axis a without considering the interaction energy between charges. As shown in the figure, a clear phase boundary is formed between the LiFePO4 and FePO4 phases, and the Li0.5FePO4 phase shrinks to a limited region; (c) MC simulation of constant current discharge process of Li x FePO4 olivine nanocrystals at room temperature. The gray dots represent lithium atoms in the active particles with a surface current density of 0.5 A·m-2. The four groups are 0 s: initial solid solution, 0.0001 s: forming two Li-rich phases and combining together, 0.00056 s:Li-rich phase growth, 10.81 s: the Li-poor phase was almost completely consumed; (d) the relationship between the battery voltage and the concentration of Li ions in the active material; (e) the voltage distribution of Li/LiNi x Mn2-x O4 battery calculated by MC, x = 0.1, x = 0.2, x = 0.3, x = 0.5, respectively; (f) the relationship between the open circuit potential of a full battery (LiMn2O4 anode and carbon cathode) and the Ni occupation of the cathode during discharge


除了对正极材料相变与开路电压的模拟，MC模拟同样适用于对石墨基、硅基负极材料相特征的模拟。Moon等[11]利用密度泛函理论(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等[66]也通过DFT、团簇扩展和MC模拟预测了Li x C6相图，预测结果在高Li浓度(x>0.5)与实验结果相吻合[图8(e)]。

图8

新窗口打开| 下载原图ZIP| 生成PPT

图8   Li离子排布对Li离子迁移势垒的影响[(a)～(c)][11]；(d)a-Si与c-Si中Li离子的扩散系数[11]；(e) 通过MC模拟获得的Li x C6 的第一性原理相图。图中符号为：G为石墨，Ⅱ为Ⅱ级相，ⅡD为Ⅱ级无序相，Ⅰ为Ⅰ级相，ⅠD为Ⅰ级无序相，Ⅱ'为Ⅱ级相且Li离子组合为2×2[66]

Fig. 8   Effect of Li ions arrangement on the transition barrier [(a)-(c)]; (d) diffusion coefficient of Li ion in a-Si and c-Si; (e) first-principles phase diagram obtained from Monte Carlo simulations. The phase regions are denoted in the following way: G (graphite), Ⅱ (stage Ⅱ), ⅡD (disordered stage Ⅱ), Ⅰ (stage Ⅰ), ⅠD (disordered stage Ⅰ), and Ⅱ′ (stage Ⅱ with 2×2 Li ordering)


为区分KMC与GCMC模拟在电池负极材料模拟中的差异性，Arriazu等[45]分别通过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

新窗口打开| 下载原图ZIP| 生成PPT

图9   MC模拟石墨负极中锂离子的排布：(a) 锂/石墨插层化合物模型；(b) 离子在石墨中的传导模型；(c)Li离子插入、脱出过程中的Li离子含量与时间的关系；(d) Rüdorff- Hoffmann锂/石墨插层化合物模型；(e) Dumas-Hérold锂/石墨插层化合物模型；(f) 两种模型对应的能量值[45]

Fig. 9   Distribution of Li ions by MC simulation in graphite anode: (a) lithium/graphite intercalation compound model; (b) ionic conduction model in graphite; (c) the relationship between Li ion content and time in the process of Li ion insertion and removal; (d) Rüdorff-Hoffmann model of lithium/graphite intercalation compound; (e) dumas-Herold lithium/graphite intercalation compound model; (f) changes in energy values corresponding to the two models


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

除了对开路电压与相变特性等的模拟，MC也可用于对电极材料的离子扩散特性的分析。Ouyang等[67]利用第一性原理计算得到锂离子和Cr离子在纯LiFePO4和掺Cr的LiFePO4中沿一维扩散路径迁移的迁移能垒[图10(a)]，采用MC模拟研究了势垒对锂离子二次电池正极材料LiFePO4电化学性能的影响。模拟中，Cr离子被证明不能在晶体中迁移，基于此，作者搭建了如图10(b)所示的一维离子输运模型。模拟结果显示，随着Cr掺杂量的增加，薄膜的容量不断减小[图10(c)]，随着超胞尺寸的增大，所有掺杂量下的容量也都大大降低。这说明减小粉末正极材料的粒径可以大大提高其容量。

图10

新窗口打开| 下载原图ZIP| 生成PPT

图10   (a) Cr离子与Li离子的迁移能垒示意图；(b) 在每条一维通道中，锂离子被Cr离子隔开的示意图；(c) 对于不同的模拟超胞尺寸，容量与掺杂量之间的依赖关系。数字200、300、500、1000和2000指的是超胞的大小[67]

Fig. 10   (a) schematic diagram of transition barrier between Cr ion and Li ions; (b) schematic diagram of lithium ions separated by Cr ions in each one-dimensional channel; (c) variation relation of the simulated capacity with the dopant amount and super-cell size. The numbers 200, 300, 500, 1000 and 2000 refer to the size of the super-cell being used


对于固态电解质，其离子输运的动力学特性更加受到关注，KMC在其中的应用更为广泛。由于典型的固态电解质存在三大特点[68]：①相邻位点之间的离子迁移容易发生；②可供迁移离子占据的位点的数量大于迁移离子的数量；③这些可用的位点连接成一个连续的扩散通路。这些特点使得格子气KMC可以很好地适用于固态电解质材料的模拟。Morgan等[20]根据石榴石构型固态电解质的结构特征[图11(a)～(c)]，采用晶格气体MC模拟预测了石榴石结构固态电解质Li x La3Zr2O12的相关系数变化趋势，发现当xLi = 3(来自位置能量)和xLi = 6(来自最近邻斥力)时，其存在特别强的相关效应[图11(d)～(e)]。

图11

新窗口打开| 下载原图ZIP| 生成PPT

图11   石榴石框架中离子迁移扩散网络示意图与相关系数随迁移离子浓度的变化：(a) 石榴石结构中由ZrO6 八面体和LaO8 十二面体连接的离子输运通道的三维图。迁移离子随机分布在四面体和八面体间隙中。相互连通的四面体隙和八面体间隙形成三维离子输运网络；(b) 石榴石结构中环形结构示意图，其中四面体和八面体空位被部分锂离子占据。应该注意的是，并不是所有的空位都需要被锂离子占据；(c) 石榴石离子输运通道的二维几何连接。一个八面体间隙与两个四面体间隙相连，一个四面体间隙与四个八面体间隙相连(四个八面体空位与四面体的每个面相连，图中仅展示两个八面体做示例)。在(a)、(b)和(c)中用蓝色表示八面体空位，用黄色表示四面体空位；(d) 集体相关系数随着Li离子浓度(xLi)的变化[20]；(e) 自相关系数随着Li离子浓度(xLi)的变化[20]， xLi 指的是Li x La3Zr2O12 中Li离子的化学计量数，取值0～9

Fig. 11   Schematic diagram of ion diffusion network in garnet frame and the change of correlation coefficient with Li ion concentration: (a) three-dimensional drawings of ion transport channels connected by ZrO6 octahedron and LaO8 dodecahedron in garnet structures. Migrated ions are randomly distributed in tetrahedral and octahedral spaces. The interconnected tetrahedral gap and octahedral gap form a three-dimensional ion transport network; (b) schematic diagram of partial ring structure, in which tetrahedral and octahedral vacancies are partially occupied by lithium ions. It should be noted that not all vacancies need to be occupied by lithium ions; (c) two-dimensional geometric connections of garnet ion transport channels. One octahedral gap is connected to two tetrahedral gaps, and one tetrahedral gap is connected to four octahedral gaps (four octahedral gaps are connected to each face of the tetrahedron, two of which are not drawn in the two-dimensional topology). In (a), (b) and (c), octahedral vacancies are represented by blue and tetrahedral vacancies by yellow; (d) the change of collective correlation coefficient with Li ion concentration (xLi); (e) change of autocorrelation coefficient with Li ion concentration (xLi), xLi refers to the stoichiometric number of Li ions in Li x La3Zr2O12, ranging from 0 to 9


考虑到键价和理论等静态计算方法的局限性，Chen等[8]提出将键价和理论计算的迁移能垒数据与KMC模拟相结合的方法，得到了指定温度下的绝对电导率[图12(b)]，这显然优于仅考虑几何效应计算得到的电导率变化趋势[图12(a)]。由于迁移离子的具体排布对离子迁移能垒具有极大的影响，所以上述的近似也不够精确。为此，Ven等[46]基于第一性原理团簇展开的方式，对KMC模拟的势垒参数进行了计算。他们发现在无序材料中，应用团簇展开可以近似获得任何构型下的迁移势垒。这意味着在KMC模拟中，任何时刻都能近似地保持细致平衡，并能在热力学平衡状态下计算出所需的动力学量。作者进一步给出了这一研究范式在Li x CoO2中的应用。结果表明，迁移机制和激活能垒在很大程度上取决于迁移锂离子周围近邻锂离子与锂空位的排布[图12(c)]。通过模拟，作者发现在所有锂浓度下，锂在层状Li x CoO2中的扩散依赖于双空位机制。此外，由于激活能垒受离子浓度影响巨大，扩散系数随锂离子浓度x的变化可以达到几个数量级。综上所述，KMC的模拟一般依赖于其他计算方法的输入(包括模型结构，迁移能垒等)，这对计算结果精度的影响是巨大的。

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12   几个KMC计算过程的对比：(a) LLZO的通道连接关系与仅考虑几何效应下的离子电导率示意图[20]；(b) BVSE计算得到的离子输运通道与通过KMC模拟得到的电导率与实验值的比较。实心蓝色曲线和填充圆表示KMC所模拟的电导率。实红线表示有效计算的电导率(来自高温KMC模拟)，虚线表示对较低温度下离子电导率的推断。三角形绿色曲线表示实验电导率数据[8]；(c) 几种典型的离子排布状态与基于局域结构团簇展开的KMC计算得到的离子扩散系数[46]

Fig. 12   Comparison of several KMC calculation processes: (a) schematic diagram of channel connection and ionic conductivity of LLZO only considering geometric effects; (b) ionic transport channel calculated by BVSE and the comparison of ion conductivity calculated by KMC with experimental observation. Solid blue curves and filled circles represent the ion conductivity simulated by KMC. The solid red line represents the ion conductivity effectively calculated (from the KMC simulation at high temperatures) and the dashed line represents the inferred ionic conductivity at lower temperatures. The triangular green curve represents the experimental observation conductivity; (c) several typical ion distribution states and the ion diffusion coefficient calculated by KMC based on cluster expansion[46]


考虑到KMC模拟可以有效地预测和分析晶体材料离子电导率，Grieshammer等[57]开发了一款名为MOCASSIN的固态离子导体MC计算软件，利用该软件对富镧橄榄石中的氧间隙迁移以及掺杂全水合锆酸钡中的质子传导进行了模拟研究。模拟结果揭示了缺陷相互作用对离子电导率的影响。这些效应的组合可能会导致固态离子材料中意想不到的离子输运行为，特别是在具有较多可移动离子的模拟体系中。除此之外，Hoffmann等[58]也开发了一款名为Kmos的软件，在软件描述的系统中，所有的反应空间都可以抽象为一系列离散的空间格点(即为前文的格子气模型)。Kmos通过代码生成新的可用事件来不断演化，以实现较高的模拟效率。

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

除上述对离子导体材料的模拟之外，MC模拟在晶体生长以及SEI膜层生长等领域也有着重要的应用[69]。基于晶体生长理论的计算模拟是深入研究晶体生长机理，进一步指导晶体生长实验的一项重要技术。MC在模拟晶体生长的过程中分别考虑了原子的沉积、表面扩散、吸附等过程。在KMC模拟中，晶体的生长过程可被视为一个随机过程，使用某些特定的生长模型来描述其中随机事件的发生概率。通常具体模拟步骤如下：①建立微观模型；②对模型进行描述；③模型发生演变；④比较演变前后的结果确定演变概率，这与第二部分的相关总结类似。图13(a)～(g)展示了简立方晶格模型的表面在不同配位下的晶体生长过程，可以看出这几种配位对晶体长大的速度与方向均会产生极大的影响[70]。Barai等[71]基于简化的二维格子气模型通过KMC模拟建立了一个多尺度计算模型，用以捕捉晶体初生离子的形核和生长以及它们聚集成二次过渡金属氢氧化物前驱体离子。结果显示随着溶液pH的增加和氨浓度的降低，二次颗粒的尺寸减小。结合模拟给出的相图[图13(h)]，相关研究可以为实验合成提供有效指导。Chen等[72]采用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等[73]利用KMC研究了充电期间石墨电极上SEI的生长过程。结果显示，锂离子的嵌入速率和充电时间均受锂还原速率和通过SEI的锂传输速率的控制。在较低的温度下，Li的还原速率较高，但SEI的电阻也较高，所以会限制扩散过程。另一方面，在较高温度下，较低的Li还原速率会导致总充电时间增加[图13(l)]。

图13

新窗口打开| 下载原图ZIP| 生成PPT

图13   (a)～(g)为简单立方晶格表面不同配位下的晶体生长过程[70]；(h) 晶体尺寸与溶液pH值与氨含量依赖关系的相图。图中显示了在粒径分布为高斯分布的情况下，溶液的平均粒径及其对应的一阶标准差。这里显示的所有数据都是从计算模型中获取的。绿色区域表示次级颗粒较小，黄色区域表示次级颗粒较大。浅蓝点的二次粒径分布标准差较小，紫色点的二次粒径分布标准差较大。在pH值和氨含量方面，沉积较大(Dpart约8 μm)和较小(Dpart约4 μm)的二次活性颗粒的最佳操作条件也在图中突出显示(用红色圆圈表示)[71]；(i) T = 1100 K下的不同Si/C比条件下，两种晶体在沉积速率 F = 0.1 ML/s下，6H-SiC在生长中的比例[72]；(j) T = 1400 K下的不同Si/C比条件下，两种晶体在沉积速率 F = 0.1 ML/s下，6H-SiC在生长中的比例[72]；(k) SEI膜的生成过程，一般分为吸附、解吸、表面扩散和钝化等步骤；(l)SEI膜厚度、充电速度与温度的依赖关系[73]

Fig.13   (a)-(g) the crystal growth process on the surface of a simple cubic lattice under different coordination; (h) Phase diagram of the dependence of crystal size on solution pH and ammonia content. The figure shows the mean particle size of the solution and its corresponding first standard deviation in the case of gaussian particle size distribution. All the data is taken from the computational model. Green areas indicate smaller secondary particles and yellow areas indicate larger secondary particles. The standard deviation of the secondary particle size distribution of the light blue dot is smaller, and the purple dot is larger. Optimum operating conditions, in terms of pH and ammonia content, for precipitating relatively larger (Dpart ~8 μm) and smaller (Dpart ~4 μm) sized secondary active particles, have also been highlighted within the figure (by red circles); (i) under different Si/C ratio conditions at 1100 K, the ratio of 6H-SiC in the growth of the two crystals at the deposition rate F =0.1 ML/s; (j) under different Si/C ratio conditions at 1400 K, the ratio of 6H-SiC in the growth of the two crystals at the deposition rate F =0.1 ML/s; (k) the generation process of SEI film is generally divided into adsorption, desorption, surface diffusion and passivation, etc. (l) the dependence of SEI film thickness, charging speed and temperature


目前对晶体生长的MC模拟仍存在许多问题，比如Maazi等[53]模拟晶体生长时，在随机矩阵位点上需要重复数百万次晶粒生长过程。这对于较大的晶格系统来说计算时间会很长。为了提高晶粒生长模拟的计算速度，研究人员对经典MC模型进行了两种修正：①所有阵点都以相同的概率进行随机重定向，即每个蒙特卡罗步都不重复；②位置能量除以一个参数，该参数表示离子存在时，作用在晶界上的有效作用。这些举措极大地加速了大颗粒的生长，减弱了邻近小颗粒的生长，有效地提高MC模拟的速度。

4 结论与展望

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

参考文献

GROSJEAN C, MIRANDA P H, PERRIN M, et al. Assessment of world lithium resources and consequences of their geographic distribution on the expected development of the electric vehicle industry[J]. Renewable and Sustainable Energy Reviews, 2012, 16(3): 1735-1744.

[本文引用: 1]

卢侠, 李泓. 锂电池基础科学问题(Ⅱ)——电池材料缺陷化学[J]. 储能科学与技术, 2013, 2(2): 157-164.

LU X, LI H. Fundamental scientific aspects of lithium batteries (Ⅱ)—Defect chemistry in battery materials[J]. Energy Storage Science and Technology, 2013, 2(2): 157-164.

FLURI A, MARCOLONGO A, RODDATIS V, et al. Enhanced proton conductivity in Y-doped BaZrO3 via strain engineering[J]. Advanced Science, 2017, 4(12): doi: 10.1002/advs.201700467.

[本文引用: 1]

YAN J, WANG Q, WEI T, et al. Supercapacitors: Recent advances in design and fabrication of electrochemical supercapacitors with high energy densities[J]. Advanced Energy Materials, 2014, 4(4): doi: 10.1002/aenm.201300816.

[本文引用: 1]

MANTHIRAM A, YU X, WANG S. Lithium battery chemistries enabled by solid-state electrolytes[J]. Nature Reviews Materials, 2017, 2: doi: 10.1038/natrevmats.2016.103.

[本文引用: 1]

HE B, YE A, CHI S, et al. CAVD, towards better characterization of void space for ionic transport analysis[J]. Scientific Data, 2020, 7: doi: 10.1038/s41597-020-0491-x.

[本文引用: 2]

BLATOV V A, SHEVCHENKO A P, PROSERPIO D M. Applied topological analysis of crystal structures with the program package ToposPro[J]. Crystal Growth & Design, 2014, 14(7): 3576-3586.

[本文引用: 1]

CHEN H M, WONG L L, ADAMS S. SoftBV—A software tool for screening the materials genome of inorganic fast ion conductors[J]. Acta Crystallographica Section B, Structural Science, Crystal Engineering and Materials, 2019, 75(Pt 1): 18-33.

[本文引用: 4]

WONG L L, PHUAH K C, DAI R Y, et al. Bond valence pathway analyzer—An automatic rapid screening tool for fast ion conductors within softBV[J]. Chemistry of Materials, 2021, 33(2): 625-641.

[本文引用: 1]

ZOU Z Y, MA N, WANG A P, et al. Relationships between Na+ distribution, concerted migration, and diffusion properties in rhombohedral NASICON[J]. Advanced Energy Materials, 2020, 10(30): doi: 10.1002/aenm.202001486.

[本文引用: 2]

MOON J, LEE B, CHO M, et al. Ab initio and kinetic Monte Carlo simulation study of lithiation in crystalline and amorphous silicon[J]. Journal of Power Sources, 2014, 272: 1010-1017.

[本文引用: 4]

MURCH G E. The haven ratio in fast ionic conductors[J]. Solid State Ionics, 1982, 7(3): 177-198.

[本文引用: 1]

LU X Y. Application of the Nernst-Einstein equation to concrete[J]. Cement and Concrete Research, 1997, 27(2): 293-302.

[本文引用: 1]

DE KLERK N J J, VAN DER MAAS E, WAGEMAKER M. Analysis of diffusion in solid-state electrolytes through MD simulations, improvement of the Li-ion conductivity in β-Li3PS4 as an example[J]. ACS Applied Energy Materials, 2018, 1(7): 3230-3242.

[本文引用: 1]

HE X, ZHU Y, MO Y. Origin of fast ion diffusion in super-ionic conductors[J]. Nature Communications, 2017, 8: doi: 10.1038/ncomms15893.

[本文引用: 3]

ZHU P P, SMITH R W. Dynamic simulation of crystal growth by Monte Carlo method—I. Model description and kinetics[J]. Acta Metallurgica et Materialia, 1992, 40(4): 683-692.

[本文引用: 1]

BÉJAUD R, DURINCK J, BROCHARD S. Twin-interface interactions in nanostructured Cu/Ag: Molecular dynamics study[J]. Acta Materialia, 2018, 144: 314-324.

[本文引用: 1]

LYU D, WANG W, LIU J P, et al. Phase diagrams and magnetic properties of a ferrimagnetic Ising bilayer superlattice: A Monte Carlo study[J]. Journal of Magnetism and Magnetic Materials, 2018, 465: 348-359.

[本文引用: 4]

KAR P, HARINIPRIYA S. Modeling of lithium ion batteries employing grand canonical Monte Carlo and multiscale simulation[J]. Journal of the Electrochemical Society, 2014, 161(5): A726-A735.

[本文引用: 2]

MORGAN B J. Lattice-geometry effects in garnet solid electrolytes: A lattice-gas Monte Carlo simulation study[J]. Royal Society Open Science, 2017, 4(11): doi: 10.1098/rsos.170824.

[本文引用: 9]

WANG W, BI J L, LIU R J, et al. Effects of the single-ion anisotropy on magnetic and thermodynamic properties of a ferrimagnetic mixed-spin (1, 3/2) cylindrical Ising nanowire[J]. Superlattices and Microstructures, 2016, 98: 433-447.

[本文引用: 1]

MALIK R, ZHOU F, CEDER G. Phase diagram and electrochemical properties of mixed olivines from first-principles calculations[J]. Physical Review B, 2009, 79(21): doi: 10.1103/PhysRevB.79.214201.

[本文引用: 1]

RAN Y, ZOU Z, LIU B, et al. Towards prediction of ordered phases in rechargeable battery chemistry via group-subgroup transformation[J]. npj Computational Materials, 2021, 7: 1-11.

[本文引用: 1]

OZAWA K. Lithium-ion rechargeable batteries with LiCoO2 and carbon electrodes: The LiCoO2/C system[J]. Solid State Ionics, 1994, 69(3/4): 212-221.

[本文引用: 1]

KIM J H, MYUNG S T, YOON C S, et al. Comparative study of LiNi0.5Mn1.5O4- δ and LiNi0.5Mn1.5O4 cathodes having two crystallographic structures: Fd3̄m and P4332[J]. Chemistry of Materials, 2004, 16(5): 906-914.

[本文引用: 1]

WU H, PI J C, LIU Q, et al. All-inorganic lead free double perovskite Li-battery anode material hosting high Li+ ion concentrations[J]. The Journal of Physical Chemistry Letters, 2021, 12(17): 4125-4129.

[本文引用: 1]

HUANG X. Advances in Li-ion batteries and related materials[J]. Advances in Materials in China, 2010, 29(8): 46-52.

[本文引用: 1]

GEYER C J. Practical Markov chain Monte Carlo[J]. Statistical Science, 1992, 7(4): 473-483.

[本文引用: 5]

DARLING R, NEWMAN J. Dynamic Monte Carlo simulations of diffusion in LiyMn2O4[J]. Journal of the Electrochemical Society, 1999, 146(10): 3765-3772.

SHIMIZU A, TACHIKAWA H. Dynamics behavior of lithium in graphite lattice: MD calculation approach[J]. Journal of Physics and Chemistry of Solids, 2000, 61(12): 1895-1899.

DUANE S, KENNEDY A D, PENDLETON B J, et al. Hybrid Monte Carlo[J]. Physics Letters B, 1987, 195(2): 216-222.

[本文引用: 1]

HARDY J, DE PAZZIS O, POMEAU Y. Molecular dynamics of a classical lattice gas: Transport properties and time correlation functions[J]. Physical Review A, 1976, 13(5): 1949-1961.

[本文引用: 4]

CARLO C M. Markov chain Monte Carlo and Gibbs sampling[R]. Lecture notes for EEB, 2004, 581.

[本文引用: 2]

METROPOLIS N, ULAM S. The Monte Carlo method[J]. Journal of the American Statistical Association, 1949, 44(247): 335-341.

[本文引用: 2]

BROWNE C B, POWLEY E, WHITEHOUSE D, et al. A survey of Monte Carlo tree search methods[J]. IEEE Transactions on Computational Intelligence and AI in Games, 2012, 4(1): 1-43.

[本文引用: 1]

ZSIGMOND G, LIEUTENANT K, MEZEI F. Monte Carlo simulations of neutron scattering instruments by VITESS: Virtual instrumentation tool for ESS[J]. Neutron News, 2002, 13(4): 11-14.

FRANCIS Z, INCERTI S, ZEIN S A, et al. Monte Carlo simulation of SARS-CoV-2 radiation-induced inactivation for vaccine development[J]. Radiation Research, 2021, 195(3): 221-229.

TURCHIN V F. On the computation of multidimensional integrals by the Monte-Carlo method[J]. Theory of Probability & Its Applications, 1971, 16(4): 720-724.

YANG Z P, MING J, QIU C X, et al. A multigrid multilevel Monte Carlo method for stokes-darcy model with random hydraulic conductivity and beavers-Joseph condition[J]. Journal of Scientific Computing, 2022, 90(2): 1-30.

[本文引用: 1]

BENES N E, BOUWMEESTER H J M, VERWEIJ H. Multi-component lattice gas diffusion[J]. Chemical Engineering Science, 2002, 57(14): 2673-2678.

[本文引用: 1]

BROOKS S, GELMAN A, JONES G, et al. Handbook of Markov chain Monte Carlo[M]. Florida: CRC Press, 2011.

[本文引用: 1]

CHENG J R, YUAN X H, ZHAO L, et al. GCMC simulation of hydrogen physisorption on carbon nanotubes and nanotube arrays[J]. Carbon, 2004, 42(10): 2019-2024.

[本文引用: 1]

COSOLI P, FERRONE M, PRICL S, et al. Hydrogen sulphide removal from biogas by zeolite adsorption: Part I. GCMC molecular simulations[J]. Chemical Engineering Journal, 2008, 145(1): 86-92.

[本文引用: 1]

DO D D, DO H D. Pore characterization of carbonaceous materials by DFT and GCMC simulations: A review[J]. Adsorption Science & Technology, 2003, 21(5): 389-423.

[本文引用: 1]

GAVILÁN-ARRIAZU E M, PINTO O A, LÓPEZ M B, et al. The kinetic origin of the Daumas-Hérold model for the Li-ion/graphite intercalation system[J]. Electrochemistry Communications, 2018, 93: 133-137.

[本文引用: 5]

DER V A, CEDER G, ASTA M, et al. First-principles theory of ionic diffusion with nondilute carriers[J]. Physical Review B, 2001, 64(18): doi: 10.1103/PhysRevB.64.184307.

[本文引用: 7]

METHEKAR R N, NORTHROP P W, CHEN K J, et al. Kinetic Monte Carlo simulation of surface heterogeneity in graphite anodes for lithium-ion batteries: Passive layer formation[C]//Proceedings of the 2011 American Control Conference, 2011: 1512-1517.

[本文引用: 1]

AKKERMANS R L, SPENLEY N A, ROBERTSON S H. Monte Carlo methods in materials studio[J]. Molecular Simulation, 2013, 39(14/15): 1153-1164.

[本文引用: 1]

CHATTERJEE A, VLACHOS D G. An overview of spatial microscopic and accelerated kinetic Monte Carlo methods[J]. Journal of Computer-Aided Materials Design, 2007, 14(2): 253-308.

[本文引用: 1]

FRISCH U, HASSLACHER B, POMEAU Y. Lattice-gas automata for the navier-stokes equation[J]. Physical Review Letters, 1986, 56(14): 1505-1508.

[本文引用: 1]

LEE J, URBAN A, LI X, et al. Unlocking the potential of cation-disordered oxides for rechargeable lithium batteries[J]. Science, 2014, 343(6170): 519-522.

[本文引用: 1]

WANG W, CHEN D D, LYU D, et al. Monte Carlo study of magnetic and thermodynamic properties of a ferrimagnetic Ising nanoparticle with hexagonal core-shell structure[J]. Journal of Physics and Chemistry of Solids, 2017, 108: 39-51.

[本文引用: 1]

MAAZI N, BOULECHFAR R. A modified grain growth Monte Carlo algorithm for increased calculation speed in the presence of Zener drag effect[J]. Materials Science and Engineering: B, 2019, 242: 52-62.

[本文引用: 1]

LONGONE P, MARTÍN Á, RAMIREZ-PASTOR A J. Lattice-gas Monte Carlo study of sI clathrate hydrates of ethylene: Stability analysis and cell distortion[J]. Fluid Phase Equilibria, 2020, 521: doi: 10.1016/j.fluid.2020.112739.

[本文引用: 1]

METROPOLIS N, ROSENBLUTH A W, ROSENBLUTH M N, et al. Equation of state calculations by fast computing machines[J]. The Journal of Chemical Physics, 1953, 21(6): 1087-1092.

[本文引用: 2]

YANG Y, WANG W, LYU D, et al. Monte Carlo study of magnetic behaviors in a quadrangle ferrimagnetic Ising nanoisland[J]. Journal of Physics and Chemistry of Solids, 2018, 120: 109-122.

[本文引用: 2]

GRIESHAMMER S, EISELE S. Kinetic Monte Carlo simulations for solid state ionics: Case studies with the MOCASSIN program[J]. Diffusion Foundations, 2021, 29: 117-142.

[本文引用: 2]

HOFFMANN M J, MATERA S, REUTER K. Kmos: A lattice kinetic Monte Carlo framework[J]. Computer Physics Communications, 2014, 185(7): 2138-2150.

[本文引用: 1]

XIAO P H, HENKELMAN G. Kinetic Monte Carlo study of Li intercalation in LiFePO4[J]. ACS Nano, 2018, 12(1): 844-851.

[本文引用: 4]

HE B, MI P H, YE A J, et al. A highly efficient and informative method to identify ion transport networks in fast ion conductors[J]. Acta Materialia, 2021, 203: doi: 10.1016/j.actamat.2020.116490.

[本文引用: 1]

THANGADURAI V, NARAYANAN S, PINZARU D. Garnet-type solid-state fast Li ion conductors for Li batteries: Critical review[J]. Chemical Society Reviews, 2014, 43(13): 4714-4727.

[本文引用: 2]

ETACHERI V, MAROM R, ELAZARI R, et al. Challenges in the development of advanced Li-ion batteries: A review[J]. Energy & Environmental Science, 2011, 4(9): doi: 10.1039/C1EE01598B.

[本文引用: 1]

YAMADA A, CHUNG S C, HINOKUMA K. Optimized LiFePO4 for lithium battery cathodes[J]. Journal of the Electrochemical Society, 2001, 148(3): doi: 10.1149/1.1348257.

[本文引用: 1]

HIN C. Kinetic Monte Carlo simulations of anisotropic lithium intercalation into LixFePO4 electrode nanocrystals[J]. Advanced Functional Materials, 2011, 21(13): 2477-2487.

[本文引用: 2]

ZHENG T, DAHN J R. Lattice-gas model to understand voltage profiles of LiNixMn2– xO4/Lielectrochemical cells[J]. Physical Review B, 1997, 56(7): 3800-3805.

[本文引用: 2]

PERSSON K, HINUMA Y, MENG Y S, et al. Thermodynamic and kinetic properties of the Li-graphite system from first-principles calculations[J]. Physical Review B, 2010, 82(12): doi: 10.1103/PhysRevB.82.125416.

[本文引用: 2]

OUYANG C Y, SHI S Q, WANG Z X, et al. The effect of Cr doping on Li ion diffusion in LiFePO4 from first principles investigations and Monte Carlo simulations[J]. Journal of Physics: Condensed Matterial, 2004, 16(13): 2265-2272.

[本文引用: 2]

KUMAR P P, YASHONATH S. Ionic conduction in the solid state[J]. Journal of Chemical Sciences, 2006, 118(1): 135-154.

[本文引用: 1]

HUANG C M, JOANNE C L, PATNAIK B S V, et al. Monte Carlo simulation of grain growth in polycrystalline materials[J]. Applied Surface Science, 2006, 252(11): 3997-4002.

[本文引用: 1]

KOTRLA M. Numerical simulations in the theory of crystal growth[J]. Computer Physics Communications, 1996, 97(1/2): 82-100.

[本文引用: 2]

BARAI P, FENG Z G, KONDO H, et al. Multiscale computational model for particle size evolution during coprecipitation of Li-ion battery cathode precursors[J]. The Journal of Physical Chemistry B, 2019, 123(15): 3291-3303.

[本文引用: 2]

CHEN X J, ZHAO H, AI W S. Study on the competitive growth mechanism of SiC polytypes using Kinetic Monte Carlo method[J]. Journal of Crystal Growth, 2021, 559: doi: 10.1016/j.jcrysgro. 2021.126042.

[本文引用: 3]

HAO F, LIU Z X, BALBUENA P B, et al. Mesoscale elucidation of solid electrolyte interphase layer formation in Li-ion battery anode[J]. The Journal of Physical Chemistry C, 2017, 121(47): 26233-26240.

[本文引用: 2]

HUANG J, ZHONG P, HA Y, et al. Non-topotactic reactions enable high rate capability in Li-rich cathode materials[J]. Nature Energy, 2021, 6(7): 706-714.

[本文引用: 1]

SHI S Q, GAO J, LIU Y, et al. Multi-scale computation methods: Their applications in lithium-ion battery research and development[J]. Chinese Physics B, 2016, 25(1): doi: 10.1088/1674-1056/25/1/018212.

[本文引用: 1]

/

〈 〉