李建江，魏 鵬，楊少峰，賀新福，胡長(zhǎng)軍

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

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

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

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

1 MD與KMC耦合模擬

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

1.1 MD方法與KMC方法原理

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

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

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

1.2 耦合模擬流程

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

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

1.3 耦合程序難點(diǎn)分析

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

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

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

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

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

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

2.1 Wigner-Seitz原胞缺陷分析法

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

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

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

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

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

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

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

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

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

圖6 一個(gè)格點(diǎn)(9號(hào))同時(shí)對(duì)應(yīng)兩個(gè)原子的情況

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

//定義頂點(diǎn)格點(diǎn)計(jì)數(shù)變量與體心格點(diǎn)計(jì)數(shù)變量

//定義原子分別到頂點(diǎn)格點(diǎn)的距離和體心格點(diǎn)的距離

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

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

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

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

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

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

計(jì)算與原子最近的體心格點(diǎn)(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;

將原子類型、格點(diǎn)坐標(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;

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

End if

End if

End for

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

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

計(jì)算與原子最近的體心格點(diǎn)(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;

將原子類型、格點(diǎn)坐標(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;

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

End if

End if

End for

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

3 實(shí)驗(yàn)結(jié)果與分析

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

3.1 SD算法驗(yàn)證

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

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

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

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

3.2 MD和KMC耦合模擬RPV鋼的輻照損傷實(shí)驗(yàn)

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

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

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

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

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

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

4 結(jié)束語(yǔ)

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