晶格格林函數在位錯多尺度模擬中的應用*

劉鳳華，李 微，劉晚果

(1 黔南民族師范學院物理與電子科學學院，貴州 都勻 558000；2 黔南民族師范學院化學化工學院，貴州 都勻 558000)

位錯是材料中一種重要的結構缺陷，材料的力學性能及其他特性均與位錯的幾何特征和運動行為有著密切聯系。理論上位錯芯結構的精確性是影響材料力學性能計算和預測的一個關鍵點，材料的電子結構及彈性行為均對位錯有著很強的敏感性。因此，有必要在原子尺度上對位錯進行電子結構的研究。

從宏觀上，經典連續介質理論已經可以成功地處理位錯的遠程應力場問題，在這個尺度上，我們需要的計算量較小；但是在位錯核心區原子畸變嚴重，電子效應明顯，因而需要采用精確的計算方法來描述這個尺度上的物理現象，我們一般采用的是密度泛函(DFT)計算方法。由于位錯長程應力場的存在，并且在垂直于位錯線的平面內對稱性已被破壞，因此DFT計算中的周期性邊界條件已不適用于含位錯體系的模擬。為了克服這一問題，有關位錯的DFT計算一般都采用位錯偶極子的模型結構[1-2]，構造出的位錯偶極子可以將其應力抵消掉這樣就可采用周期邊界條件，但是這種方法明顯忽略了位錯遠程應力場信息。而多尺度-跨層次方法摒棄了只依賴單一尺度的計算觀點，通過耦合量子力學區域和連續介質區域的物理量得到位錯的遠程應變場，進而實現了對位錯的直接模擬。

多尺度現象及相關理論方法是復雜物質體系研究中重要的科學問題[3-4]。多尺度方法的最初應用開始于20世紀70年代，Martin Karplus等[5]將量子化學和分子動力學相結合來研究酶的反應機制，通過計算模擬得到了在復雜系統中比傳統實驗更精確的預測結果，節省了大量的時間和經濟成本。而Martin Karplus等三人也因此獲得了2013年的諾貝爾化學獎。這也說明了在材料的研究進程中，多尺度方法的應用越來越收到人們的重視。尤其近幾年在材料基因組思想的影響下，采用多尺度方法模擬含缺陷體系已成為材料設計的一種主要理論研究手段。通過多尺度模擬，可以將材料的宏觀特性和其微觀結構有機結合起來，找出影響材料宏觀特性的微觀機制，更有效促進材料的研發進程。

1 多尺度模擬方法

多尺度方法的根本目標是發展有效方法，將不同尺度上的相關物理量耦合起來。根據耦合的強弱，多尺度算法可分為串行式(sequential)與并行式(concurrent)兩大基本類別[6]。在國內，王崇愚課題組最先開始了串行多尺度方法的研究，將第一性原理計算出的物理量，通過解析傳遞的方式得代入到經典彈性理論公式中，從而得出位錯的宏觀攀移速度特性，并且該結果和x射線實驗結果符合較好[7-8]，這也說明了該方法的可行性，但是這種方式只適用于尺度間耦合較弱的情況。并行多尺度一般的處理方法是將模擬體系劃分為不同尺度層次，每個層次上采用自己特定的處理方法，多個尺度上的計算同時進行，它適合于尺度間變量耦合比較強的情況。并行多尺度算法大致可以分為三大類，其中一類是基于能量耦合的方法[9-13]，通過構造全局的能量泛函得出相關物理量。MAAD(Macroscopic atomistic ab initio dynamics)方法是此類算法中一種典型算法，其基本思想是通過構造整個系統的哈密頓量Htot同時進行量子、原子以及宏觀尺度的計算，通過對總哈密頓的求導來求得系統的速度、加速度、軌跡及系統的構型[14]。但是這種方法很難應用到金屬體系中，這限制了該方法的應用。另一類算法是Nair等[15]發展的QM-CAD (Quantum Mechanics-coupled atomic discrete dislocation)以位移為耦合條件的多尺度方法，和其他多尺度方法中的區域劃分有所不同，在該方法中QM和過渡區的原子交界處多了一層界面原子，在計算中保持過渡區原子固定，只對QM區和界面原子做第一性原理計算，以得到的界面原子的位移作為邊界條件，由離散位錯本征方程確定出過渡區原子位移。第三類則是基于力參量耦合的方法，像Buehler等發展的量子力學作用力場與 Tersoff 經驗力場相結合的耦合模式[16]，還有Payne 發展的LOTF(Learn on the fly)方法[17-18]，該方法主要特點是分子動力學計算中需要的勢參數，由第一性原理計算結果進行擬合，將得到的分子動力學的分子構型選取量子力學部分做DFT計算，再擬合新的勢參數做為下一次計算的勢函數；由于過渡區原子上的力不參與勢函數的擬合，因此避免了在QM/MM(Quantum Mechanics/ Molecular Mechanics)邊界上力不匹配的問題。分子動力學和第一性原理迭代弛豫的最終結果是體系合力為零且滿足動量守恒。但是基于力耦合的多尺度方法普遍的一個問題是在弛豫過程中體系很難達到平衡狀態。

考慮到材料中位錯核心區域和遠程應變場對計算精度的要求不同，Sinclair 等[19]在晶格格林函數的基礎上發展起來柔性邊界條件的多尺度方法，這是一種典型的基于力耦合的多尺度方法。柔性邊界條件的引入就是為了消除原子模擬中在邊界上所產生的不匹配力現象。其處理方法是將含缺陷體系在空間尺度上分成三部分，用一個很小的原子區域來模擬復雜的缺陷結構，并將原子區和外部彈性區之間的力弛豫掉，而格林函數則充當了原子區和彈性區之間的“橋梁”。不同于其他位錯多尺度模擬方法像有限元法[13]和經驗勢場法[20]，該方法將缺陷的局域應變場和長程應力場耦合起來正在適用于越來越多的材料模擬中。

2 晶格格林函數

為了反映晶體中由于點缺陷的存在而引起的位移場變化情況，Madrdudin[21]于1963年提出了晶格格林函數(LGF)的概念。當初對晶格格林函數定義為：于Born-von Karman邊界條件基礎上，取聲子格林函數在頻率趨于零時的極限值，點缺陷所產生的晶格畸變位移依賴于該缺陷對晶體所施加的力和晶體的力常數矩陣。

在簡諧近似的理論下，晶體中原子在離開平衡位置位移u(l) 時體系能量為[22]：

(1)

式中Φ0為體系處在平衡狀態下的能量，φαβ(l,l’)為兩原子l和l’間3×3力常數矩陣，α,β=x,y,z：

(2)

(3)

(4)

在實空間晶格格林函數表示的是力常數矩陣的逆：

G(0)=(Φ)-1

(5)

其物理意義是當單位大小的力沿β方向作用在原子l’上時，周圍l原子在α方向上所產生的位移大小。當力fβ(l’)作用在原子l’上時，該力對周圍原子l所產生的位移可通過格林函數表示為：

(6)

所以，如果基體相的晶格格林函數可以確定，那么根據缺陷所產生的力，就可以很容易地求出缺陷所引起的晶格畸變位移場。因此，晶格格林函數的精確度決定著計算模擬的效率和精確性。早期的晶格格林函數的解析求解只限制在立方晶系金屬中[21,23]，由于含缺陷特別是位錯體系的格林函數計算起來特別復雜，一般的方式是用完整晶體的晶格格林函數來近似含缺陷體系的格林函數[24]，這無疑是增加了體系收斂過程中的迭代步數。近幾年含缺陷體系的LGF[25-26]，特別是Tan等[27]計算出了含位錯體系的晶格格林函數，大大降低了在弛豫過程中的迭代次數并提高了計算精度。

3 晶格格林函數在位錯模擬中的應用

晶格格林函數最初是應用于立方晶系的點缺陷的模擬計算中[21]，通過計算含點缺陷體系的格林函數來模擬點缺陷所產生的位移場。Tewary等[28]將彈性格林函數和晶格格林函數結合，對Si-Ge(001)量子阱中的點缺陷(間隙或空位)進行了多尺度模擬計算，從數值上得出了點缺陷附近的離散芯結構以及其所產生的應變場。

早期的晶格格林函數計算體系和范圍較小，為了更好地模擬晶體中的缺陷結構，S Rao等[29]在進一步發展了此算法，在位錯核心區采用原子間相互作用勢，實現了在較少原子數的晶胞中位錯的計算模擬。為了得到更到的精度，C Woodward[30]采用柔性邊界的格林函數方法成功模擬了體心立方金屬Mo和Ta中的位錯芯結構，分別對不同尺寸的晶胞進行模擬，發現采用柔性邊界條件后尺寸較小的晶胞仍能完整體現體系的位錯應變場。隨后，依據同樣的方法，C Woodward得出了具有L10結構的金屬間化合物TiAl合金中α/2<110>螺位錯的位錯芯附近的應變場[31]，從螺位錯的位錯芯結構推測出螺位錯易于發生交滑移，并且還發現該螺位錯可以分解出重要的刃位錯分量，由此計算出的晶格阻力是測得的單晶TiAl合金中屈服強度的兩倍。在面心立方金屬Al中，采用格林函數對Al中的螺位錯和刃位錯也都有相關計算[24]，結果顯示位錯芯寬度均和實驗符合很好，在此位錯芯基礎上，Curtin[32]進一步計算了溶質和位錯的相互作用，從而實現了理論上Al的屈服強度的預測，這是材料強度預測的一大進步。隨著算法的發展，格林函數的求解也由完整晶格格林函數向含缺陷結構的晶格格林函數轉變，從而使格林函數多尺度方法更普適于一般缺陷體系，如在含界面結構的體系中，Trinkle等[33]就計算出了在Ti中螺位錯和孿晶界的相互作用；以及在具有HCP結構的 Mg金屬中，根據柔性邊界條件得到了Mg中位錯芯結構以及和溶質原子的相互作用能[34]。通過將格林函數應用于不同金屬不同結構中的位錯線模擬[35-37]，均得到和實驗符合較好的結果，這也說明了該方法的可靠性。并且對于核外電子數較多的金屬鎳，通過格林函數方法[38]也可得到其平衡位錯芯結構，且位錯芯寬度和其他計算符合很好[35-36]。

4 晶格格林函數在模擬中的計算流程

對于含位錯體系，一般將整個模擬晶胞就分成三個區域：位錯芯區1，格林函數“橋接”區2和經典彈性理論的連續區3。一般在沿位錯線方向采用周期性邊界條件，離位錯線最近的區域取為位錯核心區，最外面為彈性連續區，取真空邊界條件。在一區位錯芯區，原子錯排較大，這時胡克定律不再適用，利用第一性原理計算對該區原子進行弛豫計算；三區為彈性連續區，該部分的處理可采用位錯遠程場的經典彈性理論；由于在一區和三區原子間，分別采用不同尺度的處理方法，這就會使得在兩區的交界處的原子上產生力不匹配的現象，這部分的原子屬于二區緩沖區，為了耦合一區和三區之間的力的不匹配，在這一區域上采用晶格格林函數的方法進行弛豫。

圖1 格林函數方法中模型的分區圖示[29]

如果采用完整晶格的晶格格林函數來處理含位錯體系模型，一般的方法是將位錯看成是一種微擾，由于完整晶格格林函數和含位錯體系的格林函數之間的差別，這就要求在使用完整格林函數弛豫過程中需不斷迭代計算過程，圖2為其迭代計算流程[39]。

圖2 格林函數和第一性原理迭代計算流程圖

5 結 語

多尺度算法在材料性能研究中的作用越來越重要。材料基因組計劃的宗旨就是在建立一個新的以計算模擬和理論預測優先、實驗驗證在后的新材料研發文化，從而取代現有的以經驗和實驗為主的材料研發的模式[40]。以晶格格林函數為“橋接”的多尺度力匹配方法，實現了對位錯芯結構的第一性原理計算，為理論預測提供了模型構型基礎，正在適用于越來越多的金屬材料，發揮著突出作用。