原子動(dòng)力學(xué)蒙特卡洛程序Open KMC在反應(yīng)堆壓力容器鋼缺陷損傷研究中的優(yōu)化與應(yīng)用*

尚子豪，商紅慧，王東杰，張?jiān)迫?，賀新福，陳澤華，王 棟，張廣婷

（1.太原理工大學(xué)大數(shù)據(jù)學(xué)院，山西 太原030024；

2.中國(guó)科學(xué)院計(jì)算技術(shù)研究所計(jì)算機(jī)體系結(jié)構(gòu)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京100190；

3.中國(guó)原子能科學(xué)研究院，北京102413）

1 引言

核能是一種可以大規(guī)模使用的工業(yè)能源，核能發(fā)電具有發(fā)電成本較低、發(fā)電小時(shí)數(shù)長(zhǎng)、電網(wǎng)接入穩(wěn)定以及零碳排放等優(yōu)勢(shì)。核反應(yīng)堆壓力容器RPV（Reactor Pressure Vessel）是防止放射性物質(zhì)泄漏的最主要屏障，并且是反應(yīng)堆中不可拆換的部件，RPV因此成為核電站最關(guān)鍵的部件之一。RPV材料的使用壽命直接決定了核電站的服役壽命，在RPV鋼中加入少量的銅能顯著提高鋼的耐腐蝕性能，同時(shí)增加鋼的強(qiáng)度和疲勞抗力等。但是，長(zhǎng)時(shí)間暴露在高能中子輻照?qǐng)鲋袝?huì)導(dǎo)致RPV鋼中的銅元素以富銅團(tuán)簇的形式析出，這些高數(shù)量密度的納米富銅團(tuán)簇是引起RPV鋼產(chǎn)生脆化的主要原因[1，2]。

正常狀態(tài)下，由于銅的自擴(kuò)散系數(shù)較低，未受到輻照時(shí)銅原子基本不能通過(guò)擴(kuò)散聚集而形成團(tuán)簇，然而高能中子輻照產(chǎn)生的點(diǎn)缺陷會(huì)大大提高溶質(zhì)原子的擴(kuò)散聚集能力，從而產(chǎn)生原子團(tuán)簇。利用現(xiàn)代實(shí)驗(yàn)技術(shù)，如自擴(kuò)散示蹤法、電阻率法和正電子湮沒(méi)法等[3-5]對(duì)點(diǎn)缺陷相關(guān)性質(zhì)的實(shí)驗(yàn)測(cè)量，可以推斷出許多有關(guān)缺陷的信息。然而，由于分析方法和儀器分辨率的限制，用這些技術(shù)直接觀察和研究點(diǎn)缺陷團(tuán)簇的原子結(jié)構(gòu)和性質(zhì)很困難。利用計(jì)算機(jī)模擬的方法可以對(duì)各種點(diǎn)缺陷進(jìn)行深入的研究。

目前對(duì)RPV鋼的富銅團(tuán)簇的研究主要集中在輻照注量、輻照溫度以及元素含量對(duì)其性能的影響，而通過(guò)引入點(diǎn)缺陷的方式來(lái)觀察富銅團(tuán)簇析出過(guò)程的研究較少。動(dòng)力學(xué)蒙特卡洛KMC（Kinetic Monte Carlo）方法[6]可以有效地模擬富銅團(tuán)簇析出的微觀結(jié)構(gòu)演化。許多研究人員致力于開(kāi)發(fā)不同的KMC方法。然而，這類軟件多為串行版本，這使得對(duì)富銅團(tuán)簇析出過(guò)程的模擬只能在較小的體系中進(jìn)行，也因此只能引入單個(gè)或是少量的點(diǎn)缺陷，限制了模擬的精確程度。為此我們自主開(kāi)發(fā)了可以并行模擬上億原子體系的Open KMC軟件[7]。本文基于Pair勢(shì)[8]和嵌入原子勢(shì)EAM（Embedded-Atom Method）[9]，通過(guò)引入空位點(diǎn)缺陷的方式，模擬熱時(shí)效下Fe-Cu二元合金中團(tuán)簇的析出過(guò)程。并對(duì)程序中的空位點(diǎn)缺陷選取問(wèn)題進(jìn)行了計(jì)算優(yōu)化，通過(guò)多組測(cè)試驗(yàn)證了優(yōu)化算法的正確性和有效性，同時(shí)對(duì)優(yōu)化后的程序做了大規(guī)模并行模擬的性能分析。

2 程序優(yōu)化

2.1 基于空位躍遷的KMC方法

KMC模擬方法是一種通過(guò)構(gòu)造隨機(jī)過(guò)程來(lái)模擬體系長(zhǎng)時(shí)間演化的方法。經(jīng)典的原子動(dòng)力學(xué)蒙特卡洛AKMC（Atomic Kinetic Monte Carlo）方法[10]利用晶格點(diǎn)陣來(lái)映射原子與格點(diǎn)的位置關(guān)系，應(yīng)用于提供控制微觀系統(tǒng)演化的精確解。α-Fe的晶格為體心立方（bcc）結(jié)構(gòu)，空位躍遷事件被認(rèn)為是空位與其周圍第1近鄰原子的交換。圖1分別給出了bcc晶格中心原子的第1、第2和第3近鄰原子（晶格邊長(zhǎng)為2）。模擬通常包括4個(gè)步驟：（1）計(jì)算列表中空位的躍遷概率；（2）計(jì)算時(shí)間增量；（3）執(zhí)行空位躍遷，增加時(shí)間；（4）是否達(dá)到時(shí)間閾值，若達(dá)到則退出，否則循環(huán)執(zhí)行以上步驟。

Figure 1 Three kinds of nearest neighbors for a bcc lattice body center atom圖1 bcc晶格體心原子的三重近鄰

第1步根據(jù)當(dāng)前狀態(tài)計(jì)算空位可能發(fā)生的躍遷事件X的概率，這個(gè)概率實(shí)際上并非真正意義上的0～1的概率，而是一種頻率，其計(jì)算方式如下所示：

其中，Γ0代表嘗試頻率，嘗試頻率被設(shè)置為6×1012Hz，k是玻爾茲曼因子，T是絕對(duì)溫度，EXa表示事件發(fā)生需要的激活能?？紤]到計(jì)算時(shí)間和結(jié)果精度，我們采用環(huán)境相關(guān)模型[11]，得到了滿足式（2）所示的平衡規(guī)則的激活能EXa：

其中，EX0表示參考激活能，僅取決于與空位交換位置的原子的化學(xué)性質(zhì)。本文中使用的包括Fe：0.65 eV，Cu：0.56 e V。Ei和Ef分別表示空位躍遷前后系統(tǒng)的總能量值。

第2步采用駐留時(shí)間算法[12]計(jì)算時(shí)間增量，它與所有可能發(fā)生的事件頻率之和的倒數(shù)成正比，以便將事件的傳播與每個(gè)獨(dú)立事件相關(guān)聯(lián)。計(jì)算公式如式（3）所示：

其中隨機(jī)數(shù)r∈（0，1）用于選擇事件。

第3步是依據(jù)概率隨機(jī)確定空位的跳躍方向，發(fā)生跳躍并更新周圍原子的能量和躍遷概率。同時(shí)，在模擬時(shí)間上累積時(shí)間增量Δt。最后一步，重復(fù)上述步驟，直到模擬時(shí)間達(dá)到預(yù)設(shè)閾值。

在計(jì)算第1步中的能量時(shí)，原子間相互作用勢(shì)采用Pair勢(shì)能模型[8]或EAM勢(shì)能模型[9]，Pair勢(shì)由于其簡(jiǎn)潔的形式、較小的計(jì)算量和有效的計(jì)算結(jié)果，在原子層次的計(jì)算機(jī)模擬中得到了廣泛的應(yīng)用。Pair勢(shì)中原子間的相互作用最遠(yuǎn)只考慮到第2近鄰位置。組成系統(tǒng)的所有實(shí)體間的化學(xué)相互作用可以通過(guò)式（4）描述：

當(dāng)方程中a等于1或2時(shí)，分別對(duì)應(yīng)于第1或第2近鄰相互作用。Va表示空位，b為Fe-Fe鍵的數(shù)目，c為Va-Va鍵的數(shù)目，l為Fe-Va鍵的數(shù)目，m為Fe-Cu鍵的數(shù)目，n為Va-Cu鍵的數(shù)目，p為Cu-Cu鍵的數(shù)目。

EAM勢(shì)中原子間相互作用由Pair勢(shì)和多體嵌入勢(shì)共同貢獻(xiàn)，并且至多考慮到第3近鄰。其勢(shì)能計(jì)算方式如式（5）所示：其中，μ和v為原子表示，F(xiàn)（）為多體嵌入能函數(shù)，ρμ為μ處插值的總電子密度。相比于EAM勢(shì)能模型，Pair勢(shì)能模型的局限性在于缺乏多體嵌入勢(shì)能，不能考慮局部環(huán)境的影響。此外，Pair勢(shì)能模型需要通過(guò)實(shí)驗(yàn)或熱力學(xué)數(shù)據(jù)仔細(xì)調(diào)整。本文所用的Pair勢(shì)能模型參數(shù)取自文獻(xiàn)[13]，可以得到接近EAM模型的結(jié)果，也可以推廣到更復(fù)雜的多元合金研究當(dāng)中。

2.2 計(jì)算優(yōu)化

KMC算法中模型的最主要部分包括躍遷概率和系統(tǒng)能量的計(jì)算。如前面提到的，在模擬的第3步空位發(fā)生躍遷之后，需要對(duì)該空位周圍原子的能量和躍遷概率進(jìn)行更新。相當(dāng)于在下次空位躍遷事件發(fā)生之前，我們要對(duì)體系內(nèi)所有受到上一次躍遷影響的空位進(jìn)行更新。以往的KMC版本是將上次發(fā)生躍遷的原子附近理想半徑下的全部空位作為需要被更新的對(duì)象。事實(shí)上，這些空位當(dāng)中有很多因?yàn)榫嚯x較遠(yuǎn)，而并沒(méi)有受到上次躍遷事件的影響，導(dǎo)致模擬過(guò)程出現(xiàn)了很多冗余的模擬計(jì)算。因此，如何選取需要被更新的空位成為一個(gè)關(guān)鍵問(wèn)題。

被更新空位的選取與采用的勢(shì)能函數(shù)的計(jì)算方式有關(guān)。對(duì)于Pair勢(shì)能模型來(lái)說(shuō)，原子間的相互作用最遠(yuǎn)考慮到第2近鄰位置。某一空位是否受到上次躍遷影響取決于這個(gè)空位的第1近鄰原子是否與上次發(fā)生躍遷的原子存在相互作用。如果能找到發(fā)生躍遷的原子可以影響到的最遠(yuǎn)距離的空位位置，并將這個(gè)最遠(yuǎn)距離作為截?cái)喟霃剑看诬S遷事件結(jié)束后只更新截?cái)喟霃揭詢?nèi)的空位，就能避免大量的冗余計(jì)算。

截?cái)喟霃降耐茖?dǎo)過(guò)程在三維坐標(biāo)系中完成。假設(shè)發(fā)生躍遷交換的原子分別為i和j，并在坐標(biāo)系中固定i和j的位置坐標(biāo)。確定Pair勢(shì)能模型下最遠(yuǎn)處需要被更新的空位位置可以分為下面3個(gè)步驟：（1）找到j(luò)原子附近2層近鄰內(nèi)的全部原子并記錄坐標(biāo)；（2）分別計(jì)算這些原子與i原子的距離，找到其中一個(gè)與i原子距離最遠(yuǎn)的原子記為s；（3）記錄s原子的8個(gè)1nn近鄰原子坐標(biāo)，并找到其中與i原子距離值最大的對(duì)應(yīng)原子記為k。其中近鄰中的nn代表nearest neighbor，表示第幾近鄰的含義。通過(guò)以上3步得到最遠(yuǎn)處需要被更新的空位k位置。所有原子間位置關(guān)系如圖2所示，將立體圖展開(kāi)得到圖3中的三視圖，其中晶格邊長(zhǎng)為2。

Figure 2 Stereogram of Pair truncation radius derivation圖2 Pair截?cái)喟霃酵茖?dǎo)立體圖

Figure 3 Three views of Pair truncation radius derivation圖3 Pair截?cái)喟霃酵茖?dǎo)三視圖

Table 1 Nearest neighbor distance comparison table表1 近鄰距離對(duì)照表

通過(guò)相同的思路，可以確定EAM勢(shì)能模型下空位更新的截?cái)喟霃?。與Pair勢(shì)能模型的區(qū)別在于，EAM勢(shì)能模型的原子間相互作用考慮到第3近鄰位置，因此在尋找空位k步驟的第1步里要找到j(luò)原子附近3層近鄰內(nèi)的全部原子并記錄坐標(biāo)，接著執(zhí)行相同的第2步和第3步。計(jì)算出的EAM勢(shì)能模型下的截?cái)喟霃絛EAM=6。在不改變?cè)娱g相對(duì)距離的條件下能夠確定3種情況的原子間位置關(guān)系，這里給出其中一組所有原子的位置關(guān)系如圖4和圖5所示。

Figure 4 Stereogram of EAM truncation radius derivation圖4 EAM截?cái)喟霃酵茖?dǎo)立體圖

Figure 5 Three views of EAM truncation radius derivation圖5 EAM截?cái)喟霃酵茖?dǎo)三視圖

算法1為空位選取優(yōu)化算法的偽代碼實(shí)現(xiàn)。

算法1空位選取優(yōu)化算法

輸入：lattice_constant，vacant，input_i，events，d。

輸出：空位更新后的系統(tǒng)。

1.unitdis←2/lattice_constant；/*晶格常數(shù)轉(zhuǎn)換因子，用于獲取原子坐標(biāo)時(shí)將晶格常數(shù)轉(zhuǎn)換為2，從而與截?cái)喟霃降耐茖?dǎo)過(guò)程關(guān)聯(lián)起來(lái)*/

2.i←vacant[input_i]；//獲取發(fā)生躍遷事件的空位i

3.proball←sum（events[i].propensity）；/*空位i對(duì)應(yīng)的事件列表內(nèi)所有躍遷事件的概率和*/

4.ifproball＞0then

5.threshhold←proball*random（）；

6.j←threshhold；//選取發(fā)生躍遷交換的原子j

7.do_jump（i，j）；//發(fā)生躍遷事件

8.endif

9.d1，d2←0，0；

10.forkinvacantdo

11.round（i，j，k）；/*獲得躍遷原子i，j和查詢空位k的實(shí)際坐標(biāo)*/

12.round（i，j，k）*unitdis；/*獲得晶格常數(shù)轉(zhuǎn)換為2后各原子對(duì)應(yīng)的坐標(biāo)*/

13.d1←compute_linedis（i，k）；/*計(jì)算i與k的距離d1*/

14.d2←compute_linedis（j，k）；/*計(jì)算j與k的距離d2*/

15.ifd1≤dord2≤dthen

16.site_propensity（k）；/*更新空位k的躍遷概率*/

17. endif

18.endfor

在403規(guī)模的盒子中進(jìn)行模擬，分別驗(yàn)證Pair勢(shì)能模型和EAM勢(shì)能模型下優(yōu)化算法的正確性。首先用無(wú)優(yōu)化的程序，即不考慮截?cái)喟霃?，每次躍遷結(jié)束后都對(duì)全局空位進(jìn)行更新的方法對(duì)不同空位濃度進(jìn)行了7組測(cè)試，并得到了系統(tǒng)的初始能量和末態(tài)能量。接著用空位選取算法優(yōu)化后的程序在相同的輸入環(huán)境下同樣進(jìn)行7組測(cè)試，并得到了算法優(yōu)化后的系統(tǒng)始末能量，2種模型下的能量對(duì)比分別如表2和表3所示。

Table 2 Comparison of system energy before and after Pair potential model optimization表2 Pair勢(shì)能模型優(yōu)化前后系統(tǒng)能量對(duì)比

Table 3 Comparison of system energy before and after EAM potential model optimization表3 EAM勢(shì)能模型優(yōu)化前后系統(tǒng)能量對(duì)比

將2種模型無(wú)優(yōu)化時(shí)的末態(tài)能量畫成曲線的形式，優(yōu)化后的末態(tài)能量畫成點(diǎn)的形式得到圖6，可以看出，優(yōu)化后的點(diǎn)均落在無(wú)優(yōu)化時(shí)的曲線上面，即2種模型在初始能量相同的情況下，不同空位濃度下算法優(yōu)化后的末態(tài)能量均與無(wú)優(yōu)化時(shí)的末態(tài)能量保持一致，優(yōu)化后的程序并沒(méi)有影響模擬的演化進(jìn)程，保證了算法優(yōu)化后模擬結(jié)果的正確性。

Figure 6 Final energy comparison before and after optimization of two models圖6 2種模型各自優(yōu)化前后末態(tài)能量對(duì)比

接著為了測(cè)試優(yōu)化算法的效率提升，本文分別對(duì)Pair勢(shì)能模型和EAM勢(shì)能模型在無(wú)優(yōu)化和空位選取優(yōu)化算法下進(jìn)行了不同空位濃度的測(cè)試，并統(tǒng)計(jì)了模擬用時(shí)。另外，本文選用了動(dòng)力學(xué)模擬軟件Lakimoca[9]當(dāng)中的截?cái)喟霃綌?shù)值（d=10.5）進(jìn)行空位選取，同樣對(duì)2種勢(shì)能模型下不同空位濃度的體系進(jìn)行測(cè)試并得到相應(yīng)的模擬用時(shí)。最終得到了如圖7和圖8所示的效率提升對(duì)比圖。由圖7和圖8可知，隨著空位濃度的增加模擬用時(shí)越來(lái)越長(zhǎng)，EAM勢(shì)能模型下的無(wú)優(yōu)化模擬用時(shí)大幅度多于空位選取優(yōu)化的模擬用時(shí)，當(dāng)空位濃度為1.6%時(shí)，無(wú)優(yōu)化的模擬用時(shí)最高達(dá)到了空位選取優(yōu)化模擬用時(shí)的26.75倍。同時(shí)，Lakimoca空位選取方式下的模擬用時(shí)相比于無(wú)優(yōu)化時(shí)有很大的減少，但當(dāng)空位濃度達(dá)到1.6%時(shí)，仍有超出空位選取優(yōu)化一倍的時(shí)間開(kāi)銷。Pair勢(shì)能模型下3種方法的模擬用時(shí)與EAM勢(shì)能模型下的用時(shí)變化具有相同的增大趨勢(shì)，但由于Pair勢(shì)的計(jì)算方式本身比EAM勢(shì)更加簡(jiǎn)單，因此在相同空位濃度環(huán)境下比EAM勢(shì)用時(shí)更少。顯然，無(wú)論是Pair模型還是EAM模型，空位選取優(yōu)化算法下的模擬都具備更高的效率，尤其當(dāng)空位濃度升高時(shí)，能夠避免很大一部分由于空位更新帶來(lái)的冗余計(jì)算，使得模擬花費(fèi)的時(shí)間開(kāi)銷進(jìn)一步減小。

Figure 7 Comparison of Pair potential efficiency improvement圖7 Pair勢(shì)效率提升對(duì)比

Figure 8 Comparison of EAM potential efficiency improvement圖8 EAM勢(shì)效率提升對(duì)比

2.3 性能分析

為了得到優(yōu)化后的程序在大規(guī)模并行模擬當(dāng)中表現(xiàn)出的性能數(shù)據(jù)，本文從強(qiáng)可擴(kuò)展性和弱可擴(kuò)展性2方面對(duì)程序進(jìn)行了測(cè)試與分析。測(cè)試環(huán)境的相關(guān)配置參數(shù)如表4所示。

Table 4 Test environment related parameters表4 測(cè)試環(huán)境相關(guān)參數(shù)

強(qiáng)可擴(kuò)展性的測(cè)試是指固定模擬體系規(guī)模，只增加CPU核數(shù)，并觀察加速比和并行效率的變化。本文將模擬體系規(guī)模固定為10003，初始化模擬時(shí)間為t=1×10-4s，并設(shè)置用來(lái)控制自適應(yīng)同步算法的時(shí)間間隔參數(shù)nstop=10（nstop參數(shù)為每個(gè)子域設(shè)定一個(gè)可發(fā)生事件的平均值，結(jié)合所有子域發(fā)生事件的概率均值來(lái)計(jì)算子域的時(shí)間增量）。在溶質(zhì)Cu濃度為1.34%，空位濃度為0.08%，模擬溫度T=573 K的體系環(huán)境中，將CPU核數(shù)從1 000增加到5 000，分別對(duì)Pair勢(shì)和EAM勢(shì)進(jìn)行算例測(cè)試，得到各組算例的模擬用時(shí)，并以1 000核算例為基準(zhǔn)計(jì)算了其他各組算例的加速比和并行效率，測(cè)試結(jié)果分別如表5和表6所示。

Table 5 Test results of Pair potential strong scalability（Scale of simulation system：10003）表5 Pair勢(shì)強(qiáng)可擴(kuò)展性測(cè)試結(jié)果（模擬體系規(guī)模：10003）

Table 6 Test results of EAM potential strong scalability（Scale of simulation system：10003）表6 EAM勢(shì)強(qiáng)可擴(kuò)展性測(cè)試結(jié)果（模擬體系規(guī)模：10003）

2種勢(shì)能模型并行效率的變化過(guò)程如圖9所示。結(jié)合表5和表6可以看出，當(dāng)CPU核數(shù)從1 000增加到3 000時(shí)2種勢(shì)能模型的模擬用時(shí)下降可觀，呈現(xiàn)出不錯(cuò)的加速比和并行效率。當(dāng)CPU核數(shù)進(jìn)一步增加到5 000時(shí)加速比增長(zhǎng)速度繼續(xù)放緩，其中EAM勢(shì)在5 000核時(shí)的加速比已經(jīng)小于4 000核時(shí)的加速比，并且整體的并行效率也呈現(xiàn)下滑。這是因?yàn)镃PU核數(shù)增加使得進(jìn)程間的通信工作量變大，通信時(shí)間開(kāi)銷逐漸占有更高比重，使得并行效率呈現(xiàn)下滑的趨勢(shì)。5 000核的測(cè)試算例中，Pair勢(shì)的并行效率為52%，EAM勢(shì)的并行效率為49%，在大規(guī)模體系并行模擬中總體表現(xiàn)出良好的強(qiáng)可擴(kuò)展性。

Figure 9 Parallel efficiencies of strong scalability of two models圖9 2種模型的強(qiáng)可擴(kuò)展性的并行效率

弱可擴(kuò)展性的測(cè)試是指保持單個(gè)CPU核的工作負(fù)載不變，在增大模擬體系規(guī)模的同時(shí)以相同的比例增大CPU核數(shù)，并觀察并行效率的變化過(guò)程。弱可擴(kuò)展性的測(cè)試采用與強(qiáng)可擴(kuò)展性測(cè)試當(dāng)中相同溶質(zhì)濃度、空位濃度和模擬溫度的體系環(huán)境，將模擬體系規(guī)模從2003逐漸增大到10003，同時(shí)保持CPU核數(shù)以相同比例增長(zhǎng)，分別對(duì)Pair勢(shì)和EAM勢(shì)進(jìn)行測(cè)試，得到了各組算例的模擬用時(shí)，并以2003模擬規(guī)模的算例為基準(zhǔn)計(jì)算了其他各組算例的并行效率，測(cè)試結(jié)果分別如表7和表8所示。

Table 7 Test results of Pair potential weak scalability表7 Pair勢(shì)弱可擴(kuò)展性測(cè)試結(jié)果

弱可擴(kuò)展性測(cè)試中2種勢(shì)能模型的并行效率隨CPU核數(shù)增加的變化過(guò)程如圖10所示。結(jié)合表7和表8可以看出，隨著模擬體系規(guī)模和CPU核數(shù)的同比例增加，模擬用時(shí)并沒(méi)有保持相對(duì)恒定而是逐漸增加，通信時(shí)間增加帶來(lái)的負(fù)面影響變大。在5 000核測(cè)試當(dāng)中Pair勢(shì)的并行效率為0.60，EAM勢(shì)的并行效率為0.68，仍保持著較好的并行性能，因此程序在大規(guī)模體系并行模擬中總體具有良好的弱可擴(kuò)展性。

Table 8 Test results of EAM potential weak scalability表8 EAM勢(shì)弱可擴(kuò)展性測(cè)試結(jié)果

Figure 10 Parallel efficiencies of weak scalability of two models圖10 2種模型的弱可擴(kuò)展性的并行效率

3 數(shù)值模擬與分析

模擬采用自主開(kāi)發(fā)的Open KMC軟件，相比于傳統(tǒng)的原子動(dòng)力學(xué)蒙特卡洛（AKMC）串行方法的效率在大規(guī)模仿真中會(huì)受到模擬時(shí)間過(guò)長(zhǎng)和內(nèi)存溢出的限制，Open KMC在AKMC的基礎(chǔ)上將模擬體系分解成多個(gè)進(jìn)程，同時(shí)將每個(gè)進(jìn)程的工作負(fù)載劃分為更細(xì)粒度的扇區(qū)，以此避免進(jìn)程間的邊界沖突。利用眾核處理器，Open KMC可以實(shí)現(xiàn)和加速上億個(gè)原子的大規(guī)模KMC并行仿真。本文在文獻(xiàn)[7]的Open KMC程序基礎(chǔ)上實(shí)現(xiàn)了2.2節(jié)中的計(jì)算優(yōu)化，在不影響模擬結(jié)果的情況下使得模擬速度得到進(jìn)一步提高，并將優(yōu)化后的Open KMC程序應(yīng)用于接下來(lái)的大規(guī)模并行模擬當(dāng)中。另外，文獻(xiàn)[7]采用Open KMC基于Pair勢(shì)能模型進(jìn)行了大規(guī)模并行模擬，本文除了采用Pair勢(shì)能模型外，還將加入模擬結(jié)果更為可靠的EAM勢(shì)能模型做分析和比較。并行模擬所采用的測(cè)試環(huán)境的各項(xiàng)參數(shù)如表4所示。

3.1 并行模擬

首先，構(gòu)造規(guī)模為5003的bcc晶格Fe基體，原子總數(shù)為2.5×108。α-Fe的晶格常數(shù)為2.855?。溶劑原子是Fe，溶質(zhì)原子是占比為1.34%的Cu，體系中引入的空位點(diǎn)缺陷濃度為4×10-4%，模擬溫度T=573 K。在保證內(nèi)存不會(huì)溢出的前提下，為了能夠充分利用計(jì)算資源，本次模擬采用單節(jié)點(diǎn)64核進(jìn)行并行模擬。根據(jù)Vincent等人[14]的工作，本文將模擬時(shí)間進(jìn)行了重新標(biāo)度，從而獲得相應(yīng)的物理時(shí)間尺度。其時(shí)間標(biāo)度轉(zhuǎn)換公式如下所示：

在模擬剛開(kāi)始階段體系中孤立的Cu原子數(shù)目為3 006 580，體系引入的1 000個(gè)空位隨機(jī)分布在并行模擬的64個(gè)進(jìn)程當(dāng)中。圖11顯示了Pair勢(shì)能模型與EAM勢(shì)能模型下體系內(nèi)孤立銅原子的下降過(guò)程，模擬對(duì)應(yīng)的實(shí)際物理時(shí)間已經(jīng)通過(guò)式（6）和式（7）縮放為蒙特卡洛模擬時(shí)間。由于Pair勢(shì)能模型下的曲線在模擬初期下降很陡，圖11中右上角的小圖更好地顯示了2種模型在模擬的最早時(shí)間內(nèi)孤立銅原子的變化過(guò)程。2種模型體系中孤立Cu原子的數(shù)目均在一開(kāi)始迅速下降，接著下降速度逐漸放緩，最終保持著個(gè)位數(shù)的下降趨勢(shì)。顯然Pair勢(shì)相比EAM勢(shì)的模擬結(jié)果仍有差距，但考慮到2種模型下的結(jié)果呈現(xiàn)相同的趨勢(shì)，并且由于Pair勢(shì)能更小的計(jì)算量使得模擬用時(shí)大大減少，因此利用Pair勢(shì)能模型進(jìn)行數(shù)值模擬與分析在多數(shù)情況下也是較好的選擇。同時(shí)，根據(jù)本次實(shí)驗(yàn)的模擬結(jié)果也可以對(duì)Pair勢(shì)能模型參數(shù)做進(jìn)一步的調(diào)整，從而使得模擬結(jié)果更加能接近EAM勢(shì)能模型的結(jié)果。

Figure 11 Decline curve of the number of isolated Cu atoms圖11 孤立Cu原子數(shù)量下降曲線

利用EAM勢(shì)能模型下的模擬結(jié)果對(duì)實(shí)驗(yàn)現(xiàn)象進(jìn)行分析。經(jīng)統(tǒng)計(jì)在整個(gè)模擬過(guò)程的前0.1 s蒙特卡洛模擬時(shí)間里，體系內(nèi)已有超過(guò)50%的孤立Cu原子形成了團(tuán)簇。用時(shí)間標(biāo)度公式換算后相當(dāng)于Fe-Cu合金服役工作65年。利用粒子可視化軟件Ovito對(duì)0.1 s蒙特卡洛模擬時(shí)間下的形核結(jié)果進(jìn)行輸出。此時(shí)體系中已經(jīng)析出大量如圖12和圖13所示的富Cu團(tuán)簇和Cu-Va的復(fù)合體，其中淺色的圓球代表空位點(diǎn)缺陷，深色的圓球代表Cu原子。

Figure 12 Pure Cu clusters in the system圖12 體系內(nèi)的純Cu團(tuán)簇

Figure 13 Cu-Va complex in the system圖13 體系內(nèi)的Cu-Va復(fù)合體

為了驗(yàn)證體系中Cu和空位的富集現(xiàn)象，計(jì)算了原子間短程有序SRO（Short Range Order）參數(shù)值[15]，其中有序參數(shù)定義為：

Table 9 SRO value____表9 SRO值

將體系中每幀狀態(tài)下對(duì)應(yīng)的最大純Cu團(tuán)簇尺寸、最大Cu-Va復(fù)合體（空位數(shù)不少于10）尺寸統(tǒng)計(jì)輸出如圖14所示。因?yàn)轶w系中空位的數(shù)量遠(yuǎn)少于溶質(zhì)Cu原子的數(shù)量，在模擬的前期純Cu團(tuán)簇更易形成。但是，隨著模擬的進(jìn)行，空位逐漸聚集并與Cu原子組合形成Cu-Va復(fù)合體團(tuán)簇，復(fù)合體尺寸反超純Cu團(tuán)簇尺寸，成為了體系內(nèi)尺寸更大的團(tuán)簇。這是由于單空位向空位團(tuán)簇遷移可以使系統(tǒng)能量降低，因此形成更大更穩(wěn)定的團(tuán)簇；同時(shí)，Cu-Va間的結(jié)合能大于Cu-Cu間的結(jié)合能，組成的Cu-Va復(fù)合體可以進(jìn)一步降低體系的能量，這也是Cu-Va復(fù)合體團(tuán)簇在體系中可以穩(wěn)定存活的主要原因。

Figure 14 Comparison of maximum cluster size圖14 最大團(tuán)簇尺寸比較

在溶質(zhì)的擴(kuò)散聚集過(guò)程當(dāng)中綜合了空位形成能和遷移能的作用，以往的AKMC模擬大多在體系中只放入一個(gè)空位，并且體系的規(guī)模都很小，例如1×106個(gè)原子。而本文模擬結(jié)果表明，當(dāng)體系中空位數(shù)量不止1個(gè)時(shí)，體系中出現(xiàn)了很多較大的Cu-Va復(fù)合體團(tuán)簇，這說(shuō)明傳統(tǒng)的單空位模擬方案未必合理，我們不能忽略空位增加對(duì)團(tuán)簇形成產(chǎn)生的影響。本文模擬的體系中包含1 000個(gè)空位，為了研究空位數(shù)量對(duì)團(tuán)簇形成產(chǎn)生的影響，本文在與之前相同的體系環(huán)境中分別引入500，800，1 200，1 600和2 000個(gè)空位進(jìn)行了5組對(duì)比測(cè)試。首先是體系中孤立Cu原子的數(shù)目變化曲線對(duì)比，如圖15所示，可以看出，在沉淀析出過(guò)程趨于平穩(wěn)之前空位數(shù)量的增多加快了體系內(nèi)Cu原子的聚集過(guò)程。進(jìn)一步地，為了研究空位數(shù)量對(duì)團(tuán)簇的數(shù)量密度和尺寸產(chǎn)生的影響，提取了蒙特卡洛模擬時(shí)間0.1 s下各組測(cè)試中最大的團(tuán)簇，發(fā)現(xiàn)它們的尺寸均不超過(guò)40個(gè)原子。

Figure 15 Decline curve of isolated Cu atoms in different vacancy number systems圖15 不同空位數(shù)量體系下的孤立Cu原子下降曲線

本文假設(shè)體系中尺寸≥25個(gè)原子的團(tuán)簇為較大團(tuán)簇，將各組測(cè)試環(huán)境中全部團(tuán)簇（尺寸≥5個(gè)原子）的數(shù)量密度和尺寸≥25個(gè)原子的團(tuán)簇的數(shù)量密度統(tǒng)計(jì)輸出如圖16和圖17所示。圖16和圖17的縱坐標(biāo)表示單位體積內(nèi)的團(tuán)簇?cái)?shù)量?？梢钥吹剑S著空位數(shù)量的增加，體系內(nèi)總團(tuán)簇的數(shù)量密度相對(duì)穩(wěn)定且處于同一個(gè)數(shù)量級(jí)，沒(méi)有較大差異。與此同時(shí)空位越多，體系內(nèi)的較大團(tuán)簇?cái)?shù)量密度有明顯提高。此時(shí)體系內(nèi)的細(xì)小團(tuán)簇會(huì)出現(xiàn)分解，相對(duì)較大的團(tuán)簇通過(guò)吸收細(xì)小團(tuán)簇分解出的Cu原子來(lái)保持尺寸增長(zhǎng)。而空位增多能夠促進(jìn)這種團(tuán)簇粗化現(xiàn)象，進(jìn)而形成更多較大尺寸的沉淀。

3.2 實(shí)驗(yàn)比較

為了與真實(shí)的實(shí)驗(yàn)過(guò)程進(jìn)行比較，本文在403的體系下對(duì)663 K～773 K的溫度進(jìn)行了一系列熱時(shí)效模擬。其中溶劑為Fe，溶質(zhì)為占比1.34%的Cu，引入占比8×10-4%的空位濃度。并與Lê等人[5]的實(shí)驗(yàn)結(jié)果進(jìn)行了比較。本文對(duì)模擬的演化過(guò)程進(jìn)行了蒙特卡洛模擬時(shí)間10 s的跟蹤，利用式（6）和式（7）轉(zhuǎn)換得到不同溫度下實(shí)際演化的物理時(shí)間real_time=105s，并分析了這段時(shí)間的溶質(zhì)析出過(guò)程。105s對(duì)應(yīng)的蒙特卡洛模擬時(shí)間分別為：溫度為663 K時(shí)tMC_Pair=3.2×10-4s和tMC_EAM=9.53×10-3s；溫度為693 K時(shí)tMC_Pair=1.3×10-3s和tMC_EAM=3.3×10-2s；溫度為733 K時(shí)tMC_Pair=7.04×10-3s和tMC_EAM=0.15 s；溫度為773 K時(shí)tMC_Pair=3.2×10-2s和tMC_EAM=0.59 s。按照式（9）計(jì)算沉淀推進(jìn)系數(shù)，其中NX（0）是體系初始化時(shí)的孤立銅原子數(shù)量，NX（t）是時(shí)間t時(shí)體系內(nèi)孤立銅原子的數(shù)量，NX（∞）為無(wú)窮遠(yuǎn)處的體系內(nèi)孤立銅原子數(shù)量。

Figure 16 Evolution curve of total clusters number density under different vacancy number systems圖16 不同空位數(shù)量體系下的總團(tuán)簇?cái)?shù)量密度演化曲線

Figure 17 Evolution curve of number density of clusters with large size（≥25）圖17 較大尺寸（≥25）團(tuán)簇?cái)?shù)量密度演化曲線

ζ（t）=（NX（0）-NX（t））/（NX（0）-NX（∞））

（9）

圖18顯示了4種溫度下的Pair勢(shì)能模型和EAM勢(shì)能模型分別經(jīng)歷蒙特卡洛模擬時(shí)間10 s演化的體系內(nèi)孤立銅原子的下降過(guò)程。本文中NX（∞）選用蒙特卡洛模擬時(shí)間10 s時(shí)的孤立銅原子數(shù)量。將計(jì)算后的結(jié)果與電阻率測(cè)量得到的真實(shí)實(shí)驗(yàn)值進(jìn)行了比較，結(jié)果如圖19所示。

顯然，在熱時(shí)效的模擬過(guò)程中，Cu的析出是漸進(jìn)的。在這一部分工作中，文獻(xiàn)[7]中的模擬結(jié)果與實(shí)驗(yàn)值的對(duì)比顯示出了很好的吻合程度。而本文中的模擬結(jié)果與實(shí)驗(yàn)值曲線則有明顯的區(qū)別。這是由于之前的計(jì)算并沒(méi)有將NX（∞）取為理想的無(wú)窮遠(yuǎn)處的孤立銅原子數(shù)量，而是在模擬的體系尚未達(dá)到相對(duì)穩(wěn)定的狀態(tài)時(shí)取最后時(shí)刻的孤立銅原子數(shù)量作為NX（∞）并代入計(jì)算，導(dǎo)致計(jì)算出的推進(jìn)因子并不能實(shí)際反映模擬過(guò)程中的沉淀推進(jìn)過(guò)程。本文中的沉淀推進(jìn)過(guò)程與實(shí)驗(yàn)值的偏差來(lái)自簡(jiǎn)化模型的定量偏差，例如系統(tǒng)內(nèi)只有一種類型的缺陷，即空位點(diǎn)缺陷，以及原子間缺乏長(zhǎng)程相互作用。但是，模擬曲線至少可以定性地描述體系隨時(shí)間增加的演化過(guò)程，并且當(dāng)溫度升高時(shí)，EAM勢(shì)能模型的擬合結(jié)果與實(shí)驗(yàn)值更加接近，在溫度為773 K時(shí)，所得模擬結(jié)果與實(shí)驗(yàn)值基本吻合。簇粗化，形成尺寸較大的沉淀。最后，本文計(jì)算了4種不同溫度下模擬的沉淀推進(jìn)結(jié)果，并與實(shí)驗(yàn)值計(jì)算出的結(jié)果進(jìn)行了比較，在一定程度上本文的沉淀推進(jìn)過(guò)程與實(shí)驗(yàn)結(jié)果具有相同的趨勢(shì)，并且當(dāng)溫度升高時(shí)具備更好的擬合效果。

Figure 18 Evolution of isolated Cu atoms in Pair potential model and EAM potential model at different temperatures圖18 不同溫度下Pair勢(shì)能模型與EAM勢(shì)能模型的孤立Cu原子演化過(guò)程

Figure 19 Evolution of precipitation advance factor in Pair potential model and EAM potential model at different temperatures圖19 不同溫度下Pair勢(shì)能模型與EAM勢(shì)能模型沉淀推進(jìn)因子演化過(guò)程

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

本文基于Pair勢(shì)能模型和EAM勢(shì)能模型，針對(duì)Open KMC軟件中原子躍遷事件發(fā)生后的空位更新問(wèn)題進(jìn)行了算法優(yōu)化，在保證模擬結(jié)果與未優(yōu)化時(shí)一致的情況下，縮短了模擬用時(shí)，并且當(dāng)空位濃度升高時(shí)具有更好的效率提升。同時(shí)，對(duì)優(yōu)化后的程序做了大規(guī)模并行模擬的性能分析，其結(jié)果顯示程序具備良好的強(qiáng)可擴(kuò)展性和弱可擴(kuò)展性。進(jìn)而利用優(yōu)化后的Open KMC軟件并行模擬了熱時(shí)效下大規(guī)模體系Fe-Cu合金中引入多空位后的團(tuán)簇析出過(guò)程。模擬結(jié)果表明：（1）本文方法預(yù)測(cè)了體系中富Cu沉淀和Cu-Va復(fù)合體的析出，并通過(guò)計(jì)算SRO參數(shù)驗(yàn)證了原子聚集現(xiàn)象。（2）統(tǒng)計(jì)了體系內(nèi)最大團(tuán)簇尺寸的變化過(guò)程。發(fā)現(xiàn)盡管引入空位點(diǎn)缺陷的數(shù)量遠(yuǎn)小于溶質(zhì)Cu原子的數(shù)量，但體系內(nèi)形成的Cu-Va復(fù)合體相比于純Cu團(tuán)簇更易成為體系尺寸更大的團(tuán)簇。（3）在相同的體系環(huán)境下，增大引入合金中的空位數(shù)量，可以加快沉淀的析出過(guò)程。（4）空位數(shù)量的增加不會(huì)對(duì)體系內(nèi)總的團(tuán)簇?cái)?shù)量密度產(chǎn)生明顯的影響，但是能夠促進(jìn)團(tuán)