【文章內(nèi)容簡(jiǎn)介】 對(duì)物質(zhì)的組成、微觀結(jié)構(gòu)、性能以及服役性能進(jìn)行計(jì)算模擬與設(shè)計(jì)的學(xué)科。分子動(dòng)力學(xué)模擬[44]是指用計(jì)算機(jī)模擬原子或分子的運(yùn)動(dòng)過程, 以揭示材料內(nèi)在的本質(zhì)規(guī)律。近年來，分子動(dòng)力學(xué)模擬已成為納米領(lǐng)域研究最常用的計(jì)算方法。與蒙特卡洛方法、第一性原理，以及有限元一起成為計(jì)算機(jī)模擬技術(shù)的主流技術(shù)，被稱為計(jì)算機(jī)實(shí)驗(yàn)。計(jì)算機(jī)實(shí)驗(yàn)在當(dāng)今科學(xué)研究領(lǐng)域扮演著重要角色。在計(jì)算機(jī)實(shí)驗(yàn)中，通過建立理論模型，并借助計(jì)算機(jī)以及一定的算法和程序來計(jì)算。其中，蒙特卡洛方法是研究隨即性問題最有效的方法，可以研究材料介觀尺度問題。第一性原理方法可以無需任何實(shí)驗(yàn)數(shù)據(jù)，完全從材料組成的種類以及排列方式出發(fā)計(jì)算材料性能，其研究體系比分子動(dòng)力學(xué)小得多。有限元?jiǎng)t是研究宏觀尺寸的有效方法。這些計(jì)算機(jī)模擬方法被廣泛的應(yīng)用于科學(xué)研究，為理論研究和實(shí)驗(yàn)研究搭建橋梁，成為除理論研究和實(shí)驗(yàn)研究外探索自然界的另一種有效方法。模擬有時(shí)是一種理論，但有時(shí)又是一種實(shí)驗(yàn)。一方面，我們處理的是模型，而不是真實(shí)的物體，這表明模擬屬于一種理論方法；另一方面，通過計(jì)算機(jī)模擬證實(shí)一個(gè)模型的程序又非常類似于實(shí)驗(yàn)：我們運(yùn)行一個(gè)程序，然后分析結(jié)果，其方法非常類似于實(shí)驗(yàn)過程。計(jì)算機(jī)模擬、理論研究和實(shí)驗(yàn)研究三者相互聯(lián)系，相互補(bǔ)充，相互促進(jìn)，共同推動(dòng)著科學(xué)技術(shù)的發(fā)展。計(jì)算機(jī)模擬可以方便的觀察到許多實(shí)驗(yàn)無法觀察到的細(xì)節(jié)，甚至實(shí)驗(yàn)很難達(dá)到或是無法達(dá)到的條件，揭示材料的本質(zhì)規(guī)律。鑒于計(jì)算機(jī)模擬的優(yōu)點(diǎn)，本文使用分子動(dòng)力學(xué)方法來研究位錯(cuò)和Fe/Ni界面的相互作用機(jī)理，以及Cu/Ni界面的結(jié)合強(qiáng)度。分子動(dòng)力學(xué)模擬通過用經(jīng)典的牛頓運(yùn)動(dòng)方程來描述粒子的運(yùn)動(dòng)過程。 積分方法常用的積分方法有Velet算法、Leap Frog算法、VelocityVerlet算法和Beeman算法等。(1) Verlet算法在分子動(dòng)力學(xué)模擬中，最常用的時(shí)間積分方法為Verlet[45] 算法。： () ()式中r代表位移，v代表速度，a代表加速度，t代表時(shí)間。+δt時(shí)刻的位置。 ()。 () 或者 ()Verlet算法具有算法簡(jiǎn)單，存儲(chǔ)要求適度等優(yōu)點(diǎn)，但容易造成精度損失，且非自啟動(dòng)算法。（2） Leap Frog算法LeapFrog算法[46]，即蛙跳算法。其表達(dá)式如下： () 即得到 () () Leapfrog算法比Verlet 算法計(jì)算效率高、計(jì)算量小。但其缺點(diǎn)是位置和速度不是同步的。（3）VelocityVerlet算法VelocityVerlet算法[47]其優(yōu)點(diǎn)是計(jì)算量適中，且給出了顯式速度項(xiàng)，目前得到了廣泛的應(yīng)用。其表達(dá)式如下： () () 在VelocityVerlet算法中前一個(gè)時(shí)刻的速度、力和位置需要被儲(chǔ)存。（4）Beeman算法Beeman算法[48]也是在Verlet算法基礎(chǔ)上發(fā)展的另一種計(jì)算精度更高算法，但增加了其計(jì)算量。 原子間的相互作用勢(shì)在分子動(dòng)力學(xué)模擬中，勢(shì)函數(shù)是影響計(jì)算結(jié)果的關(guān)鍵因素之一，其決定著計(jì)算結(jié)果能否準(zhǔn)確地反映實(shí)際材料的性質(zhì)。原子間的相互作用力也可由勢(shì)函數(shù)求得，若勢(shì)函數(shù)存在誤差，則所計(jì)算的相互作用力也會(huì)受到影響。勢(shì)函數(shù)的準(zhǔn)確性決定著分子動(dòng)力學(xué)模擬的精度。(a) 對(duì)勢(shì)對(duì)勢(shì)通常有LennardJones勢(shì)、Morse勢(shì)和Johnson勢(shì)等勢(shì)模型。(1) LennardJones（LJ）勢(shì)[49]： ()式中ε為能量參數(shù)，rij為原子i和原子j之間的距離，σ為長(zhǎng)度參數(shù)。當(dāng)rij=σ時(shí)，勢(shì)能為零；當(dāng)rij?σ時(shí)，表現(xiàn)為引力勢(shì)，由式中第二項(xiàng)來表示，, 當(dāng)距離較遠(yuǎn)時(shí)，原子間的引力會(huì)越來越弱，趨近于零，通常為了計(jì)算方便，選取截?cái)喟霃絩c, 當(dāng)rij?rc時(shí)，取勢(shì)能為零，在LennardJones勢(shì)模型中，截?cái)喟霃絩c通常取21/6σ；當(dāng)rij?σ時(shí)，表現(xiàn)為斥力勢(shì)，由式中第一項(xiàng)來表示，當(dāng)原子間的距離越近，部分電子云發(fā)生重疊，原子間的斥力會(huì)越來越大。(2) Morse勢(shì)[50]： ()式中A、α、r0分別是結(jié)合能，調(diào)節(jié)參數(shù)和平衡間距，這三個(gè)值都是取經(jīng)驗(yàn)參數(shù)。金屬固體常用這種勢(shì)函數(shù)來描述，常用于金屬銅。(3) Johnson勢(shì)[51] ()αFe常用Johnson勢(shì)函數(shù)來描述。(b) 多體相互作用勢(shì)實(shí)際上，在多粒子系統(tǒng)中，就不能簡(jiǎn)單的應(yīng)用以上對(duì)勢(shì)模型，而必須考慮多體間的相互作用。常見的多體勢(shì)有嵌入原子法、StillingerWeber勢(shì)、FinnisSinclair勢(shì)等。(1) 嵌入原子勢(shì)（EAM）[52]在嵌入原子勢(shì)（Embedded Atom Method,簡(jiǎn)稱EMA）中，材料中的每一個(gè)原子被看作如同一個(gè)雜質(zhì)原子鑲嵌在包含所有其它原子的基體中。EMA勢(shì)可以處理各種各樣的金屬體系，其中包括斷裂、表面、雜質(zhì)，以及合金等體系。系統(tǒng)的總能量如下： ()EAM勢(shì)是將晶體的總勢(shì)能分成兩部分組成，上式中右端第一項(xiàng)фij代表對(duì)勢(shì)相互作用，rij為原子i和原子j之間的距離。第二項(xiàng)Fi即嵌入能，是電子密度ρi的函數(shù)，代表多體相互作用。式中1/2表示對(duì)勢(shì)為兩原子共有，式中ρi可以表達(dá)為： () 對(duì)于不同的金屬，EAM勢(shì)需通過擬合金屬的宏觀參數(shù)而得到。(2) StillingerWeber勢(shì)[53]StillingerWeber勢(shì)常用于描述半導(dǎo)體材料。 ()式中θjik 為ij鍵和ik鍵之間的夾角，g(r)為截?cái)喟霃皆诘谝唤徍偷诙徶g的衰變函數(shù)。截?cái)喟霃降倪x取是為計(jì)算方便，因?yàn)楫?dāng)距離非常遠(yuǎn)時(shí)，原子間的相互作用勢(shì)非常小，趨近于零，當(dāng)截?cái)喟霃饺〉迷酱?，結(jié)果越精確，但計(jì)算量也越大。當(dāng)cosθjik=1/3時(shí)，這種結(jié)構(gòu)類似金剛石四面體結(jié)構(gòu)，其結(jié)構(gòu)非常穩(wěn)定。(3) FinnisSinclair勢(shì)[54]以合金屬能帶緊束縛理論為基礎(chǔ)而發(fā)展的FinnisSinclair勢(shì)，給出了詳細(xì)的多體作用勢(shì)。其表達(dá)式如下： () ()式中，當(dāng)x0時(shí)，H(x)=1；當(dāng)x0時(shí)，H(x)=0。Ak、Rk、ak、rk為常數(shù)，且有R1R2，r1r2…r6。它們的值隨隨具體應(yīng)用的材料而定。 邊界條件選取合適的邊界條件，對(duì)分子動(dòng)力學(xué)模擬非常重要。但通常主要考慮以下兩點(diǎn)：首先，為了減小計(jì)算量，應(yīng)盡可能減小模擬的單元，但是為了排除可能的動(dòng)力學(xué)擾動(dòng)對(duì)結(jié)果造成影響，以及為了滿足統(tǒng)計(jì)學(xué)處理的可靠性要求，又要求模擬原胞應(yīng)該足夠大；其次，要從物理角度考慮應(yīng)變相容性、體積變化及環(huán)境的應(yīng)力平衡等實(shí)際耦合問題。當(dāng)模擬大塊體時(shí)，需使用三維周期性邊界條件；當(dāng)模擬薄膜材料時(shí)，需使用二維周期性邊界條件，如取X、Y方向?yàn)橹芷谛赃吔鐥l件，Z方向?yàn)樽杂蛇吔鐥l件；當(dāng)模擬納米線時(shí)，可選取一維周期性邊界條件，如取Z軸方向?yàn)橹芷谛赃吔鐥l件，X、Y方向?yàn)樽杂蛇吔鐥l件；有時(shí)又需使用混合邊界條件，即幾種邊界條件的組合。實(shí)際應(yīng)用時(shí)應(yīng)根據(jù)實(shí)際情況選擇合適的邊界條件。由于實(shí)際的宏觀物質(zhì)由無可計(jì)數(shù)的原子、分子組成，而計(jì)算機(jī)的計(jì)算能力又有限，所以為了進(jìn)行模擬，就必須采用周期性邊界條件。當(dāng)使用周期性邊界條件時(shí)，實(shí)際的粒子被假想在一個(gè)盒子里，這個(gè)盒子在空間三維方向被重復(fù)無窮多次，填充滿整個(gè)空間。即如果在盒子中位置r處有一粒子。 ()這里l, m, n為任意整數(shù)，a, b,c 為大小對(duì)應(yīng)盒子邊長(zhǎng)的矢量。所有這些鏡像粒子一起運(yùn)動(dòng)，但實(shí)際在計(jì)算程序中的只有一個(gè)粒子。這樣大大的減小了計(jì)算量，而且也不失精度。在基本分子動(dòng)力學(xué)元胞中，當(dāng)有原子離開分子動(dòng)力學(xué)元胞時(shí)就有對(duì)應(yīng)原子進(jìn)入其中。 分子動(dòng)力學(xué)元胞三維周期性邊界的二維圖示 溫度、壓力控制方法在分子動(dòng)力學(xué)模擬中，我們往往要考慮模擬系統(tǒng)是否與外界發(fā)生作用，如溫度、壓力等。對(duì)于不同的模擬對(duì)象，往往要對(duì)原子數(shù)N、體積V、壓強(qiáng)P、溫度T等量中的一些進(jìn)行控制，以實(shí)現(xiàn)不同的系綜。如在NVT系綜下需要控制原子數(shù)N、體積V、溫度T保持不變。而在NPT系綜下則需要對(duì)原子數(shù)N、壓強(qiáng)P、體積V進(jìn)行控制。 控溫方法在NVT、NPT系綜下，甚至NVE系綜下，我們期望模擬體系溫度保持恒定不變。而系統(tǒng)的溫度T與動(dòng)能K及速度有直接的關(guān)系，如下， ()kB、υi分別為Boltzmann常數(shù)和第i個(gè)原子的速度，N為原子個(gè)數(shù)， Nc為約束自由度數(shù)。溫度控制方法通常有速度標(biāo)度法、Berendsen熱浴法和Nos233。Hoover控溫方法等。(1) 速度標(biāo)度法[55]，故可通過對(duì)速度進(jìn)行標(biāo)度進(jìn)而達(dá)到控制溫度的目的。若速度乘以標(biāo)度因子λ后，t時(shí)刻的溫度的變化為 () ()Treq，λ分別為期望的參考溫度和速度標(biāo)度因子，給每一步乘上標(biāo)度因子λ來達(dá)到控制溫度的目的。(2) Berendsen熱浴法[56]Berendsen熱浴法是假想系統(tǒng)和一熱浴相接觸。對(duì)速度每一步都進(jìn)行標(biāo)度，標(biāo)度因子為 ()式中，?t，tT 分別為時(shí)間步長(zhǎng)和耦合參數(shù)，Tbath，T(t)分別為t時(shí)刻系統(tǒng)的期望溫度和瞬時(shí)溫度。（3）Nos233。Hoover控溫方法[5758]Nos233。Hoover控溫方法中，一恒定溫度的熱浴與體系相耦合。 控壓方法在NPT, NPH系綜下，我們需要模擬定壓情況下的分子動(dòng)力學(xué)問題，下面主要介紹Berendsen壓浴法和Andersen方法。（1）Berendsen壓浴法[59]Berendsen壓浴法是通過對(duì)體積乘以標(biāo)度因子Cp，來通過控制系統(tǒng)體積來實(shí)現(xiàn)對(duì)壓力的控制。如同讓一壓力為期望值的壓浴與系統(tǒng)接觸。 () 式中，P(t) 和Pbath分別為系統(tǒng)在t時(shí)刻的瞬時(shí)壓力和期望壓力；k和tP為耦合參數(shù)。(2) Andersen方法[60]，通過活塞調(diào)節(jié)體積變化來控制壓力。 活塞法調(diào)節(jié)壓力 分析方法 中心對(duì)稱參數(shù)（centrosymmetry parameter）為了可視化晶體缺陷，特別是位錯(cuò)的形核和發(fā)射過程，我們需要在晶體中可靠的確定這種缺陷結(jié)構(gòu)的位置?；诓煌臉?biāo)準(zhǔn)有許多不同的方法來研究晶體缺陷，如原子的勢(shì)能、配位數(shù)，以及中心對(duì)稱參數(shù)等。本文采用Kelchner等人[61]介紹的中心對(duì)稱參數(shù)方法。在固態(tài)物質(zhì)中，中心對(duì)稱參數(shù)對(duì)于區(qū)分原子為完整晶體，局部缺陷（位錯(cuò)或堆垛層錯(cuò)）或表面原子非常有用。其定義為： ()這里 α 為近鄰原子對(duì)數(shù)。對(duì)于 BCC結(jié)構(gòu)，α等于4；FCC結(jié)構(gòu)，α為6。Ri 和Ri+α為晶體中對(duì)稱的矢量。當(dāng)P值接近零時(shí)，對(duì)應(yīng)無缺陷的完整晶體。當(dāng)晶體中存在缺陷時(shí)，晶體的對(duì)稱性被打亂，此時(shí)P值就大于零，晶體的對(duì)稱性破壞越嚴(yán)重其P值就越大。在金屬材料中，不同的P值范圍，代表不同的晶體缺陷（部分位錯(cuò)、堆垛層錯(cuò)，或表面原子等），但對(duì)于不同的金屬材料會(huì)有些差異。如對(duì)于低溫下FCC結(jié)構(gòu)的Cu晶體中，其對(duì)應(yīng)原子為部分位錯(cuò)和堆垛層錯(cuò)，其對(duì)應(yīng)原子為表面原子。本文中運(yùn)用分子動(dòng)力學(xué)軟件LAMMPS計(jì)算原子的中心對(duì)稱參數(shù)，然后運(yùn)用Fortran程序?qū)τ?jì)算得出的數(shù)據(jù)進(jìn)行處理，提取出感興趣的量，然后通過VMD軟件來可視化晶體的缺陷結(jié)構(gòu)。 徑向分布函數(shù)徑向分布函數(shù)[62](Radial distribution function, RDF),也稱g(r),是描述與時(shí)間無關(guān)的粒子間關(guān)聯(lián)性的量度。 ()式中，P0為理想晶體的原子密度；n(r)為半徑在r到r+?r的球殼內(nèi)的原子數(shù)。由g(r)曲線，可對(duì)材料進(jìn)行結(jié)構(gòu)分析。如對(duì)于體心立方晶體和面心立方晶體徑向分布函數(shù)中波峰的第二近鄰即為晶格常數(shù)。，即為Fe納米晶體的晶格常數(shù)。 Fe晶體的徑向分布函數(shù) 模擬軟件與可視化本文采用LAMMPS模擬軟件運(yùn)算,并運(yùn)用Fortran，Matlab處理數(shù)據(jù)，且結(jié)合VMD軟件可視化。LAMMPS[64]是一款功能非常強(qiáng)大的分子模擬軟件。其優(yōu)點(diǎn)主要有兼容當(dāng)前大多數(shù)的勢(shì)能模型。 計(jì)算效率高，編程水平高?？梢阅M固體、氣體和液體等系統(tǒng)，可借助不同的邊界條件和力場(chǎng)模擬原子、分子、聚合物、金屬、顆粒等體系。且具有免費(fèi)開源的代碼，我們可以根據(jù)需要應(yīng)用或者修改其代碼。LAMMPS由美國Sandia國家實(shí)驗(yàn)室發(fā)布。其缺點(diǎn)是不可以直接可視化，要借助別的軟件來可視化，如VMD軟件、AtomEye軟件等。除了LAMMPS模擬軟件外，常用的分子模擬軟件還有Materials Explorer, NAMA, AM