1、.動力學蒙特卡洛方法(KMC)及相關討論動態模擬在目前的計算科學中占據著非常重要的位置。隨著計算能力和第一原理算法的發展，復雜的動態參數（擴散勢壘、缺陷相互作用能等）均可利用第一原理計算得出。因此，部分復雜的體系動態變化，如表面形貌演化或輻射損傷中缺陷集團的聚合-分解演變等，已可以較為精確的予以研究。KMC動力學蒙特卡洛方法（kinetic Monte Carlo）原理簡單，適應性強，因此在很多情況下都是研究人員的首選。此外，KMC在復雜體系或復雜過程中的算法發展也非常活躍。本文試圖介紹KMC方法的基礎理論和若干進展。KMC方法基本原理在原子模擬領域內，分子動力學(molecular dyna

2、mics, MD)具有突出的優勢。它可以非常精確的描述體系演化的軌跡。一般情況下MD的時間步長在飛秒(s)量級，因此足以追蹤原子振動的具體變化。但是這一優勢同時限制了MD在大時間尺度模擬上的應用。現有的計算條件足以支持MD到10 ns，運用特殊的算法可以達到10 s的尺度。即便如此，很多動態過程，如表面生長或材料老化等，時間跨度均在s以上，大大超出了MD的應用范圍。有什么方法可以克服這種局限呢？當體系處于穩定狀態時，我們可以將其描述為處于維勢能函數面的一個局域極小值（阱底）處。有限溫度下，雖然體系內的原子不停的進行熱運動，但是絕大部分時間內原子都是在勢能阱底附近振動。偶然情況下體系會越過不同勢

3、阱間的勢壘從而完成一次“演化”，這類小概率事件才是決定體系演化的重點。因此，如果我們將關注點從“原子”升格到“體系”，同時將“原子運動軌跡”粗化為“體系組態躍遷”，那么模擬的時間跨度就將從原子振動的尺度提高到組態躍遷的尺度。這是因為這種處理方法擯棄了與體系穿越勢壘無關的微小振動，而只著眼于體系的組態變化。因此，雖然不能描繪原子的運動軌跡，但是作為體系演化，其“組態軌跡”仍然是正確的。此外，因為組態變化的時間間隔很長，體系完成的連續兩次演化是獨立的，無記憶的，所以這個過程是一種典型的馬爾可夫過程(Markov process)，即體系從組態到組態，這一過程只與其躍遷速率有關。如果精確地知道，我們

4、便可以構造一個隨機過程，使得體系按照正確的軌跡演化。這里正確''的意思是某條給定演化軌跡出現的幾率與MD模擬結果完全一致(假設我們進行了大量的MD模擬，每次模擬中每個原子的初始動量隨機給定)。這種通過構造隨機過程研究體系演化的方法即為動力學蒙特卡洛方法(kinetic Monte Carlo, KMC) 1。指數分布與KMC的時間步長在KMC模擬中，構造呈指數分布的隨機數是一個相當重要的步驟。這一節中我們對此進行討論。因為體系在勢能面上無記憶的隨機行走，所以任意單位時間內，它找到躍遷途徑的概率不變，設為。因此在區間內，體系不發生躍遷的概率為類似的，在區間內，體系不發生躍遷的概率

5、為以此類推，當時，在區間內，體系不發生躍遷的概率為因此，當趨于時，體系不發生躍遷的概率為      (1)這一行為類似于原子核的衰變方程。從方程(1)我們可以得到單位時間內體系躍遷概率。從方程(1)的推導過程可以看出體系的躍遷概率是一個隨時間積累的物理量，因此對時間積分到某一時刻必然等于，也即。因此我們立即可以得到 1       (2)是體系處于態時所有可能的躍遷途徑的速率之和，即       (3)對于每個具體的躍遷途徑，

6、上述討論均成立。因此，我們可以定義單位時間內體系進行躍遷的概率為       (4) 單位時間內體系的躍遷概率呈指數分布這一事實說明KMC的時間步長也應是指數分布。因此我們需要產生一個指數分布的隨機數序列。這一點可以非常容易的通過一個(0,1平均分布的隨機數序列轉化得到：從而       (5)最后一步是因為和的分布相同。也可以通過上述步驟從方程(4)得到。計算躍遷速率過渡態理論(TST)決定了KMC模擬的精度甚至準確性。為避開通過原子軌跡來確定的做法(這樣又回到了

7、MD的情況)，一般情況下采用過渡態理論(transition state theory, TST)進行計算 2。在TST中，體系的躍遷速率決定于體系在鞍點處的行為，而平衡態(勢阱)處的狀態對其影響可以忽略不計。如果大量的相同的體系組成正則系綜，則在平衡狀態下體系在單位時間內越過某個垂直于躍遷途徑的縱截面的流量即為。簡單起見，假設有大量相同的一維雙組態(勢阱)體系，平衡狀態下鞍點所在的假想面(對應于流量最小的縱截面)為，則TST給出該體系從組態A遷出到B的速率為 5,6       (6)方程(6)中表示在組態A所屬態空間里對正則系綜的

8、平均。表示只考慮體系從組態A遷出而不考慮遷入A的情況(后一種情況體系也對通過縱截面的流量有貢獻)。根據普遍公式設體系的哈密頓量為，即可分解為動能和勢能，同時設粒子坐標時體系處于組態A。則方程(6)可寫為                      (7)上式中無限小量是為了將函數全部包含進去。最后一項對于函數的系綜平均可以直接通過Metropolis Monte Carlo方法計算出來：計算粒子落在范圍內的次數相對于M

9、etropolis行走總次數的比例。方程(7)最后等于        (8)將上述討論擴展到3維情況非常直接，這里只給出結果，詳細討論請參閱文獻 5:       (9)其中是縱截面方程，代表3維情況中粒子流動方向與截面法向不平行對于計數的影響。簡諧近似下的過渡態理論(hTST)雖然上一節已經給出了TST計算躍遷速率的方法，但是在具體工作中，更多地是利用簡諧近似下的過渡態理論(harmonic TST, hTST)通過解析表達式給出。根據TST，躍遷速率為 3

10、60;      (10)其中為在躍遷中體系在鞍點和態處的自由能之差將上式代入方程(10)，可以得到       (11)hTST認為體系在穩態附近的振動可以用諧振子表示，因此其配分函數是經典諧振子體系的配分函數。分別寫出體系在態和鞍點處的配分函數和：根據Boltzmann公式，       (12)并將配分函數代入，則方程(11)得       (13)方程(13)

11、在通常的文獻上經常可以見到。聲子譜可以通過Hessian矩陣對角化或者密度泛函微擾法(DFPT)求出，而就是的勢壘，可以通過NEB或者drag方法求出。因此，方程(13)保證了可以通過原子模擬(MD或者DFT方法)解析地求出。事實上這個方程有兩點需要注意。首先雖然方程(10)中出現了普朗克常數，但是在最終結果中被抵消了。這是因為TST本質上是一個經典理論，所以充分考慮了統計效應后不會出現 1。其次，方程(13)表明對于每一個躍遷過程，鞍點處的聲子譜應該單獨計算。這樣會大大增加計算量，因此在絕大部分計算中均設前置因子為常數，不隨躍遷過程而變化。具體數值取決于體系，對于金屬而言，一般取 Hz。KM

12、C幾種不同的實現算法 點陣映射 到目前為止，進行KMC模擬的所有理論基礎均已具備。但是前面所進行的討論并沒有聯系到具體的模型。KMC在固體物理中的應用往往利用點陣映射將原子與格點聯系起來。從而將躍遷(事件)具象化為原子格點關系的變化。比如空位(團)/吸附原子(島)遷移等等。雖然與實際情況并不完全一致，但這樣做在很多情況下可以簡化建模的工作量，而且是非常合理的近似。很多情況下體系中的原子雖然對理想格點均有一定的偏離，但是并不太大()，因此這種原子點陣映射是有效的。這種做法的另一個好處是可以對躍遷進行局域化處理。每條躍遷途徑只與其近鄰的體系環境有關，這樣可以極大的減少躍遷途徑的

13、數目，從而簡化計算 1。需要指出的是，這種映射對于KMC模擬并不是必須的。比如化學分子反應爐或者生物分子的生長等等，這些情況下根本不存在點陣。無拒絕方式 KMC的實現方法有很多種，這些算法大致可以分為拒絕(rejection)和無拒絕(rejection-free)兩種范疇。每種范疇之下還有不同的實現方式。本文只選擇幾種最為常用的方法加以介紹。I. 直接法直接法(direct method)是最常用的一種KMC算法，其效率非常高。每一步只需要產生兩個在之間平均分布的隨機數和。其中被用來選定躍遷途徑，確定模擬的前進時間。設體系處于態，將每條躍遷途徑想象成長度與躍遷速率成正比的線段。將這

14、些線段首尾相連。如果落在線段中，這個線段所代表的躍遷途徑就被選中，體系移動到態，同時體系時間根據方程(5)前進。總結其算法如下：1. 根據方程(4)計算體系處于態時的總躍遷速率； 2. 選擇隨機數； 3. 尋找途徑，滿足； 4. 體系移動到態，同時模擬時間前進； 5. 重復上述過程。 需要指出的是，雖然一般步驟4中的根據方程(5)生成，但是如果將其換為并不會影響模擬結果。在文獻5和6中均采用這種方式。II. 第一反應法第一反應法(first reaction method, FRM)在思路上比直接法更為自然。前面說過，對于處于穩態的體系而言，它可以有不同的躍遷途徑可以選擇。每條途徑均可以根據方

15、程(4)給出一個指數分布的"發生時間"，也即從當前算起第一次發生的時間。然后從中選出最小值(最先發生的"第一反應")，體系躍遷到相應的組態，模擬時間相應地前進。總結其算法如下：1. 設共有條反應途徑，生成個隨機數； 2. 根據公式，給出每條路徑的預計發生時間； 3. 找出的最小值； 4. 體系移動到態，同時模擬時間前進； 5. 重復上述過程。 可以看出，這種算法的效率比直接法低下，因為每一步KMC模擬需要生成個隨機數。通常情況下KMC模擬需要步來達到較好的統計性質，如果每一步都需要生成個隨機數，則利用這種方法需要一個高質量的偽隨機數發生器，這一點在比較大

16、時尤為重要。III. 次級反應法次級反應法(next reaction method, NRM)是FRM方法的一種衍生方法，其核心思想是假設體系的一次躍遷并不會導致處于新態的體系對于其他躍遷途徑的舍棄(比如充滿可以發生種化學反應的分子，第一種反應發生并不會造成別的反應物的變化)，這樣體系還可以選擇中的次小值，從而躍遷到態，模擬時間前進。如果這次躍遷還可以滿足上述假設條件，再重復上述過程。理想情況下，平均每一步KMC模擬只需要生成1個隨機數。這無疑會大大提高效率以及時間跨度。但是實際上NRM的假設條件很難在體系每次躍遷之后都得到滿足，在固體物理的模擬中尤其如此，因此其應用范圍集中于研究復雜化學環

17、境下的反應過程。試探-接受/拒絕方式 這一大類算法雖然在效率上不如直接法，但是它們所采用的試探-接受/拒絕在形式上更接近Metropolis MC方法，而且可以很方便的引入恒定步長，即固定。因此有必要進行詳細的介紹。IV. 選擇直接法選擇直接法在決定體系是否躍遷方面和Metropolis MC方法形式上非常相像，均是通過產生隨機數和預定的閾值比較決定事件是否被采納。具體算法如下：1. 設共有條反應途徑，選擇反應速率最大值，設為。生成在均勻分布的隨機數； 2. 設； 3. 如果 < ，則體系躍遷至新態，否則保持在態； 4. 模擬時間前進； 5. 重復上述過程。 這種方法的長處在于

18、每一步只需要生成一個隨機數。但是缺點也很明顯，對于反應速率相差太大，尤其是只有一個低勢壘途徑(與其他途徑相比過大)的體系來講，這種方法的效率會非常低下。某些情況下，這種低效率問題可以通過如下方法改進：將全部途徑按照的大小分為幾個亞組，每個亞組選定一個上限。但是這一步驟在整個KMC模擬過程中可能需要重復很多次，因此并不能完全解決問題。事實上低勢壘在KMC中是個普遍的問題。這一點在后面還要簡要提及。V. 恒定步長法與上述四種方法不同，恒定步長法(constant time step method, CTSM)中體系的前進時間是個給定的參數citedawnkaski。在理想情況下，CTSM與直接法效

19、率相同，每一步只需產生兩個隨機數。具體算法如下：1. 給定恒定時間步長； 2. 將所有途徑(共有個)設為長度恒為的線段，生成在均勻分布的隨機數，選擇途徑； 3. 生成在均勻分布的隨機數，如果 < ，則體系躍遷至新態，否則保持在態； 4. 模擬時間前進； 5. 重復上述過程。 實際模擬中，需要滿足(1)小于(見"第一反應法")，以及(2)對于最大的途徑，接受率大致在0.5。其中第一個條件保證了所有的遷移途徑發生概率都小于1，第二個條件則保證體系演化的效率不會過于低下。CTSM是非常行之有效的一類KMC算法，但是選擇時需要特別的注意以保證效率。決定于具體體系以及模擬溫度。

20、這在一定程度上增加了CTSM的實現及使用難度。低勢壘問題前面已經指出，低勢壘的途徑需要特別注意。如果體系在演化過程中一直存在著勢壘較其他途徑低很多的一個或幾個途徑，會對模擬過程產生不利的影響。這個問題被稱之為低勢壘問題。低勢壘途徑對于KMC模擬最直接的影響就是大大縮短了模擬過程所涵蓋的時間跨度。這一點可以從方程(5)中看出。更為深刻的影響在于，這些由低勢壘的途徑聯系起來的組態會組成一個近似于封閉的族。體系會頻繁的訪問這些態，而其他的對于體系演化更為重要的高勢壘途徑被選擇的概率非常低，這顯然會降低KMC的模擬效率。例如，吸附原子在高指數金屬表面擴散，其沿臺階的遷移所對應的勢壘要遠低于與臺階分離的

21、移動。這樣，KMC模擬的絕大部分時間內吸附原子都在臺階處來回往復，而不會選擇離開臺階在平臺上擴散。這顯然不是我們希望看到的情形。一種解決辦法是人為地將這些低勢壘加高以降低體系訪問這些組態的幾率，但是無法預測這種干擾是否會造成體系對于真實情況的嚴重偏離。另一種選擇是利用NRM或者CTSM進行模擬，但是其效果如何尚待檢測。如果考察體系的勢能面，這類低勢壘的途徑一般處在一個"超勢阱"之中。體系在這個超勢阱中可以很快的達到熱平衡，所需時間要短于從其中逸出的時間。如果可以明確的知道超勢阱所包含的組態以及從超勢阱逸出的所有途徑，我們就可以按照Boltzmann分布合理的選擇其中一條途徑

22、，使得體系向前演化。但是如何確定哪些組態包含在超勢阱之中以及體系是否已在其中達到熱平衡本身就是兩個難題。對于第一點，Mason提出可以利用Zobrist密鑰法標定訪問過于頻繁的組態 7；Novotny則提出通過建立及對角化一個描述體系在這些組態間演化的傳遞矩陣來解決第二點 8。對這個問題的詳細討論已超出了本文的討論范圍，請參閱文獻7以及8。實體動力學蒙特卡洛方法 OKMC上述的KMC都假設任何時候原子均處于其理想點陣格子上。但是很多情況下這種點陣映射是無效的，比如間隙原子或者位錯。這類結構缺陷的運動在材料的輻射損傷和老化過程中扮演著非常重要的角色。而且與單個原子或者空位的運動相比，這類缺陷的運

23、動時間跨度更長，也更為復雜，比如間隙原子團和空穴的湮沒，間隙原子團的解構/融合，或者位錯的攀移/交滑移等等。傳統的KMC算法很難有效的處理這類問題，一方面是因為時間跨度太大，另一方面這類缺陷各自均可視為獨立的實體(object)，其運動更近似于系統激發，因此單個或幾個原子運動的積累效果很多情況下并不能有效地反應這些實體的整體運動。實體動力學蒙特卡洛方法(Object kinetic Monte Carlo, OKMC)就是為了處理這類問題而被提出的。OKMC在算法上與普通的KMC完全一樣。需要注意的地方是在OKMC中并不存在原子點陣。所有的實體在一個真空的箱子中按照其物理實質離散化運動，比如位

24、錯環的最小移動距離是其Burgers矢量大小，方向則為Burgers矢量方向；空位的移動距離為第一近鄰或第二近鄰的原子間距，等等。模擬過程中我們需要追蹤該實體的形心，從而決定其位置、移動距離等等。此外，OKMC中對于躍遷速率的確定也和普通的KMC有所區別。本文前面已經指出，可以表達為的形式。普通的KMC假定為常數，不同途徑的由決定。但是在OKMC的模擬中，的直接確定非常困難，因此一般的策略是對于特定的事件(包括實體自身的運動以及不同實體間的反應等)，躍遷勢壘保持恒定，而將前置因子視為實體規模(所包含的原子/空位數目)的函數，通過MD模擬得出，一般而言可以表示為形如的表達式，其中和是擬合參量，是

25、實體規模。最后需要注意的是在OKMC的模型中，實體有空間范圍，因此需要一個額外的參數來表征其空間半徑(假設為球形分布，否則的數目多于一個)。在模擬不同實體間的反應時，需要特別考慮其形心的間距，如果小于"反應距離"，即，反應一定進行，否則認為兩個實體互相獨立。Domain利用OKMC研究了Fe-Cu合金的輻射損傷9，在模擬中考慮了間隙原子(空位)的聚合、間隙原子(空位)團的發射、間隙原子-空位湮沒、空位團對雜質的捕獲、表面對于空位(團)的捕獲、甚至輻射轟擊引起的間隙原子(空位)萌生、增殖等等事件。從中可以看出，對于OKMC，一個棘手的問題是需要預先想到所有的事件。此外，OKM

26、C所需要的所有參量基本上不可能通過原子模擬直接獲得，人為的設定參數不可避免。這些參數會在多大程度上決定OKMC的準確程度無法預先得知。需要根據現有的實驗數據進行修改、調試。這些困難都限制了OKMC的普及。但是如前所述，這種方法可以有效地進行大尺度的時間(天)和空間模擬(m以上)，而且對于缺陷的描述更為直接和符合直觀，因此在材料研究中同樣占有重要的地位。KMC的若干進展等時蛙跳算法(-leap KMC) 引入這類算法前，我們先簡要介紹兩個常用的離散分布：泊松分布(Poisson Distribution, PD)以及二項式分布(Binomial Distribution, BD)。泊松

27、隨機數定義為給定事件發生率以及觀測時間下事件發生的數目。如果用代表給定的發生數目，則恰好等于的概率是一個泊松分布：       (14)也即如果產生一個泊松隨機數序列，則這個序列符合泊松分布PD。需要指出，是無界的，范圍是任意非負整數。與其類似，二項式隨機數定義為重復次獨立的成功率均為的伯努利實驗的成功數。如果給定成功數，則恰好等于的概率是一個二項式分布：       (15)為了和本文中的標號一致，我們將躍遷的成功率表示為，將方程(15)重新寫為  &#

28、160;    (16)與PD不同，BD中的是有界的，為0到之間的任意整數。可以看出，如果將這兩種隨機數理解為給定躍遷路徑(發生率為)在一定的時間步長()內發生的次數，則可以立即運用于粒子數空間內的KMC中，其時間范圍可以得到很大提高。這就是等時蛙跳算法-leap KMC 10,11。-leap KMC方法最早由Gillespie提出，通過PD方程(14)，在給定時間步長下決定每個躍遷途徑發生的次數，然后將體系移到這些躍遷累計發生后產生的新態。因為每一步模擬體系不止發生一次躍遷，所以模擬的速度可以大大加快。我們以多種反應物在化學反應爐中的演化為例加以詳細說明。設在

29、爐內共有種分子，在時刻各自的個數為，則在粒子數空間中構成一個矢量，或稱為一個組態。總共有種反應路徑。對于給定的，反應速率是占據態的函數。此外，我們單獨定義一個矢量，其中由通過反應而得，即。因此的元素代表反應所引起的種分子的數目變化。由此建立算法如下：VI. PD-leap KMC 101. 給定恒定時間步長； 2. 對于每條反應途徑按照方程(14)生成泊松隨機數序列，按照模擬步數從序列中找出每種反應發生的次數； 3. 按照更新體系； 4. 模擬時間前進； 5. 重復上述過程。 Gillespie仔細考慮了的選擇條件，稱為蛙跳條件(leap condition)：  

30、0;    (17)其中 如前所述，沒有上限，因此即使滿足方程(17)，在模擬過程中也可能會出現某種分子總數為負數的情況，這顯然不符合實際，也是PD-leap KMC的一個弱點。Tian和Burrage提出可以用二項式分布BD取代PD，因為有上限，所以可以有效的解決這個問題。此外，他們對于某種分子參與多種反應的情況也進行了考慮，從而提高了-leap KMC的穩定性和普適性。其算法如下：VII. BD-leap KMC 11 · 給定恒定時間步長，滿足； · 對于每條反應途徑按照方程(16)生成二項式隨機數序列，按照模擬步數從

31、序列中找出每種反應發生的次數；如果有某種分子同時參與了和，則首先生成然后通過確定的發生次數； · 按照更新體系； · 模擬時間前進； · 重復上述過程。 步驟1、2中出現的是參與反應的各類分子的個數的最小值，即此外Gillespie，Tian和Burrage還考慮用預測時刻體系狀態的方法來進一步提高精度。具體請參閱文獻10,11。如果-leap算法和OKMC結合起來可以進一步加大模擬的時間尺度，但是目前還沒有這方面工作的介紹。基于即時動態分析的KMC方法(on-the-fly KMC)到目前為止，所有的KMC都是在模擬之前建立好所有可能的躍遷途徑。但是實際上&qu

32、ot;所有"是很難達到的目標。因為很多途徑遠離一般的直覺，而且在演化過程中體系有可能尋找到新的途徑。因此，躍遷途徑應該隨著體系的演化而不斷更新，是動態的過程。Henkelman和Jönsson將途徑搜索和KMC結合起來，提出了即時動態的KMC方法 on-the-fly KMC 12：在每一個穩態(勢阱)處，選定一個激活原子(一般是近鄰不飽和的原子)，在以其為中心的局部區域內引入呈高斯分布的隨機位移，即加入擾動，然后利用dimer方法 13尋找所有可能的躍遷途徑。建立起即時的途徑庫之后再通過普通KMC算法進行模擬。顯然，這種方法的計算量非常大，需要一個有效的標識方法來識別所有

33、已經遇到過的途徑以避免重復計算。Trushin提出可以利用包括至激活原子第三殼層的所有格點(順時針排列)的占據與否(分別標記為1和0)來構建二進制數，從而根據始態和終態的標號來唯一地標識某條途徑 14，例如，激活原子標為"1"，其第一殼層的原子標記為"2","3",""，依此類推，然后將原子的標號""作為二進制的數位，這樣，每一個穩態都有唯一的一個二進制數與之對應。雖然仍不完善，但是這種方法具有非常清晰的邏輯結構，具有良好的擴展性。和KMC方法一般情況下KMC的大部分時間花費在選擇途徑上。如果采

34、用普通的方法，即循環疊加直至從而選擇，這種情況下計算用時與途徑數目呈線性增長，即算法。按照二叉樹安排不同數目的之和可以改進到 15：將所有作為樹葉(不足2整數次冪的葉子由0填補)，每兩片葉子之和作為父節點，依次類推直至樹根。一株二叉樹構建完畢后，生成一個隨機數，由樹根開始尋找，若不大于左子節點，沿左分支向下尋找；否則設，沿右分支向下尋找，直至樹葉，體系按途徑演化。Slepoy和Thompson等進一步提出分流-拒絕(composition-rejection, CR)方法以實現搜索用時與途徑總數無關的算法 16：(1)先找出和，按照()將條途徑分為個組，(2)然后生成隨機數，按照上述二叉樹尋找所落入的組別，(3)再生成兩個隨機數和，設，其中為該組中包含的途徑數，如果，則選擇途徑，否則重復步驟(3)，直至有一條途徑被選中為止。可以看出，CR算法雖然搜索速度很快，但是每一步KMC需要產生至少4個隨機數(用于確定前進時間)，因此需要高質量的隨機數發生器。不過對于躍遷途徑復雜的體系演化而言，CR的效率無疑是很有吸引力的。 1 A.F. Voter, it Radiation Effects in Solids (Springer 2006) p. 1