MD與KMC的耦合模擬研究與實現(xiàn)

李建江，魏 鵬，楊少峰，賀新福，胡長軍

(1.北京科技大學(xué)計算機(jī)科學(xué)與技術(shù)系 北京 海淀區(qū) 100083；2.中國原子能科學(xué)研究院 北京 房山區(qū) 102413)

材料輻照損傷是當(dāng)前材料領(lǐng)域和計算機(jī)領(lǐng)域研究的熱點之一。在輻照條件下的材料模擬中，入射粒子與材料中的晶格原子相互碰撞并傳遞能量，晶格原子獲得能量后將離開晶格點陣并引發(fā)級聯(lián)碰撞，形成點缺陷及團(tuán)簇。在輻射條件下，點缺陷將復(fù)合、遷移、聚集成團(tuán)，最終導(dǎo)致材料微觀結(jié)構(gòu)演化以及宏觀力學(xué)性能退化。多尺度模擬是研究材料輻照效應(yīng)的有效手段，可從不同的時間、空間尺度研究材料的損傷機(jī)理，預(yù)測材料的損傷程度。對材料輻照損傷的研究是一個跨越飛秒和年的大時間尺度以及跨越納米和米的大空間尺度問題。實現(xiàn)從原子尺度到宏觀尺度的耦合模擬，是一個龐大而復(fù)雜的工程，而每一個尺度又可劃分為多個更細(xì)的尺度，分別對應(yīng)不同的模擬方法。本文只研究微觀多尺度上的耦合模擬。在微觀尺度上，分子動力學(xué)(MD)和動力學(xué)蒙特卡洛(KMC)為常用的兩種經(jīng)典的材料模擬方法。MD方法[1]通過求解動量守恒、質(zhì)量守恒和能量守恒這三大守恒定律中的運動方程，來描述系統(tǒng)某一時刻的狀態(tài)，再通過積分運算獲得系統(tǒng)的性質(zhì)。KMC方法[2]是一種基于概率統(tǒng)計理論的方法，通過計算事件的發(fā)生概率和使用隨機(jī)數(shù)獲得系統(tǒng)組態(tài)的躍遷。由于KMC方法的計算對象是組態(tài)而非原子，對原子的特征描述維度(如：速度、力、質(zhì)量等)比MD少，故其模擬的時間尺度和空間尺度都比MD大。雖然，MD方法對系統(tǒng)的描述比KMC方法更加準(zhǔn)確，但由于模擬的時間步長太短(約為10-15s)，使得計算量和計算開銷巨大，對計算機(jī)的運行性能要求也較高[3]。與之相反，KMC方法則可突破MD方法在時間尺度上的限制，可達(dá)到秒或者更長的模擬時間，此方法可以對系統(tǒng)的動力學(xué)特征進(jìn)行描述，但由于它使用隨機(jī)數(shù)而降低了對系統(tǒng)描述的準(zhǔn)確性。因此，目前常用的做法是將MD方法和KMC方法耦合起來，在微觀尺度模擬材料輻照損傷中缺陷的形成及演化過程。

多尺度模擬的關(guān)鍵之一在于各層次之間信息傳遞的準(zhǔn)確性，同一系綜下多尺度模擬方法可以實現(xiàn)信息的傳遞，提高模擬精度。多尺度模擬方法在材料輻照損傷研究方面已開展大量工作，并取得了很好的結(jié)果，如文獻(xiàn)[4]針對材料輻照損傷問題，開發(fā)了預(yù)測性計算工具來模擬輻照材料的微觀結(jié)構(gòu)變化。描述了如何使用來自數(shù)千個MD模擬的碰撞級聯(lián)的統(tǒng)計平均值來為KMC模擬提供輸入，該模擬可以處理更大的尺寸，更多的缺陷和更長的持續(xù)時間。文獻(xiàn)[5]結(jié)合分子動力學(xué)方法和動力學(xué)蒙特卡羅方法，研究了單個粒子入射硅引起的位移損傷缺陷的產(chǎn)生和演化過程。文獻(xiàn)[6]提出了自適應(yīng)動力學(xué)蒙特卡洛方法，該方法將KMC方法的簡單性與基于MD的鞍點搜索算法相結(jié)合來模擬亞穩(wěn)系統(tǒng)。文獻(xiàn)[7]拓展了多尺度耦合模擬工具M(jìn)aMiCo，該工具將連續(xù)流體動力學(xué)求解器與離散粒子動力學(xué)耦合，引入了對分子動力學(xué)模擬集合的采樣，最后通過使用超級計算機(jī)Shaheen II的多達(dá)65 536個計算核心來驗證了分子連續(xù)模擬的并行可擴(kuò)展性。文獻(xiàn)[8]采用KMC獲得富銅析出物，然后將KMC的模擬結(jié)果傳給相場，用相場方法模擬富銅析出物的長大過程，同時采用MD方法對材料的結(jié)構(gòu)進(jìn)行檢查，確保模擬過程中材料結(jié)構(gòu)的一致性，有效地模擬了鐵中富銅析出物的粗化過程。文獻(xiàn)[9]通過MD的模擬和KMC順序耦合的方法(MD的模擬結(jié)果作為KMC模擬的輸入)分別模擬了銅和鐵的輻照損傷過程。文獻(xiàn)[10]利用MD方法計算系統(tǒng)的能量，將計算獲得的能量用于KMC計算來研究鐵中氦-空位團(tuán)簇的增長和收縮機(jī)制。

本文通過對動力學(xué)方法和動力學(xué)蒙特卡洛方法的特點分析，將分子動力學(xué)方法和動力學(xué)蒙特卡洛方法進(jìn)行耦合，可進(jìn)一步突破現(xiàn)有時空尺度的限制，從而實現(xiàn)長時間、多空間尺度的原子層次直接模擬。

1 MD與KMC耦合模擬

本文中，MD和KMC分別模擬材料輻照損傷的兩個階段，MD模擬缺陷的產(chǎn)生及團(tuán)簇的形成過程為KMC模擬提供初始的原子排布和缺陷信息，KMC模擬完成缺陷從小到大演化的模擬過程。

1.1 MD方法與KMC方法原理

MD方法[1]的基本思想是把所有物質(zhì)都看成是由電子、原子或分子組成的粒子系統(tǒng)，每個粒子在其他粒子所提供的經(jīng)驗勢場的作用下運動，這種運動被認(rèn)為遵循牛頓運動定律。若已知體系中各粒子的初始位置，通過勢能函數(shù)計算出粒子所對應(yīng)的勢能，利用勢能求得粒子在體系中所受的力，再通過求解牛頓運動方程，計算出粒子的加速度，從而可以通過對時間的積分計算得到粒子在任意時刻的速度以及位置。如果時間足夠長，即經(jīng)過的循環(huán)迭代次數(shù)足夠多，最終獲得材料的宏觀性能就越明顯。

KMC模擬方法[2]是一種通過構(gòu)造隨機(jī)過程來模擬體系長時間演化的方法。KMC一般用來模擬由組態(tài)躍遷主導(dǎo)的體系的長時間(通常是秒或者更長)的演化。

本文中的KMC模擬程序采用最常用的一種KMC算法直接法(direct method)，效率非常高，屬于無拒絕KMC。每一步只需要產(chǎn)生兩個平均分布在(0,1]之間的隨機(jī)數(shù)r1和r2。其中，r1被用來選定躍遷途徑，r2確定模擬的前進(jìn)時間。設(shè)體系處于態(tài)i，將每條躍遷途徑j(luò)想象成長度與躍遷速率kij成正比的線段。將這些線段首尾相連，如果r1×ktotal落在線段jk中，則這個線段所代表的躍遷途徑j(luò)k就被選中，體系移動到態(tài)jk，同時由隨機(jī)數(shù)r2計算體系前進(jìn)時間。

1.2 耦合模擬流程

由于MD和KMC這兩種方法的計算原理和過程不同，所關(guān)注的原子屬性也不一樣。其中，在MD方法中原子的位置(三維坐標(biāo))可以為模擬空間中的任意位置。而在KMC模擬中，由于KMC對間隙原子的躍遷情況計算復(fù)雜、耗時較多，在目前的系統(tǒng)中暫時不對間隙原子做KMC模擬。且由于KMC關(guān)注的是體系從一個狀態(tài)到另一個狀態(tài)的“躍遷”，故在KMC的模擬體系中，原子位置為其晶格的格點位置，原子不會出現(xiàn)在晶格中除格點外的其他位置上。因為存在上述這些不同，所以，在MD模擬結(jié)束后，MD模擬為KMC模擬提供的原子排布和缺陷信息無法在KMC程序中直接使用。為了解決這一問題，且遵循軟件開發(fā)高內(nèi)聚低耦合準(zhǔn)則，本文開發(fā)了一個基于MD和KMC耦合模擬的耦合程序。它的主要功能是過濾間隙原子，實現(xiàn)原子位置重置(置于對應(yīng)格點)，以及完成數(shù)據(jù)的格式轉(zhuǎn)換工作。MD和KMC耦合模擬的示意圖如圖1所示。其中，t1為MD模擬的當(dāng)前累計模擬時間，tMD為MD模擬的總時間，t2為KMC模擬的當(dāng)前累計模擬時間，tKMC為KMC模擬的總時間。

圖1 MD和KMC耦合模擬示意圖

1.3 耦合程序難點分析

要想實現(xiàn)該耦合程序，必須要知道體系中數(shù)據(jù)所表示的含義，了解模擬對象的晶格結(jié)構(gòu)和空間排布。本文研究是以鐵為基質(zhì)，混有少量銅雜質(zhì)的金屬合金RPV鋼(RPV鋼中的雜質(zhì)除了銅，還有鎳、錳等，由于其他雜質(zhì)的含量少于銅，因此本文只研究銅雜質(zhì)在輻照損傷中對RPV鋼的影響)。雖然，純銅固體的晶格結(jié)構(gòu)為面心立方晶格(face-centered cubic,FCC)結(jié)構(gòu)，但由于合金中銅含量極少，且在研究范圍內(nèi)的銅團(tuán)簇半徑很小，因此整個體系仍是體心立方晶格(body-centered cubic, BCC)結(jié)構(gòu)。圖2就是BCC的結(jié)構(gòu)示意圖，圓球表示格點。

圖2a展示的就是體心立方晶格的一個堆積，圖2b為對應(yīng)的結(jié)構(gòu)單元，圓球所在的位置是晶格點位置。當(dāng)出現(xiàn)多個原子對應(yīng)一個格點時，只有一個原子是正常原子，其他原子稱為間隙原子。如果格點上沒有原子，則該格點為空位。如有間隙原子出現(xiàn)，則在體系內(nèi)必然存在同樣數(shù)目的空位，空位和間隙原子都屬于缺陷。在MD模擬中原子可以出現(xiàn)在晶格的任意位置，而在KMC模擬中原子都在圓點位置，不會發(fā)生偏移。因此，要想實現(xiàn)MD和KMC的順序耦合，首先必須確定MD中每個原子屬于哪個格點，并進(jìn)行原子位置重置，將原子置于對應(yīng)的格點位置。

圖2 體心立方晶格的堆積與結(jié)構(gòu)單元

目前，在MD的模擬中可以處理間隙原子，但在KMC模擬中不能處理間隙原子(計算太過復(fù)雜)，所以MD和KMC耦合模擬的耦合程序要對所有原子進(jìn)行過濾。需要注意的是，由于雜質(zhì)原子在模擬過程中對結(jié)果的影響較大，所以在處理除去間隙原子的過程中本文采取雜質(zhì)原子優(yōu)先保留的策略，即如果一個格點同時對應(yīng)一個銅原子和一個鐵原子，將會選擇鐵原子為間隙原子，并將其從體系中除去。

由上可知，原子的位置重置和間隙原子過濾都需要先判定原子屬于哪個格點。如何確定原子對應(yīng)的格點，以及判定原子的類型(間隙原子、空位或正常原子)是耦合程序?qū)崿F(xiàn)中的難點。在第2章中將詳細(xì)介紹這些問題的解決方法。

2 基于Wigner-Seitz原胞法的SD算法實現(xiàn)

2.1 Wigner-Seitz原胞缺陷分析法

判斷偏離格點位置的原子屬于哪個格點是實現(xiàn)原子位置重置和過濾間隙原子的關(guān)鍵。在可視化軟件OVITO[11]中就采用了Wigner-Seitz原胞缺陷分析法對可視化的原子進(jìn)行分析。Wigner-Seitz原胞是晶格中比較對稱的一種原胞，以一個格點為原點，作原點與其他格點連接的中垂面(二維空間中為中垂線)，由這些中垂面(二維空間中為中垂線)所圍成的最小體積(二維空間中為面積)，圖3和圖4分別為在二維空間和三維空間中BCC結(jié)構(gòu)材料的Wigner-Seitz原胞。其中，圖3中的加粗區(qū)域為格點的Wigner-Seitz原胞區(qū)域，該區(qū)域為等邊六邊形；圖4中每個格點的Wigner-Seitz原胞區(qū)域為一個截角八面體。

圖3 二維空間中BCC結(jié)構(gòu)材料的Wigner-Seitz原胞

圖4 三維空間中BCC結(jié)構(gòu)材料的Wigner-Seitz原胞

Wigner-Seitz原胞缺陷分析法的思想是：由于每個Wigner-Seitz原胞中只包含一個原子，即它的原點處的格點對應(yīng)的原子如果出現(xiàn)多個，則其他原子就是間隙原子，就是MD和KMC中間處理過程需過濾掉的原子。如果某個Wigner-Seitz原胞中不含任何原子，則對應(yīng)格點即為空位。在材料模擬和分析中通常也將空位看成是一種特殊的原子。

2.2 最短距離算法原理及其實現(xiàn)

可視化軟件OVITO[11]的處理過程需要初始原子位置文件和分析時刻的原子位置文件這兩個文件，通過位置對照，利用Wigner-Seitz原胞缺陷分析法對要分析的原子中的缺陷進(jìn)行識別。考慮到在MD和KMC耦合過程中若用OVITO中的方法實現(xiàn)Wigner-Seitz原胞缺陷分析法，需要存儲體系的初始原子分布，且本文的研究對象為BCC結(jié)構(gòu)的固體材料，原子的原始排布(完美晶格)已確定，故本文依據(jù)Wigner-Seitz原胞法[11]的原理，提出了通過距離判斷尋找原子對應(yīng)晶格點的算法，以此來實現(xiàn)間隙原子、空位和正常原子的判定，該算法被稱為最短距離(SD)算法。在SD算法中不需要存儲系統(tǒng)的初始原子分布，從而大大節(jié)約了內(nèi)存開銷。

在整個物理空間中，距離任意一點最近的格點為其所在的Wigner-Seitz原胞對應(yīng)的原點。因此， SD算法的核心思想就是利用距離找出距離該原子最近的格點即可。圖5為MD模擬中的原子處于非標(biāo)準(zhǔn)格點處的情況，本文給每個格點位置(空心圓圈)編號為1、2、3、4、5、6、7、8、9，其中實心圓圈表示原子。

圖5 MD模擬中原子處于非標(biāo)準(zhǔn)格點處的情況

采用SD算法計算原子與格點的映射關(guān)系的過程如下：體系中的每個原子都一定處于某個確定的晶格中，則距離原子最近的格點必定是它所在的晶格立方體的體心處格點或頂點處的格點(簡稱為體心格點、頂點格點)。若某個頂點格點在各個維度上的坐標(biāo)與對應(yīng)的原子坐標(biāo)相減的絕對值都是最小，則該格點為距離原子最近的頂點格點(假設(shè)為圖5中的2號格點)，計算兩者間的距離設(shè)為d1；然后計算晶格中的體心格點與原子的距離設(shè)為d2；如果d1＞d2，則原子對應(yīng)的格點為9號格點(體心格點)，否則對應(yīng)的格點為2號格點(頂點格點)。如果有多個原子同時對應(yīng)9號格點，則分為兩種情況處理：一種情況是多個原子中有鐵原子也有兩個以上雜質(zhì)原子，此時將鐵原子全部視為間隙原子。按照遍歷順序，將第一個雜質(zhì)原子保留，其他視為間隙原子。另一種情況，只有鐵原子或只有雜質(zhì)原子，則按照遍歷順序，保留第一個處理的原子，其他都視為間隙原子。圖6所示為兩個原子同時對應(yīng)9號格點的情況，即d3＞d4且d1＞d2。

圖6 一個格點(9號)同時對應(yīng)兩個原子的情況

在正則系綜下為了保證整個體系中模擬的原子數(shù)目不變，在MD模擬和KMC模擬中都采用了周期性邊界條件(periodic boundary conditions, PBC)[12]來處理模擬盒邊界上原子在空間里的運動問題。在處理原子與格點映射以及除去原子的過程中，周期性邊界條件也是必須要考慮的一個問題。在整個耦合模擬的過程中，模擬盒的大小是固定不變的，采用周期性邊界條件使得在原子運動過程中當(dāng)有一個或幾個原子從模擬盒的邊界跑出時，就從相反的界面進(jìn)入一個或幾個原子，從而保證在體系中原子的數(shù)目是恒定的。圖7分別為二維空間和三維空間中周期性邊界條件的處理情況。其中，黑色點和白色點分別表示兩種不同的原子，正方形和正方體用來表示模擬盒。算法1為SD算法的偽代碼實現(xiàn)。

圖7 二維空間和三維空間中的周期性邊界條件示意圖

算法1 最短距離算法

Input： 模擬盒大小、晶格常數(shù)、原子坐標(biāo)

Output： input_KMC.xyz，interstitialatom.xyz

Dim vertex[n1][n2][n3]center[n1][n2][n3]d1d2As INTEGER

//定義頂點格點計數(shù)變量與體心格點計數(shù)變量

//定義原子分別到頂點格點的距離和體心格點的距離

//n1、n2、n3分別為體系x、y、z方向上的晶胞數(shù)(BCC材料每個晶胞包含一個頂點格點和一個體心格點)。

// input_KMC.xyz為KMC的輸入，內(nèi)容為體系的格點坐標(biāo)和對應(yīng)的原子類型；

// interstitialatom.xyz為過濾出的間隙原子的坐標(biāo)；Begin：

Read 模擬盒大小、晶格常數(shù)、原子坐標(biāo)

For all 銅原子(x,y,z)

計算與原子最近的頂點格點(xv,yv,zv)之間的距離d1;

計算與原子最近的體心格點(xc,yc,zc)之間的距離d2;

Ifd1＞d2then

If center[xc][yc][zc]＞0 then

將銅原子(x,y,z)的信息輸出到interstitialatom.xyz；

Else

center[xc][yc][zc]=center[xc][yc][zc]+1;

將原子類型、格點坐標(biāo)信息輸出到input_KMC.xyz；

End if

Else

If vertex [xv][yv][zv]＞0 then

將銅原子(x,y,z)的信息輸出到interstitialatom.xyz；

Else

vertex [xv][yv][zv]=vertex [xv][yv][zv]+1;

將原子類型、格點坐標(biāo)信息輸出到input_KMC.xyz；

End if

End if

End for

For all 鐵原子(x,y,z)

計算與原子最近的頂點格點(xv,yv,zv)之間的距離d1；

計算與原子最近的體心格點(xc,yc,zc)之間的距離d2；

Ifd1＞d2then

If center[xc][yc][zc]＞0 then

將鐵原子(x,y,z)的信息輸出到interstitialatom.xyz；

Else

center[xc][yc][zc]=center[xc][yc][zc]+1;

將原子類型、格點坐標(biāo)信息輸出到input_KMC.xyz；

End if

Else

If vertex [xv][yv][zv]＞0 then

將鐵原子(x,y,z)的信息輸出到interstitialatom.xyz；

Else

vertex [xv][yv][zv]=vertex [xv][yv][zv]+1;

將原子類型、格點坐標(biāo)信息輸出到input_KMC.xyz；

End if

End if

End for

Output 格點的坐標(biāo)及格點的原子類型End

3 實驗結(jié)果與分析

本文中的所有測試都是在新一代超級計算集群“元-Era”上進(jìn)行的。該集群總共擁有270臺曙光CB60-G16雙路刀片，CPU整體性能達(dá)到120.96 T flops。每臺刀片計算節(jié)點配置2個Intel E5-2680V2(Ivy Bridge|10 C|2.8 GHz)處理器和64 GB DDR3 ECC 1866 MHz內(nèi)存，每個處理器含10個處理核心。操作系統(tǒng)為RedHat Linux 4.4.7-3版本，編譯器為Intel C++ 2013_sp1版，MPI環(huán)境為MVAPCH2-intel1.9版本。中間程序的編程語言為C語言，不需要其他外部庫的支持，可在任何一臺具有C編譯器的linux服務(wù)器上運行。

3.1 SD算法驗證

為了驗證本文提出的SD算法實現(xiàn)Wigner-Seitz原胞缺陷分析的正確性，本文使用最短距離法處理后獲得的結(jié)果與OVITO軟件的處理結(jié)果進(jìn)行了對比。由于SD算法中在處理多個原子映射到同一個格點時，采取了雜質(zhì)原子優(yōu)先保留和同等條件下先遍歷先保留的策略，所以SD算法處理后產(chǎn)生的間隙原子和OVITO處理后產(chǎn)生的間隙原子可能會出現(xiàn)位置相同但類型不同的情況，所以在評估SD算法的正確性的時候本文只對間隙原子的數(shù)量和對應(yīng)格點的位置進(jìn)行對比分析。由于體系中的空位是間隙原子進(jìn)入具有原子的原胞造成的，所以間隙原子的數(shù)目等于空位的數(shù)目。

表1和表2表示在含2×106個原子的MD輸出文件中，設(shè)置間隙原子的個數(shù)為4、20、50、100、200、400時，分別使用本文實現(xiàn)的MD與KMC耦合模擬的中間處理程序中和OVITO中的Wigner-Seitz原胞缺陷分析法處理后得到的間隙原子數(shù)和空位數(shù)。

表1 SD算法處理結(jié)果和用OVITO中的Wigner-Seitz原胞缺陷分析法處理結(jié)果對比(間隙原子數(shù)量)

表2 SD算法處理結(jié)果和用OVITO中的Wigner-Seitz原胞缺陷分析法處理結(jié)果對比(空位數(shù)量)

3.2 MD和KMC耦合模擬RPV鋼的輻照損傷實驗

本文實現(xiàn)的用MD和KMC耦合的方法模擬在溫度為600 K時，F(xiàn)e-0.3Cu(摩爾分?jǐn)?shù)，%)的RPV (reactor pressure vessel)鋼的輻照損傷實驗。實驗用MD方法模擬RPV鋼中原子級聯(lián)碰撞產(chǎn)生空位的過程；用KMC方法模擬團(tuán)簇的形核與長大過程。模擬盒大小 為 100a×100a×100a(a為晶格常數(shù)，大小為0.285 53 nm)，原子總數(shù)為2×106。模擬過程中采用周期性邊界條件，計算所用勢函數(shù)為EAM勢。

本文用并行分子動力學(xué)模擬軟件LAMMPS[13]使用MD方法模擬RPV鋼的輻照損傷產(chǎn)生缺陷(主要是空位和間隙原子)的過程[14]，模擬過程時間均是程序模擬材料演化的時間，而非計算機(jī)的運行時間，模擬時間為0.37 ps，在MD模擬結(jié)束時進(jìn)行能量最小化處理。圖8為MD模擬結(jié)束后的原子可視圖，其中，灰色圓點表示銅原子，白色圓點表示空位，黑色部分為鐵原子。圖8b是圖8a中空位聚集部分對應(yīng)的放大圖。

圖8 MD模擬級聯(lián)碰撞結(jié)果可視圖

從圖8可以看出，在Fe-Cu合金中，通過級聯(lián)碰撞產(chǎn)生了空位缺陷，缺陷的位置雖然比較集中(高能快中子打入造成的)，但還都是離散的，并沒有形成團(tuán)簇。通過耦合程序的處理，對MD輸出的原子進(jìn)行清洗，除去間隙原子，并將原子映射到對應(yīng)的格點位置，最后對數(shù)據(jù)進(jìn)行格式轉(zhuǎn)換處理，將體系格點處的原子或缺陷的排布傳遞給相應(yīng)KMC模擬程序。程序接收原子和缺陷信息，并進(jìn)入對缺陷在輻照條件下的形核與長大過程進(jìn)行模擬。本文中KMC的模擬時間為1.25×10-2mcs (mcs為KMC模擬程序輸出的模擬時間單位，實際的時間單位秒可通過公式進(jìn)行轉(zhuǎn)換)，模擬結(jié)束時的模擬結(jié)果如圖9所示。

圖9 KMC模擬形核與長大過程結(jié)果可視圖

在圖9中，所有空位聚集在一起，形成一個大的空位團(tuán)簇(空洞)。由于MD過程中由高能快中子打出的的空位分布過于集中，在KMC模擬過程中所有空位聚集到一起的所需時間較短。并且，單個空位躍遷經(jīng)過的路徑也較短，所以每個空位與周圍金屬原子交換位置的次數(shù)也就較少。這也就使得體系中雖有銅團(tuán)簇形成(聚在一起的白色小球)，但直徑很小，且數(shù)量較少。上述結(jié)果符合輻照損傷的演化過程，說明MD與KMC耦合模擬以及整個并行耦合系統(tǒng)模擬是可行的。

4 結(jié)束語

本文采用MD和KMC耦合模擬的方法，對材料核輻照損傷過程中缺陷的產(chǎn)生以及形核和長大階段進(jìn)行模擬。通過對MD方法和KMC方法進(jìn)行研究，確立了用MD和KMC方法耦合進(jìn)行輻照損傷模擬的方案，即用MD方法模擬材料在持續(xù)的中子輻照下產(chǎn)生缺陷的過程，用KMC方法模擬缺陷形核和長大的過程。由于MD和KMC模擬中的原子表示及可模擬的原子類型不同，本文開發(fā)了一個耦合程序來實現(xiàn)MD和KMC的耦合模擬，在此過程中提出并實現(xiàn)了SD算法來完成原子與格點的映射，并在最后通過實驗驗證了SD算法以及MD和KMC耦合模擬的正確性和有效性。