劉金平，蒲博偉，鄒喆乂，李銘清，丁昱清，任 元，羅亞橋，李 杰，李亞捷，王 達，何 冰，施思齊,5,6

（1上海大學材料科學與工程學院，上海 200444；2湘潭大學材料科學與工程學院，湖南 湘潭 411105；3加州大學歐文分校物理與天文系，美國 加州 92697；4上海大學計算機工程與科學學院；5上海大學材料基因組工程研究院，上海 200444；6之江實驗室，浙江 杭州 311100）

明確離子導體的離子輸運與相變特性是制備新型全固態(tài)電池、電化學轉換器件等離子型器件的關鍵[1-3]。離子導體的典型結構包括在輸運網(wǎng)絡中嵌入或遷移實現(xiàn)電荷儲存與傳輸?shù)倪w移離子以及構建離子輸運網(wǎng)絡的僅原點附近振動的骨架離子[4-5]。高離子電導率的超離子導體的篩選過程通常采用高通量計算方法[如幾何分析方法(crystal structure analysis by voronoi decomposition，CAVD)[6-7]、配位成鍵理論(bond valence site energy，BVSE)[8-9]等]分析骨架離子中的離子輸運網(wǎng)絡。同時，借助基于牛頓力學的分子動力學(molecular dynamics，MD)模擬[10]或基于離子遷移實驗的動力學蒙特卡羅(kinetic Monte Carlo，KMC)模擬[11]分析遷移離子的輸運特性(例如離子電導率[12-13]、擴散系數(shù)[8,10,14]等)。借助群論、巨正則蒙特卡羅(grand canonical Monte Carlo，GCMC)模擬等熱力學模擬手段研究離子導體的結構演化。分子動力學(MD)模擬作為典型的動力學模擬手段因其計算精度高、模擬分子運動過程清晰等優(yōu)勢，在離子擴散模式分析[15]及晶體生長模擬[16-17]等方面具有廣泛應用。然而其內(nèi)稟缺陷限制了MD適用于更廣泛的應用場景。例如，第一性原理分子動力學模擬計算速度較慢、難以應用到較大體系的快速模擬；而平均場分子動力學模擬不易獲取其所需的高精度勢函數(shù)。因此，具有類似計算功能的KMC 模擬因計算速度快、擴展性強而被廣泛應用于評估離子導體的離子輸運特性。表1對幾種典型的MD模擬與MC模擬[18-20]進行了對比，盡管兩者的應用范圍存在眾多交疊，但基于計算原理的差異性，MC 模擬具有更廣泛的時空尺度。此外，在熱力學模擬領域(如相圖計算[21]、電壓平臺計算等)，相較于第一性原理相圖計算[22]以及空間群分解[23]等模擬手段，MC 模擬同樣在計算速度與適用性等方面存在巨大優(yōu)勢。

表1 MC模擬與MD模擬的對比[18-20,28-32]Table 1 Comparison between MC simulation and MD simulation[18-20,28-32]

目前，MC 模擬還存在計算精度不高、計算過程依賴于較多的輸入?yún)?shù)、蒙特卡羅時間與真實時間之間存在差異、無法捕捉細小的離子振動等不足之處。為促進MC模擬的發(fā)展與應用，本文介紹了其發(fā)展歷史，提煉出MC模擬的基本模擬范式與具體計算步驟，以期為科研工作者開展離子導體的相關物化特性研究鋪平道路。本研究給出了一套具體的KMC 計算程序，并對石榴石型固態(tài)電解質(zhì)的離子輸運特性進行了分析，供讀者參考。為清晰地了解MC模擬的適用情形，本文就幾種典型的離子導體材料(鋰離子電池正極材料[24-25]、負極材料[26]以及固態(tài)電解質(zhì)材料[27]等)介紹了MC 模擬的計算結果。最后，本文對當前MC模擬所面臨的挑戰(zhàn)以及未來可能突破的方向進行了總結，期待該方法的進一步發(fā)展。

1 蒙特卡羅模擬方法基礎

1.1 蒙特卡羅模擬方法

MC 模擬作為一種基于隨機抽樣統(tǒng)計的模擬手段[33]，由美國數(shù)學家馮·諾伊曼等[28]于20世紀40年代在美國的“曼哈頓計劃”中提出。到20世紀50年代，進一步提出了MC模擬中可適用于大多數(shù)物理平衡系統(tǒng)的著名算法：Metropolis 算法[34]。20 世紀60 年代，KMC 模擬開始出現(xiàn)。與此同時，混合蒙特卡羅(hybrid Monte Carlo，HMC)采樣算法[31]等也得到了廣泛地研究。進入21 世紀以來，隨著人工智能的快速崛起，蒙特卡羅樹搜索算法等更為先進實用的計算方案也得到了快速發(fā)展[35-39](圖1)。根據(jù)研究對象的不同，MC 模擬可分為兩類：一類是針對熱力學平衡過程中非時間關聯(lián)性問題的模擬，如常見的GCMC 模擬；另一類是對時間關聯(lián)的動力學問題的研究，如KMC 模擬。這兩種MC 模擬的最終目標都是盡可能地尋找徘徊在某一低勢能面附近的結構組態(tài)。從原理上講，GCMC 與KMC之間存在應用場景的重合，即巨正則體系依舊可以進行KMC模擬，反之亦然。多數(shù)KMC模擬與熱力學MC模擬均可以利用格子氣[32,40]模型建模，故而上述多數(shù)模擬均屬于格子氣蒙特卡羅(lattice-gas Monte Carlo，LGMC)模擬。而格子氣模型的MC 模擬多屬于典型的馬爾科夫鏈(Markov chain)[28,33,41]過程，所以這種MC模擬也可以叫作馬爾科夫鏈蒙特卡羅(Markov chain Monte Carlo，MCMC)。這幾種模擬手段之間的關系可以通過圖3(e)所示的維恩圖表示?；谏鲜鰩追NMC 模擬的計算過程大多相似，本文將在格子氣MC模擬的介紹中對體系描述、體系演變及新體系接受概率等細節(jié)問題進行重點介紹。

圖1 蒙特卡羅模擬的發(fā)展Fig.1 The development of Monte Carlo simulation

1.2 蒙特卡羅模擬過程

GCMC[42-43]模擬作為一種適用于巨正則系綜限制下經(jīng)典的熱力學MC模擬，被廣泛應用于電化學儲能材料的合成、多孔吸附等過程的研究[44]。在電化學儲能材料計算領域，GCMC常用于模擬石墨等碳基嵌鋰材料的結構特征。在模擬過程中，隨機選擇幾種基本的離子運動模式(比如置換、插入、脫出)，并對相應離子運動模式發(fā)生的概率進行計算[45]，一般可以通過式(2)所示的計算模式求得。值得指出的是，式中的能量變化不僅可以使用體系的始末態(tài)能量表示，也可以使用過渡態(tài)能量表示[46]。對于忽略動力學特性的熱力學模擬，僅考慮簡單的始末態(tài)能量是常見的。作為GCMC 模擬適用情形的補充，KMC 模擬則可以很好地描述離子導體結構演變的動力學過程。KMC 模擬的研究最早可以追溯到20 世紀60 年代對輻射損傷退火問題的模擬。目前，KMC 在表面吸附、離子擴散、晶體生長等領域得到了迅速發(fā)展和應用[28]。由于模擬機理不同，KMC 相對于分子動力學(MD)模擬可以對更大系統(tǒng)的動態(tài)演化過程給出更長時間的模擬[34,47]。離子導體中遷移離子的成功遷移并不頻繁，其狀態(tài)轉變會經(jīng)歷長時間的“嘗試”，只有極少的“嘗試”會成功。KMC 模擬不考慮嘗試過程只關注可以演變成功的狀態(tài)，這種狀態(tài)到狀態(tài)的轉變不受原子振動的限制，不局限于對每個原子(離子)的計算，因此可以實現(xiàn)更長時間尺度的模擬[48]。同時，KMC的離散化思想更適合對微觀離散體系進行模擬分析，這也使得離散化模型的搭建在KMC 模擬中占據(jù)著十分重要的地位[49]。LGMC[32,50]模擬即為一種典型的離散化處理方式。

GCMC 模擬與KMC 模擬均可以借鑒格子氣模型而發(fā)展成為一種特殊的LGMC模擬來描述離子在材料體系中的擴散與占位情況。在格子氣模型的搭建過程中，時間與空間被分割，按照離子之間的相互作用以及遷移、翻轉規(guī)則，將模擬體系簡化成一系列固定的空間格點與可遷移或翻轉的離子。圖2展示了幾種典型格子氣模型的示意圖。該模型中，離子位于格點上，沿著邊運動，離子與離子之間存在相互作用(比如吸引或者排斥作用)。Ising 模型[18]、滲流模型[51]等均是對格子氣模型的某種處理與變形。

圖2 MC模擬中常見的幾種二維格子氣模型：(a)二維正方格子氣模型；(b)二維六邊格子氣模型；(c)二維隨機格子氣模型(紫色的球代表格點上的離子)Fig.2 Several two-dimensional lattice-gas model which used in MC simulation commonly:(a)twodimensional square lattice-gas model;(b)twodimensional hexagonal lattice-gas model;(c)random lattice-gas model(The purple spheres represent the ions at the lattice sites)

式中，m為總的位點數(shù)目；Esite_j為第j個位點被占據(jù)所引起的能量變化；Enn為最近鄰離子間的相互作用能；cij和cj為前置系數(shù)(其值取決于結構中離子的排布狀態(tài))。如果位點“j”被占據(jù)，則cj= 1，反之cj= 0。如果位點“j”和位點“i”同時被占據(jù)，并且它們互為最近鄰位點，則cij= 1，反之cij= 0。另外，由Metropolis 等[55]提出的系統(tǒng)演變概率計算公式具有精度較高、可有效避免結構在演變過程中進入局域能量極小值等優(yōu)點，所以在MC模擬中得到廣泛應用(利用此算法的MC模擬也可以叫做Metropolis Monte Carlo，MTMC)。依照此算法，組態(tài)a到組態(tài)b的演變概率(Pa→b)為

式中，k為玻爾茲曼常數(shù)；T為溫度，ΔEa→b為系統(tǒng)從組態(tài)a 演變到組態(tài)b 的能量變化。ΔEa→b計算公式如下

式中，Ea(Eb)為組態(tài)a(b)的能量，值得指出的是，ΔEa→b在動力學模擬中也會考慮到過渡態(tài)，如遷移能壘的引入等。僅依托于上述內(nèi)容來計算離子的擴散系數(shù)、離子的遷移率等時間相關量是困難的。為了進一步描述材料結構隨時間的一系列演變過程進而計算相鄰組態(tài)之間的時間關聯(lián)性，常常通過數(shù)學近似的方式獲得真實的時間步長[20,46]

式中，Kij表示“總遷移概率”；μ是隨機數(shù)[μ∈(0,1)]。由于總概率的計算十分耗時，而在MC模擬過程中一個蒙特卡羅步(MCS)內(nèi)的離子遷移可認為是同步發(fā)生的，所以采用MCS 取代真實時間是一種有效的簡化計算方案[56]。這種簡化已在傳統(tǒng)的熱力學計算中被廣泛采用。對于MC理論在材料學中的應用，其基本思想是將時-空尺度進行分離，因而其主要用于模擬有序系統(tǒng)，很難應用于無定形系統(tǒng)。對于MC 模擬中相關名詞的概念與定義在表2進行了總結。

表2 MC模擬中的常用概念Table 2 Common concepts in MC simulation

2 蒙特卡羅模擬的計算范式歸納與示例

MC 模擬已被應用到眾多領域，但是如何提高其模擬精度與效率一直是眾多科研工作者的重要研究方向[57]。為了清晰了解蒙特卡羅模擬的計算過程，此章節(jié)就其計算流程和步驟進行了總結與提煉。一般來說，模型的搭建應盡可能地復現(xiàn)材料的結構特征，比如通過幾何分析方式或BVSE來分析材料框架以及離子傳輸路徑，通過格子氣模型分析常見的晶態(tài)離子導體，通過Ising 模型來分析磁性材料等。此外，合適的模型描述也是決定相關模擬精度的關鍵。通常研究者會通過對模型能量(如內(nèi)能)的計算，進而以能量的變化來確定演變方向以及對應構型出現(xiàn)的概率[56-59]。然而，目前還沒有一個完美的計算方案可以準確地計算材料構型轉變過程中的能量變化。為此，眾多科研工作者作了許多努力。比如Arriazu等[45]定義了如圖3(a)所示的晶格氣體模型，考慮了層內(nèi)Li離子、層間Li離子以及Li離子與骨架結構的相互作用能，利用GCMC 模擬了Li離子在石墨中的插層過程。在模擬中，作者采用一種純靜態(tài)的能量計算模式，即只根據(jù)系統(tǒng)的總能來確定此構型出現(xiàn)的概率。因其不考慮過渡態(tài)能量，所以該模擬方法在動力學模擬中是不夠精確的。對于固定的遷移或反轉，激活能壘更能準確反映其能量的變化，因此Ven 等[46]提出了如圖3(b)所示的離子遷移能壘。為避免因離子往返遷移能壘不同而造成的計算困難，其引入了一個動力學解析的激活勢壘(ΔEKRA)作為能壘的近似，同時還考慮了目標離子周圍離子的排布，采用團簇展開的策略，近似計算了離子的遷移概率。然而該近似不僅忽略了離子遷移的方向性，也忽略了多離子遷移模式。遷移模式的選擇(材料結構演變途徑)是MC 模擬需要解決的第三個關鍵問題。He 等[15]證明了超離子導體中的離子快速擴散并不像典型固體中那樣通過孤立的離子遷移發(fā)生，而是通過具有低能量勢壘的多離子協(xié)同遷移進行[圖3(c)～(d)]。實際上，遷移模式的選擇是激活能計算的前提，而目前大多數(shù)報道都是基于單離子遷移模式，這會造成模擬結果的偏差。綜上所述，合理搭建構型，確定體系能量描述并選擇合適的構型演變模式，是決定模擬精度的關鍵。本團隊通過模塊化的思想提出了MC模擬離子導體材料離子輸運特性、相變特性的一般模式。這將大大方便后續(xù)開發(fā)者用來參考并快速修改，進而擴展MC模擬到具體材料體系。相關模擬的框架主要涵蓋以下幾個部分(圖4)：模擬目標的確定、模擬參數(shù)與模型的設定、結構演化方式的選擇、轉換過程中的能量計算、新結構接受的概率計算以及樣本數(shù)據(jù)的獲得。在具體問題中可能還需要增加某些新的模塊來處理特殊的需求?；诖肆鞒?，本文給出了KMC 計算石榴石結構離子導體離子輸運特性的一般范式，如圖5(a)所示，其中模型搭建依托于BVSE 與CAVD 的分析結果。模擬過程中的輸入輸出如表3所示。此外，基于對滲流模擬過程(模型搭建→局域連通性分析→滲流團簇分析→遷移網(wǎng)絡統(tǒng)計)的分析，提出了通過KMC模擬搭建滲流模擬初始構型，通過KMC給出的離子遷移模式判斷局域通堵，進而對離子輸運網(wǎng)絡分析的一般范式[圖5(b)]。KMC模擬程序框圖如圖6所示，相應的代碼與算例可從代碼庫https://gitee.com/jinpingliu/pKMC-code.git獲得。KMC模擬主程序計算步驟如下。

圖6 MC模擬在離子導體中的計算程序詳細流程框圖Fig.6 Flow diagram of MC calculation program in ionic conductor

表3 MC模擬常用的輸入與輸出Table 3 Some inputs and outputs of MC simulation

圖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

圖4 MC模擬的一般步驟，圖中列舉了MC樣本數(shù)據(jù)的生成過程，藍色一般為動力學模擬需要的模塊，橙色一般為熱力學模擬需要的模塊，綠色為兩者通用模塊，灰色是其他模擬中可能存在的模塊(此處粒子包含帶電離子與其他粒子，具體到離子導體中將用離子描述)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 KMC模擬流程與計算案例：(a)KMC對Li7La3Zr2O12(LLZO)的建模與模擬過程示意圖(模型構建方式或計算內(nèi)容依托于具體研究的問題，此處以石榴石結構Li7La3Zr2O12為例，采用BVSE與幾何分析方法搭建模型)；(b)KMC模擬與滲流模擬的結合示意圖(KMC模擬為滲流模擬提供結構模型、遷移離子排布、離子遷移模式等數(shù)據(jù)，實現(xiàn)對可達離子，有效載流子濃度等的計算)；(c)通過本KMC模擬程序計算得到LLZO的占據(jù)率變化趨勢與實驗結果[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

（1）讀取外部輸入信息，如位點連接關系及位點坐標、遷移能壘、反轉(遷移)模式等，根據(jù)設定的溫度、遷移離子濃度等物理參量搭建初始模擬系統(tǒng)。

（2）記錄此時系統(tǒng)的狀態(tài)，記為W0，計算此時系統(tǒng)的能量，為E0。

闊葉樹種，在京城及周邊木材市場上進入4季度銷路仍然暢通。與針葉原木市場相同的是，在京城以及周邊木材市場上經(jīng)營東北原木的商家普遍認同的仍是俄產(chǎn)木材。這一塊闊葉原木由于需求不減，資源品質(zhì)有保證，價格水平下行機會幾乎全無。另外，從俄方進口的北洋闊葉樹種原木像榆木、楸木、樺木、楊木、柞木、椴木和水曲柳，不僅需求仍然保持著前兩個月的強勁勢頭，其銷售價位也繼續(xù)堅挺上揚，例如北方市場最認可的水曲柳大徑級優(yōu)質(zhì)新材售價最強能夠沖高到5 000元/m3以上，一般材也就能賣到4 500元/m3左右。

（3）對系統(tǒng)進行隨機演變，得到新的組態(tài)，記為W1，計算此時系統(tǒng)的能量，記為E1。

（4）對比演變前后組態(tài)能量的變化，依據(jù)概率P接受新的組態(tài)，并更新時間。

（5）重復(3)～(4)過程，直到體系能量基本維持不變，開始記錄樣本數(shù)據(jù)。

（6）改變溫度，遷移離子濃度等物理參數(shù)，重復上述過程，記錄多組樣本數(shù)據(jù)，判斷是否結束KMC模擬。

（7）對樣本數(shù)據(jù)進行計算，得到遷移率、占據(jù)率等物理量。

需要補充的是，上述程序的開發(fā)過程依托于對固態(tài)電解質(zhì)材料的計算，其中相關的輸入?yún)?shù)使用到了本團隊前期開發(fā)的BVSE[60]與CAVD 的相關代碼[6]。這些輸入?yún)?shù)包括骨架結構的離子輸運通道，離子在相鄰位點之間遷移的平均勢壘等。圖5(c)～(d)給出了部分計算與實驗的結果對比，由圖可知，計算得到的LixLa3Zr2O12的離子占據(jù)情況與電導率變化均與實驗結果相符。

3 蒙特卡羅模擬在離子導體中的應用

多數(shù)MC 模擬均可對應到前文對MC 模擬的理論基礎與計算范式的介紹與梳理。為了描述MC模擬的應用場景，本節(jié)將從離子導體的相變過程、電壓平臺的變化過程、離子輸運特性以及晶體生長動力學特性的分析等方面進行詳細的介紹與總結。

3.1 相變分析與電壓平臺的蒙特卡羅模擬

在鋰離子電池中，正極材料的離子輸運特性直接決定了電池的容量與循環(huán)性能[62]。LiFePO4[63]作為一種被廣泛商用的正極材料，在平衡狀態(tài)下，會自發(fā)地分裂出貧鋰和富鋰兩相。如何快速獲取相的分離過程十分重要。Xiao 等[59]通過DFT 計算得到KMC 模擬的能壘輸入，搭建了如圖7(a)所示的模型。模擬結果顯示，LiFePO4和FePO4之間形成有序的Li0.5FePO4相[圖7(b)]。Hin等[64]通過KMC模擬了LixFePO4橄欖石納米晶體中鋰離子脫出與嵌入動力學特性。結果顯示，鋰離子嵌入過程中電池電壓的變化包含以下幾個階段：①鋰離子插入貧鋰相導致電池電勢下降；②富鋰相Li1-βFePO4成核后，電勢幾乎恒定；③貧鋰相完全演變?yōu)楦讳囅嘁院?，電池的電勢繼續(xù)下降[圖7(c)～(d)]。除了對LiFePO4不同Li 離子濃度與電壓分布關系的計算，Zheng等[65]還對LiNixMn2-xO4/Li 電池系統(tǒng)Ni 摻雜下的電壓分布進行了計算。作者首先構建了一個5×5×5的格子氣模型，包含1000 個可用位點。在模擬過程中，前500 個MC模擬步作為弛豫過程，接下來的500個模擬步用于熱力學量的計算，最終得到電壓分布與x值之間的依賴關系，如圖7(e)所示?？梢钥吹诫S著x的增大，4.7 V 電壓平臺的長度線性增大，4.1 V電壓平臺的長度線性減小，而總的電池容量則保持不變。當x=0.5時，只觀察到4.7 V的電壓平臺。有序到無序的轉變和4.1、4.7 V兩個電壓平臺的分裂都可以歸因于鋰離子與骨架結構之間的復雜相互作用。不同于上述對單個半電池的模擬，Kar 等[19]通過對兩個半電池(LiMn2O4正極和碳負極)進行建模，利用GCMC 模擬給出了溫度等參數(shù)變化對Li+嵌入與脫出開路電位、電池電流、電池電壓、自由能等電池參數(shù)的影響。在模擬過程中，當固定溫度T、體積V和化學勢μ時，每個位置的能量可由式(5)給出

式中，J是鋰離子之間的相互作用能，下標NN和NNN 分別表示最近鄰和次近鄰，如果位點“i”(“j”)被占據(jù)，ci(cj)取值為1，否則為0。JNN和JNNN分別為37.5 和25.0 meV。位點占據(jù)接受概率可通過式(6)進行描述

其中，e為電子電荷量。計算得到的開路電位結果如圖7(f)所示。這組模擬結果有效地解釋了測量開路電位剖面實驗中所觀察到的遲滯現(xiàn)象。不難發(fā)現(xiàn)，上述計算式(5)與(6)來自式(1)與(2)的變形，這種計算模式也被大量的MC模擬所采用。

圖7 MC模擬對正極材料相變與開路電壓的分析：(a)MC模擬模型，上圖為LixFePO4中Li的分布；下圖為對應的能量圖。綠色的球代表了鋰離子[59]；(b)在不考慮電荷間作用能的情況下，Li離子沿著a軸的排布。如圖所示，LiFePO4(LFP)和FePO4(FP)相之間形成了清晰的相邊界，Li0.5FePO4相收縮到一個有限的區(qū)域[59]；(c)LixFePO4橄欖石納米晶體在室溫下恒電流放電過程的MC模擬。灰點代表活性顆粒中的鋰原子，表面電流密度為0.5 A/m2。四組圖分別為0 s：初始固溶體，0.0001 s：形成兩個富鋰相并結合在一起，0.00056 s：富鋰相生長，10.81 s：貧鋰相幾乎完全被消耗；(d)電池電壓與活性材料中Li濃度(mol)的關系[64]；(e)MC模擬計算出的Li/LiNixMn2-xO4電池電壓分布，分別為x=0.1、x=0.2、x=0.3、x=0.5[65]；(f)全電池(LiMn2O4正極和碳負極)開路電位與放電期間正極材料中Ni占據(jù)率的關系[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 LixFePO4,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 LixFePO4 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/LiNixMn2-xO4 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計算的形成能揭示了晶態(tài)硅與非晶硅在鋰化過程中的相分離機制。晶態(tài)硅和非晶硅在鋰化作用下的體積膨脹和相變趨勢是相似的，而Li的擴散動力學在c-Si和a-Si之間則存在較大差異。在c-Si中，Li遷移勢壘為0.6 eV，隨著Li濃度的增加，遷移勢壘迅速減小到至0.4 eV[圖8(a)～(c)]。為了利用KMC模擬非晶硅中鋰的擴散，首先要利用體積函數(shù)推導出非晶硅中鋰遷移勢壘。通過KMC 模擬發(fā)現(xiàn)鋰在a-Si 中的擴散系數(shù)比在c-Si 中的擴散系數(shù)大一個數(shù)量級[圖8(d)]。這些研究有助于在原子尺度理解鋰化作用的機理，進一步闡明硅化鋰的相分離過程。此外，Persson 等[66]也通過DFT、團簇擴展和MC模擬預測了LixC6相圖，預測結果在高Li濃度(x＞0.5)與實驗結果相吻合[圖8(e)]。

圖8 Li離子排布對Li離子遷移勢壘的影響[(a)～(c)][11]；(d)a-Si與c-Si中Li離子的擴散系數(shù)[11]；(e)通過MC模擬獲得的LixC6的第一性原理相圖。圖中符號為：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)

為區(qū)分KMC 與GCMC 模擬在電池負極材料模擬中的差異性，Arriazu 等[45]分別通過KMC 與GCMC對鋰離子嵌入石墨的過程進行了分析。圖9(a)所示為Dumas-Hérold鋰/石墨插層化合物模型，其中藍色代表鋰離子，灰色代表石墨層，鋰離子的嵌入過程呈現(xiàn)出多階段(臺階)現(xiàn)象。離子的運動模式近似為二維運動[圖9(b)]，通過KMC模擬及GCMC模擬離子嵌入脫出隨時間的變化[圖9(c)]，他們發(fā)現(xiàn)Daumas-Hérold 模型[圖9(e)]中提出的結構是由于動力學限制造成的。如果系統(tǒng)的平衡是通過人為設置一個高交換電流密度(GCMC)的方式來實現(xiàn)，結構就變成了一個具有更低能量值[圖9(f)]的經(jīng)典Rüdorff-Hoffmann 模型[圖9(d)]。因此，可以通過改變?nèi)軇┑男再|(zhì)來增加實驗系統(tǒng)中的交換電流密度進而抑制DH 結構的產(chǎn)生。此外，使用超小的石墨薄片也有利于離子與電解質(zhì)的交換，進一步抑制DH結構的產(chǎn)生。

圖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 (a)Cr離子與Li離子的遷移能壘示意圖；(b)在每條一維通道中，鋰離子被Cr離子隔開的示意圖；(c)對于不同的模擬超胞尺寸，容量與摻雜量之間的依賴關系。數(shù)字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

對于固態(tài)電解質(zhì)，其離子輸運的動力學特性更加受到關注，KMC 在其中的應用更為廣泛。由于典型的固態(tài)電解質(zhì)存在三大特點[68]：①相鄰位點之間的離子遷移容易發(fā)生；②可供遷移離子占據(jù)的位點的數(shù)量大于遷移離子的數(shù)量；③這些可用的位點連接成一個連續(xù)的擴散通路。這些特點使得格子氣KMC 可以很好地適用于固態(tài)電解質(zhì)材料的模擬。Morgan等[20]根據(jù)石榴石構型固態(tài)電解質(zhì)的結構特征[圖11(a)～(c)]，采用晶格氣體MC 模擬預測了石榴石結構固態(tài)電解質(zhì)LixLa3Zr2O12的相關系數(shù)變化趨勢，發(fā)現(xiàn)當xLi=3(來自位置能量)和xLi=6(來自最近鄰斥力)時，其存在特別強的相關效應[圖11(d)～(e)]。

圖11 石榴石框架中離子遷移擴散網(wǎng)絡示意圖與相關系數(shù)隨遷移離子濃度的變化：(a)石榴石結構中由ZrO6八面體和LaO8十二面體連接的離子輸運通道的三維圖。遷移離子隨機分布在四面體和八面體間隙中。相互連通的四面體隙和八面體間隙形成三維離子輸運網(wǎng)絡；(b)石榴石結構中環(huán)形結構示意圖，其中四面體和八面體空位被部分鋰離子占據(jù)。應該注意的是，并不是所有的空位都需要被鋰離子占據(jù)；(c)石榴石離子輸運通道的二維幾何連接。一個八面體間隙與兩個四面體間隙相連，一個四面體間隙與四個八面體間隙相連(四個八面體空位與四面體的每個面相連，圖中僅展示兩個八面體做示例)。在(a)、(b)和(c)中用藍色表示八面體空位，用黃色表示四面體空位；(d)集體相關系數(shù)隨著Li離子濃度(xLi)的變化[20]；(e)自相關系數(shù)隨著Li離子濃度(xLi)的變化[20]，xLi指的是LixLa3Zr2O12中Li離子的化學計量數(shù)，取值0～9Fig.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 LixLa3Zr2O12,ranging from 0 to 9

考慮到鍵價和理論等靜態(tài)計算方法的局限性，Chen等[8]提出將鍵價和理論計算的遷移能壘數(shù)據(jù)與KMC 模擬相結合的方法，得到了指定溫度下的絕對電導率[圖12(b)]，這顯然優(yōu)于僅考慮幾何效應計算得到的電導率變化趨勢[圖12(a)]。由于遷移離子的具體排布對離子遷移能壘具有極大的影響，所以上述的近似也不夠精確。為此，Ven 等[46]基于第一性原理團簇展開的方式，對KMC 模擬的勢壘參數(shù)進行了計算。他們發(fā)現(xiàn)在無序材料中，應用團簇展開可以近似獲得任何構型下的遷移勢壘。這意味著在KMC 模擬中，任何時刻都能近似地保持細致平衡，并能在熱力學平衡狀態(tài)下計算出所需的動力學量。作者進一步給出了這一研究范式在LixCoO2中的應用。結果表明，遷移機制和激活能壘在很大程度上取決于遷移鋰離子周圍近鄰鋰離子與鋰空位的排布[圖12(c)]。通過模擬，作者發(fā)現(xiàn)在所有鋰濃度下，鋰在層狀LixCoO2中的擴散依賴于雙空位機制。此外，由于激活能壘受離子濃度影響巨大，擴散系數(shù)隨鋰離子濃度x的變化可以達到幾個數(shù)量級。綜上所述，KMC 的模擬一般依賴于其他計算方法的輸入(包括模型結構，遷移能壘等)，這對計算結果精度的影響是巨大的。

圖12 幾個KMC計算過程的對比：(a)LLZO的通道連接關系與僅考慮幾何效應下的離子電導率示意圖[20]；(b)BVSE計算得到的離子輸運通道與通過KMC模擬得到的電導率與實驗值的比較。實心藍色曲線和填充圓表示KMC所模擬的電導率。實紅線表示有效計算的電導率(來自高溫KMC模擬)，虛線表示對較低溫度下離子電導率的推斷。三角形綠色曲線表示實驗電導率數(shù)據(jù)[8]；(c)幾種典型的離子排布狀態(tài)與基于局域結構團簇展開的KMC計算得到的離子擴散系數(shù)[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]開發(fā)了一款名為MOCASSIN 的固態(tài)離子導體MC 計算軟件，利用該軟件對富鑭橄欖石中的氧間隙遷移以及摻雜全水合鋯酸鋇中的質(zhì)子傳導進行了模擬研究。模擬結果揭示了缺陷相互作用對離子電導率的影響。這些效應的組合可能會導致固態(tài)離子材料中意想不到的離子輸運行為，特別是在具有較多可移動離子的模擬體系中。除此之外，Hoffmann 等[58]也開發(fā)了一款名為Kmos的軟件，在軟件描述的系統(tǒng)中，所有的反應空間都可以抽象為一系列離散的空間格點(即為前文的格子氣模型)。Kmos 通過代碼生成新的可用事件來不斷演化，以實現(xiàn)較高的模擬效率。

3.3 晶體生長及電解質(zhì)界面膜(SEI 膜)生成的蒙特卡羅模擬

除上述對離子導體材料的模擬之外，MC 模擬在晶體生長以及SEI膜層生長等領域也有著重要的應用[69]?；诰w生長理論的計算模擬是深入研究晶體生長機理，進一步指導晶體生長實驗的一項重要技術。MC 在模擬晶體生長的過程中分別考慮了原子的沉積、表面擴散、吸附等過程。在KMC 模擬中，晶體的生長過程可被視為一個隨機過程，使用某些特定的生長模型來描述其中隨機事件的發(fā)生概率。通常具體模擬步驟如下：①建立微觀模型；②對模型進行描述；③模型發(fā)生演變；④比較演變前后的結果確定演變概率，這與第二部分的相關總結類似。圖13(a)～(g)展示了簡立方晶格模型的表面在不同配位下的晶體生長過程，可以看出這幾種配位對晶體長大的速度與方向均會產(chǎn)生極大的影響[70]。Barai 等[71]基于簡化的二維格子氣模型通過KMC 模擬建立了一個多尺度計算模型，用以捕捉晶體初生離子的形核和生長以及它們聚集成二次過渡金屬氫氧化物前驅(qū)體離子。結果顯示隨著溶液pH的增加和氨濃度的降低，二次顆粒的尺寸減小。結合模擬給出的相圖[圖13(h)]，相關研究可以為實驗合成提供有效指導。Chen等[72]采用KMC模擬研究了4H和6H型SiC的競爭生長過程。該模型以硅和碳原子為基本遷移離子，構建了格子氣KMC 模型。該模擬不僅考慮了同種原子之間的相互作用能(Si-Si)，還考慮了不同類型原子之間的相互作用能(Si-C)?；诖四Ｐ?，作者進一步研究了生長條件對4H-SiC 和6H-SiC 競爭生長的影響。結果表明，在生長過程中較高的溫度和高比例的Si 源會促進6H相的生長。同時，隨著溫度和Si/C比的升高，晶體表面原子聚集生長的現(xiàn)象越來越明顯[圖13(i)～(j)]。KMC 模擬SEI 的生成過程一般分為如下幾個過程：吸附、解吸、表面擴散以及鈍化(非活性)材料的形成[圖13(k)]?；诖肆鞒?，Hao 等[73]利用KMC研究了充電期間石墨電極上SEI的生長過程。結果顯示，鋰離子的嵌入速率和充電時間均受鋰還原速率和通過SEI的鋰傳輸速率的控制。在較低的溫度下，Li的還原速率較高，但SEI的電阻也較高，所以會限制擴散過程。另一方面，在較高溫度下，較低的Li還原速率會導致總充電時間增加[圖13(l)]。

圖13 (a)～(g)為簡單立方晶格表面不同配位下的晶體生長過程[70]；(h)晶體尺寸與溶液pH值與氨含量依賴關系的相圖。圖中顯示了在粒徑分布為高斯分布的情況下，溶液的平均粒徑及其對應的一階標準差。這里顯示的所有數(shù)據(jù)都是從計算模型中獲取的。綠色區(qū)域表示次級顆粒較小，黃色區(qū)域表示次級顆粒較大。淺藍點的二次粒徑分布標準差較小，紫色點的二次粒徑分布標準差較大。在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]模擬晶體生長時，在隨機矩陣位點上需要重復數(shù)百萬次晶粒生長過程。這對于較大的晶格系統(tǒng)來說計算時間會很長。為了提高晶粒生長模擬的計算速度，研究人員對經(jīng)典MC模型進行了兩種修正：①所有陣點都以相同的概率進行隨機重定向，即每個蒙特卡羅步都不重復；②位置能量除以一個參數(shù)，該參數(shù)表示離子存在時，作用在晶界上的有效作用。這些舉措極大地加速了大顆粒的生長，減弱了鄰近小顆粒的生長，有效地提高MC模擬的速度。

4 結論與展望

本文總結了MC模擬的發(fā)展歷史，歸納出了一套適用于多數(shù)MC 模擬的范式，包括以下幾個部分：模擬目標的確定、模擬參數(shù)與模擬模型的設定、結構演化方式的選擇、轉換過程中判據(jù)的提取(如體系能量)、新結構接受概率的計算、樣本數(shù)據(jù)的獲得等。根據(jù)此范式，列出了本團隊前期開發(fā)的一套典型MC模擬代碼供讀者參考使用，并給出了一個典型的計算案例：預測石榴石結構固態(tài)電解質(zhì)的有效載流子濃度、遷移離子分布與遷移離子濃度的依賴關系，模擬結果與實驗規(guī)律基本一致。除此之外，本文就其在正負極材料、電解質(zhì)以及相關界面中的典型熱力學與動力學計算案例進行了剖析，包括求解離子導體中的離子擴散問題、模擬遷移離子的分布特征與相關界面層的生長過程等。這些計算不僅涵蓋了實驗現(xiàn)象的微觀機理解釋，也包含了實驗結果的預測進而提供具體的實驗方案指導。在與實驗結合方面，MC 模擬可以基于某些實驗數(shù)據(jù)(如遷移離子的分布)搭建初始模型，提出針對具體研究體系的改進策略。如Huang等[74]基于實驗觀測，提出了層狀金屬氧化物中過渡金屬在八面體位點與四面體位點的遷移會改變遷移離子的輸運特性。鑒于每種模擬方法本身的局限性[75]，亦可以采用多方法聯(lián)用的方式，擴展MC模擬的輸出，比如文中提到的KMC 與滲流的結合以獲得有效載流子濃度。為了促進MC模擬在離子導體計算中的發(fā)展，本文總結了以下幾個MC方法的發(fā)展方向，包括：①在MC 模擬中如何快速精確地獲取所有可能發(fā)生的事件以及對應的轉變概率。在材料計算領域，這些事件及其描述主要是通過實驗測試以及第一性原理計算的方式獲得，比如通過密度泛函理論(DFT)計算遷移能壘，通過分子動力學預判離子的遷移模式等。然而這些參數(shù)的計算需要解決微觀演變中所涉及到的各種近似，比如基于第一性原理模擬的團簇展開就是對能量計算的一種簡化近似。除此之外，體系轉換速率因子也直接決定了模擬體系的演變速度，該參數(shù)一般需要近似為某一個常數(shù)，但對于一些復雜體系(比如存在多種遷移模式的晶態(tài)離子導體)，轉換因子(如離子的遷移嘗試頻率等)是極為復雜的，雖然通過聲子譜等手段可以進行較為精細的計算，但對于較大體系或者較長時間的模擬，相關的計算成本是巨大的。②晶格氣體MC模擬過程中僅考慮了當前對模擬體系有利的轉移路徑，這將導致模擬會很快進入某一亞穩(wěn)態(tài)。當前常用的Metropolis 算法雖然能夠一定程度上改善此問題，但對于存在多個亞穩(wěn)態(tài)的模擬體系，以上算法還是很難精確地尋找到系統(tǒng)的演化軌跡，并且其計算成本也極為高昂。③MC 模擬過程中，模擬事件之間存在先后順序，這不僅制約了其計算效率也影響了其模擬準確性。比如在極短時間內(nèi)離子的遷移存在一定的先后順序，如何確定遷移的先后順序會直接影響離子遷移的成功率。④MC 模擬始終在“自己的”時間域內(nèi)進行，與真實時間存在差異。雖然可以通過數(shù)學手段對MC 時間(也叫做蒙特卡羅步)與真實時間進行轉換，但這種轉換僅為統(tǒng)計上的一種近似，依賴于眾多初始輸入?yún)?shù)(如嘗試頻率等)。