固體塑性實(shí)驗(yàn)室時(shí)間尺度的分子動(dòng)力學(xué)模擬概述

王云江

(1.中國科學(xué)院力學(xué)研究所 非線性力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京 100190；2.中國科學(xué)院大學(xué) 工程科學(xué)學(xué)院，北京 100049)

1 引 言

隨著超級(jí)計(jì)算機(jī)硬件的飛速發(fā)展,以及大規(guī)模并行技術(shù)的進(jìn)步,計(jì)算力學(xué)已經(jīng)深入到固體力學(xué)、流體力學(xué)以及工程力學(xué)研究的各方面,成為與實(shí)驗(yàn)和理論并駕齊驅(qū)的三大研究范式之一。如圖1所示的多尺度計(jì)算力學(xué)框架,計(jì)算力學(xué)觸及的范疇包括從電子和原子尺度的基礎(chǔ)研究到大型工程構(gòu)件的整體優(yōu)化,涉及量子力學(xué)到連續(xù)介質(zhì)力學(xué)的多層次跨越與耦合,以及與各種物理、化學(xué)和生物等外場(chǎng)的相互作用。

圖1 多尺度計(jì)算力學(xué)示意圖

多尺度計(jì)算力學(xué)包含多個(gè)層次和不同范疇的計(jì)算方法,如微觀原子尺度的密度泛函理論和分子動(dòng)力學(xué)、介觀尺度的相場(chǎng)、離散位錯(cuò)動(dòng)力學(xué)以及宏觀尺度的彈塑性有限元模擬。這些計(jì)算方法跨越從埃米(10-10m)到米、從阿秒(10-18s)到秒等10個(gè)數(shù)量級(jí)以上的空間和時(shí)間尺度,且理論基礎(chǔ)迥異。目前,不同尺度的力學(xué)問題通常需要考量計(jì)算的精度與效率,從而選擇合適的單尺度研究方法,而統(tǒng)一的跨尺度方案仍是一個(gè)棘手的問題。

其中,分子動(dòng)力學(xué)MD(Molecular dynamics)作為一種典型的原子尺度模擬技術(shù),半個(gè)多世紀(jì)以來,在物理、化學(xué)和生物等多個(gè)學(xué)科取得了重要的應(yīng)用,為原子層次觀察和解析宏觀現(xiàn)象背后所蘊(yùn)涵的微觀機(jī)理提供了一個(gè)重要的研究手段。分子動(dòng)力學(xué)模擬起始于Alder等[1,2]使用IBM 704計(jì)算機(jī)模擬了硬球模型的彈性碰撞問題。而第一個(gè)真實(shí)體系的MD模擬,為Gibson等[3]使用Born-Mayer排斥勢(shì)模型研究了固體銅的輻射損傷演化過程。之后,Rahman[4]使用Lennard-Jones經(jīng)驗(yàn)勢(shì)函數(shù)研究了液體氬的原子結(jié)構(gòu)特征和擴(kuò)散行為,并與實(shí)驗(yàn)數(shù)據(jù)對(duì)照。這些先驅(qū)工作,拉開了分子動(dòng)力學(xué)計(jì)算的序幕。

如圖2所示,從1960年開始,分子動(dòng)力學(xué)的研究論文從每年的數(shù)篇發(fā)展到2020年大約5萬篇之巨,已經(jīng)成為一個(gè)重要的研究領(lǐng)域,可見其為科學(xué)界帶來的影響。為此,對(duì)分子動(dòng)力學(xué)的發(fā)展作出重要貢獻(xiàn)的三位科學(xué)家Karplus，Levitt和Warshel在2013年獲得了諾貝爾化學(xué)獎(jiǎng),其頒獎(jiǎng)詞為“……為發(fā)展復(fù)雜化學(xué)系統(tǒng)的多尺度方法做出了卓越貢獻(xiàn)”。該結(jié)論也同樣適用于力學(xué)領(lǐng)域,尤其是固體力學(xué)與物理學(xué)和材料學(xué)相結(jié)合,使分子動(dòng)力學(xué)為理解固體宏觀塑性的微觀機(jī)理,以及液體和氣體的動(dòng)力學(xué)行為提供了原子尺度的微觀視角和模擬工具。

圖2 以molecular dynamics或atomistic simulation為關(guān)鍵詞在Web of Science數(shù)據(jù)庫中檢索到的論文隨年代的變化

分子動(dòng)力學(xué)在理性理解固體力學(xué)行為隱藏的微觀塑性機(jī)制方面具有絕對(duì)優(yōu)勢(shì),其中的一個(gè)杰出代表為李曉雁等[5]在揭示納米孿晶銅強(qiáng)化-軟化尺寸效應(yīng)方面的工作。如圖3(a)所示,基于大規(guī)模分子動(dòng)力學(xué)模擬,重現(xiàn)了實(shí)驗(yàn)中發(fā)現(xiàn)的隨孿晶間距的減小,納米孿晶銅屈服應(yīng)力先增大后減小的趨勢(shì),即發(fā)現(xiàn)從Hall-Petch效應(yīng)到逆Hall-Petch效應(yīng)的轉(zhuǎn)變。圖3(b,c)的原子圖像解釋了該轉(zhuǎn)變的微觀機(jī)制為位錯(cuò)塑性行為的轉(zhuǎn)變，即發(fā)生了與孿晶面具有一定傾角的位錯(cuò)形核(強(qiáng)化)到平行于孿晶面位錯(cuò)形核(弱化)的轉(zhuǎn)變。從而基于分子動(dòng)力學(xué)模擬提出了新的納米金屬的弱化機(jī)制,破解了逆Hall-Petch效應(yīng)的難題,體現(xiàn)了分子動(dòng)力學(xué)的魅力與重要性。

圖3 分子動(dòng)力學(xué)揭示納米孿晶銅中Hall-Petch到逆Hall-Petch效應(yīng)轉(zhuǎn)變的位錯(cuò)機(jī)制[5]

2 分子動(dòng)力學(xué)的時(shí)間尺度局限性

分子動(dòng)力學(xué)的思想為牛頓運(yùn)動(dòng)方程的數(shù)值積分。首先,根據(jù)經(jīng)驗(yàn)力場(chǎng)和系統(tǒng)中原子的位置確定每個(gè)原子的勢(shì)能；然后,將勢(shì)能對(duì)位置求導(dǎo)得出每個(gè)原子的受力；根據(jù)原子的初始速度、位置和加速度,基于牛頓運(yùn)動(dòng)方程,在比較短的時(shí)間(～10-15s)內(nèi),預(yù)測(cè)下一時(shí)刻原子的速度和位置,并且更新受力信息。如此往復(fù),直至達(dá)到模擬要求的運(yùn)行時(shí)間或其他收斂判據(jù)。由于數(shù)值積分的時(shí)間步長過短,通常百萬步的分子動(dòng)力學(xué)模擬,只能達(dá)到納秒量級(jí)(10-9s)的時(shí)間尺度[6]。而該時(shí)間尺度與實(shí)驗(yàn)室時(shí)間尺度(100s)存在數(shù)個(gè)數(shù)量級(jí)的區(qū)別,限制了分子動(dòng)力學(xué)在更大時(shí)間范圍的應(yīng)用。如分子動(dòng)力學(xué)在模擬諸如擴(kuò)散[7-9]、位錯(cuò)攀移[10]和蠕變[11]等常見的金屬材料塑性變形方面存在著顯著缺陷。

時(shí)間尺度問題為分子動(dòng)力學(xué)的可信性帶來了困擾。一個(gè)典型的問題為固體塑性變形的應(yīng)變率效應(yīng)[12-15]。如目前最先進(jìn)的分子動(dòng)力學(xué)模擬應(yīng)變率可以低至106s-1的量級(jí)(相當(dāng)于微秒級(jí)模擬)[16-18]。但是,該應(yīng)變率與典型的準(zhǔn)靜態(tài)實(shí)驗(yàn)條件(10-4s-1)仍存在10個(gè)數(shù)量級(jí)的區(qū)別。如圖4所示,實(shí)際金屬玻璃材料的強(qiáng)度隨應(yīng)變率會(huì)發(fā)生顯著的變化。注意,為了公平比較,此處的強(qiáng)度已用各自材料的模量(楊氏模擬E或剪切模量G)約化。圖4的曲線為應(yīng)變率相關(guān)的非線性擬合。可見,與實(shí)驗(yàn)值相比[19,20],分子動(dòng)力學(xué)遠(yuǎn)高估了材料的強(qiáng)度。原因?yàn)椴牧纤苄缘臒峒せ钍录l(fā)生速率具有強(qiáng)烈的應(yīng)力依賴性,高應(yīng)變率弛豫時(shí)間短,微觀結(jié)構(gòu)來不及充分優(yōu)化,導(dǎo)致應(yīng)力無法完成松弛[21-23]。

圖4 金屬玻璃的強(qiáng)度(模量約化)隨應(yīng)變率的變化

圖5 三維Ashby變形圖譜

為了在實(shí)驗(yàn)室時(shí)間尺度模擬固體材料的塑性變形行為,相應(yīng)的加速分子動(dòng)力學(xué)技術(shù)應(yīng)運(yùn)而生[14,26-28]。如果可以將分子動(dòng)力學(xué)的模擬時(shí)間擴(kuò)展到秒量級(jí),那么模擬的結(jié)果就可以同步捕捉到實(shí)驗(yàn)背后的微觀物理圖像(本文假設(shè)勢(shì)函數(shù)準(zhǔn)確)。本文提到的加速方法并不是指簡(jiǎn)單增加數(shù)值積分時(shí)間步長或者大規(guī)模并行實(shí)現(xiàn)。實(shí)際上,并行技術(shù)對(duì)于解決空間尺度問題比較有效,因?yàn)榭梢愿鶕?jù)超級(jí)計(jì)算機(jī)構(gòu)架采用分區(qū)并行策略。但是,時(shí)間序列上依次發(fā)生的動(dòng)力學(xué)過程,通常具有串行排列的特征,簡(jiǎn)單并行無法直接實(shí)現(xiàn)。為此,本文將討論幾種基于物理的加速技術(shù),在保持分子動(dòng)力學(xué)原子級(jí)精度的同時(shí),實(shí)現(xiàn)模擬時(shí)間的物理跨越。

3 加速分子動(dòng)力學(xué)策略與算例

目前流行的加速分子動(dòng)力學(xué)方法主要包括間接法與直接法。間接法是通過搜索最小能量(勢(shì)能)路徑確定激活能,基于簡(jiǎn)諧過渡態(tài)理論預(yù)測(cè)事件發(fā)生的頻率或時(shí)間尺度。此外,也可以通過繁瑣的自由能面抽樣過程,如路徑積分或基于傘形抽樣(umbrella sampling)的平均力勢(shì)(potential of mean-force)等方法計(jì)算自由能面,更精確地預(yù)測(cè)在有限溫度條件下事件發(fā)生的時(shí)間尺度。另一類加速分子動(dòng)力學(xué)方法為直接法,如元?jiǎng)恿W(xué)法(metadynamics)及其衍生方法和擴(kuò)展系綜模擬方法(replica exchange/parallel tempering)等。本文將簡(jiǎn)要介紹加速分子動(dòng)力學(xué)相關(guān)的統(tǒng)計(jì)力學(xué)基礎(chǔ),并通過具體的算例詳細(xì)闡述其基本內(nèi)核與實(shí)踐。

3.1 最小能量路徑搜索

該方法基于經(jīng)典的簡(jiǎn)諧過渡態(tài)理論[29],如圖6所示。該理論描述某一個(gè)物理、化學(xué)反應(yīng)或塑性事件,對(duì)應(yīng)于從一個(gè)構(gòu)型穩(wěn)定的能量極小出發(fā),沿某一廣義自由度(或序參量)描述的最小勢(shì)能路徑,翻越能壘,發(fā)生躍遷,并弛豫到臨近能量極小的全過程。注意,該路徑從其他任何維度觀察均為能量最小值。所以,能量最高點(diǎn)對(duì)應(yīng)于鞍點(diǎn),其勢(shì)能與初態(tài)能量極小的差為激活能ΔE。那么,基于簡(jiǎn)諧近似的過渡態(tài)理論預(yù)測(cè)此躍遷發(fā)生的頻率為

圖6 過渡態(tài)理論和元?jiǎng)恿W(xué)算法

(1)

目前,有很多有效的方法可以搜索復(fù)雜的高維勢(shì)能面，最流行的方法為微動(dòng)彈性帶法NEB(nudged elastic band)[30,31]。NEB根據(jù)路徑的復(fù)雜程度,使用插值等技術(shù)在初態(tài)與末態(tài)間插入數(shù)量不等的構(gòu)型鏡像。并將鏡像之間用彈性作用固定,其垂直分切面的位置通過勢(shì)能曲面的梯度進(jìn)行優(yōu)化,直至達(dá)到非受限方向的力收斂條件。但是,NEB方法需要事先明確初態(tài)與末態(tài)構(gòu)型,使用受限。針對(duì)無序體系中原子環(huán)境復(fù)雜和大激活能分布的特征[32-35],激活-弛豫方法ART(Activation-relaxation technique)可以實(shí)現(xiàn)大規(guī)模的高效抽樣[36,37]。并且,ART不需要指定末態(tài),從而為處理復(fù)雜體系的動(dòng)力學(xué)問題提供了更多選項(xiàng)。

圖7以純金屬晶界位錯(cuò)塑性為例,說明NEB方法如何實(shí)現(xiàn)跨時(shí)間尺度模擬晶體塑性。該算例基于銅Σ5{110}{411}傾斜晶界,建立雙晶體模型,在垂直于晶界方向上施加小于0 K屈服強(qiáng)度的拉伸應(yīng)力,即σ<σath。圖7(b)的0 K應(yīng)力-應(yīng)變曲線顯示該模型的強(qiáng)度上限為σath=4.34 GPa。圖7(c)的最小能量路徑曲線顯示,在σ=3.2 GPa的應(yīng)力條件下,位錯(cuò)形核塑性事件需要克服的激活焓為ΔH(σ)=1.82 eV。根據(jù)式(1)預(yù)測(cè),如此大的激活焓,使得該塑性事件在室溫(300 K)下發(fā)生的頻率為10-20Hz,即該事件幾乎不會(huì)發(fā)生。圖7(c)展示了激活焓隨應(yīng)力增大而減小的趨勢(shì)。只有當(dāng)施加的應(yīng)力大于3.47 eV時(shí),其對(duì)應(yīng)的激活焓才可能小于0.71 eV。那么,300 K發(fā)生該塑性事件的時(shí)間尺度大概為秒的量級(jí)。但是,對(duì)于分子動(dòng)力學(xué)來說,秒量級(jí)的過程仍然屬于稀有時(shí)間,遠(yuǎn)超過了經(jīng)典分子動(dòng)力學(xué)可以達(dá)到的時(shí)間尺度上限(微秒)。而通過NEB方法搜索動(dòng)力學(xué)過程的最小能量路徑,確定相應(yīng)的激活能,可以通過過渡態(tài)理論間接預(yù)測(cè)事件對(duì)應(yīng)的時(shí)間尺度,并與實(shí)驗(yàn)數(shù)據(jù)直接比較,佐證跨時(shí)間尺度計(jì)算的準(zhǔn)確性。

OC4 DeepCwind半潛式浮式風(fēng)機(jī)平臺(tái)有三組雙浮筒結(jié)構(gòu)、中柱、橫撐和斜撐組成。如圖10所示，即為其中的一個(gè)雙浮筒結(jié)構(gòu)。該浮筒分為上浮筒和下浮筒兩部分，中間由隔板分割。上、下浮筒為半潛式浮式風(fēng)力機(jī)系統(tǒng)提供浮力，浮筒中間為壓艙水，用來維持系統(tǒng)的穩(wěn)定性。當(dāng)浮筒發(fā)生破損時(shí)，平臺(tái)會(huì)發(fā)生傾斜，影響風(fēng)機(jī)系統(tǒng)的穩(wěn)定運(yùn)行。

圖7 搜索勢(shì)能平面預(yù)測(cè)時(shí)間尺度

此外,動(dòng)力學(xué)蒙特卡洛KMC(kinetic Monte Carlo)算法也可以克服時(shí)間尺度限制,但其強(qiáng)烈依賴于對(duì)某動(dòng)力學(xué)過程的物理認(rèn)知,即激活能為最關(guān)鍵的輸入?yún)?shù)。KMC是一種專門模擬體系長時(shí)間構(gòu)型演化的一類時(shí)間相關(guān)蒙特卡洛算法。首先,確定一類微觀過程,如空位擴(kuò)散和位錯(cuò)運(yùn)動(dòng)等,使用NEB等方法計(jì)算這些過程的激活能。然后,利用偽隨機(jī)數(shù)在能量極小態(tài)產(chǎn)生擾動(dòng),根據(jù)激活能大小,利用式(1)確定動(dòng)力學(xué)過程的反應(yīng)速率。之后,根據(jù)反應(yīng)速率隨機(jī)選擇反應(yīng)路徑,更新構(gòu)型和勢(shì)壘,累計(jì)總時(shí)間。以上過程不斷重復(fù),直至設(shè)定的時(shí)間要求,從而實(shí)現(xiàn)模擬的時(shí)間尺度跨越。Mousseau[38-40]提出的動(dòng)力學(xué)激活-弛豫方法(kinetic ART)可以根據(jù)原子環(huán)境變化適時(shí)搜索激活能,實(shí)現(xiàn)無平移和旋轉(zhuǎn)對(duì)稱性的復(fù)雜體系動(dòng)力學(xué)蒙特卡洛模擬。但是,KMC并沒有實(shí)現(xiàn)真正的熱擾動(dòng),也忽略了激活熵等不確定因素的作用,原則上只是對(duì)蒙特卡洛在時(shí)間閾上的粗?；?其精度強(qiáng)烈相關(guān)于構(gòu)型空間的復(fù)雜程度。

3.2 元?jiǎng)恿W(xué)(Metadynamics)

另一類跨時(shí)間尺度算法為直接法。其中,Parrinello等[41,42]提出的Metadynamics方法,為這一領(lǐng)域的開創(chuàng)性研究工作。鑒于國內(nèi)還沒有該方法準(zhǔn)確的中文翻譯,本文權(quán)且稱之為元?jiǎng)恿W(xué)方法。元?jiǎng)恿W(xué)完全基于經(jīng)典統(tǒng)計(jì)力學(xué),可以在統(tǒng)一理論框架內(nèi)加速復(fù)雜原子體系中的稀有時(shí)間,適時(shí)評(píng)估自由能曲面,達(dá)到延長原子模擬時(shí)間的目的。簡(jiǎn)言之,元?jiǎng)恿W(xué)使用一系列高斯函數(shù)形式的驅(qū)動(dòng)勢(shì)累加填充原始的勢(shì)能面盆地。這些高斯函數(shù)可以基于構(gòu)型的原子坐標(biāo)在運(yùn)行中適時(shí)調(diào)整位置和高度,基于抽象的廣義反應(yīng)坐標(biāo),即集體變量CV(Collective variable)定義。當(dāng)勢(shì)能盆地填平時(shí),系統(tǒng)可以在整個(gè)構(gòu)型空間等概率運(yùn)動(dòng),從而實(shí)現(xiàn)構(gòu)型空間的躍遷。此外,文獻(xiàn)[43-46]提出的ABC(autonomous basin climbing)算法,廣義上也屬于此類元?jiǎng)恿W(xué)跨時(shí)間尺度策略。

元?jiǎng)恿W(xué)的具體抽樣過程如圖6所示。首先,需要基于所有原子的坐標(biāo)矢量x定義集體變量s=s(x)。當(dāng)初始構(gòu)型處于熱力學(xué)平衡態(tài)時(shí),系統(tǒng)構(gòu)型圍繞平衡位置作熱振動(dòng)。元?jiǎng)恿W(xué)過程以s為序參量,根據(jù)即時(shí)位置添加高度為W,寬度為σ的驅(qū)動(dòng)勢(shì)

(2)

式中τ為每次抽樣,即添加驅(qū)動(dòng)勢(shì)的MD時(shí)間間隔,如1000步MD運(yùn)行時(shí)間,一共進(jìn)行k次添加高斯操作,從而運(yùn)行了時(shí)長為t的經(jīng)典MD抽樣。當(dāng)原始的勢(shì)能面填平時(shí),系統(tǒng)將從一個(gè)能量極小躍遷到臨近的能量較小。該過程對(duì)應(yīng)的自由能面可以通過回憶添加的驅(qū)動(dòng)式重構(gòu),即

F(s)=-ΔV(s)+Const.

(3)

常數(shù)項(xiàng)代表驅(qū)動(dòng)勢(shì)填滿整個(gè)相空間后,可以在不同的s處等概率巡游的結(jié)果。元?jiǎng)恿W(xué)完成后,所經(jīng)歷的分子動(dòng)力學(xué)抽樣時(shí)間記錄為tsampling。那么,考慮勢(shì)能面修正后的系統(tǒng)歷經(jīng)的等效物理時(shí)間為

(4)

圖8以非晶態(tài)合金Cu50Zr50的物質(zhì)擴(kuò)散為例,闡述元?jiǎng)恿W(xué)的計(jì)算過程[8]。因?yàn)閿U(kuò)散為單個(gè)原子或原子集團(tuán)的位置遷移,本文選擇跳躍原子相對(duì)于原始位置的跳躍距離為CV。圖8(a)為300 K條件下20納秒分子動(dòng)力學(xué)的抽樣過程,縱軸為CV隨時(shí)間的變化?？梢?CV起初圍繞0位置附近漲落；大概在11納秒附近,CV發(fā)生了一次跳躍,對(duì)應(yīng)于一個(gè)擴(kuò)散事件的發(fā)生。該擴(kuò)散事件的形貌為鏈狀形式,如圖8(c)所示。在整個(gè)20納秒的抽樣時(shí)間中,一個(gè)發(fā)生了幾次擴(kuò)散事件。圖7(b)為重構(gòu)后的原始自由能面。如第一個(gè)擴(kuò)散事件對(duì)應(yīng)的激活自由能為0.84 eV,相當(dāng)于100 s發(fā)生一次,從而元?jiǎng)恿W(xué)實(shí)現(xiàn)了在實(shí)驗(yàn)室時(shí)間尺度模擬擴(kuò)散過程,同時(shí)保留了原子細(xì)節(jié)。圖8(c)顯示擴(kuò)散事件的時(shí)間隨溫度變化的趨勢(shì)。從300 K到700 K,擴(kuò)散事件跨越了10個(gè)數(shù)量級(jí)以上的時(shí)間范圍,大大超出了經(jīng)典分子動(dòng)力學(xué)的計(jì)算能力。此外,元?jiǎng)恿W(xué)還捕捉到非晶態(tài)物質(zhì)的動(dòng)力學(xué)不均勻性特征[47],即同一個(gè)溫度條件下不同原子環(huán)境的擴(kuò)散過程對(duì)應(yīng)于非常寬頻時(shí)間譜。

圖8 跨時(shí)間尺度模擬非晶態(tài)物質(zhì)擴(kuò)散[8]

元?jiǎng)恿W(xué)不僅可以處理原子擴(kuò)散這種體積比較小的塑性機(jī)制,還可以處理涉及多原子集體運(yùn)動(dòng)的塑性形式。此處,本文以非晶態(tài)固體(Cu50Zr50)的剪切轉(zhuǎn)變塑性(shear transformation)[48]或β弛豫為例說明[49],具體結(jié)果如圖9所示。在該算例中,選取一團(tuán)近鄰原子,排除集團(tuán)的剛體轉(zhuǎn)動(dòng)和平動(dòng)后,原子相對(duì)于起始位置的均方位移RMSD作為CV。圖9(a)顯示的CV跳躍表示原子集團(tuán)發(fā)生了整體的重新排列。圖9(a)的位置信息顯示,該塑性事件為涉及4個(gè)近鄰原子的閉合環(huán)狀運(yùn)動(dòng)。由于可以在不同的溫度下深度抽樣該動(dòng)力學(xué)過程,故可得出各個(gè)溫度下的平均激活自由能與事件發(fā)生頻率,如圖9(b)所示?？梢?元?jiǎng)恿W(xué)模擬的結(jié)果趨勢(shì)與動(dòng)態(tài)力學(xué)分析儀的測(cè)量一致。并且,由于激活自由能ΔF=ΔE-TΔS,元?jiǎng)恿W(xué)還可以評(píng)估熵效應(yīng)(ΔS),這在NEB方法中是忽略的因素[25,50-52]。如對(duì)于非晶態(tài)物質(zhì)的局部激活,圖9(b)的高溫極限頻率推演出ΔS=14.8kB,對(duì)事件的發(fā)生頻率可以貢獻(xiàn)多個(gè)數(shù)量級(jí)的前因子。如此顯著的熵效應(yīng),揭示了激活路徑的復(fù)雜性[49]。

圖9 非晶剪切轉(zhuǎn)變塑性的元?jiǎng)恿W(xué)模擬[49]

3.3 自適應(yīng)加速分子動(dòng)力學(xué)(Adaptive -boost molecular dynamics)

雖然元?jiǎng)恿W(xué)實(shí)現(xiàn)了實(shí)驗(yàn)室時(shí)間尺度的分子動(dòng)力學(xué)模擬,但是其高斯形式的填充驅(qū)動(dòng)勢(shì)沒有明確的物理含義,而且高斯函數(shù)寬度和高度等參數(shù)的選擇具有一定的經(jīng)驗(yàn)性。所以,演化出各種元?jiǎng)恿W(xué)的變體。其中,以適宜緩和元?jiǎng)恿W(xué)(well-tempered metadynamics)的應(yīng)用最為廣泛[42]。其思想是將式(2)所示的填充高斯函數(shù)改造為隨時(shí)間衰減的形式,這樣可以避免勢(shì)能填充過飽和,限制構(gòu)型在物理相關(guān)的相空間抽樣。實(shí)際上,圖8和圖9所示的模擬均采用了適宜緩和元?jiǎng)恿W(xué)的形式。

除了對(duì)填充函數(shù)數(shù)學(xué)形式的改造,Ishii等[9]還從統(tǒng)計(jì)力學(xué)的角度出發(fā),發(fā)展了自適應(yīng)加速分子動(dòng)力學(xué)(Adaptive -boost molecular dynamics),提出了具有明確物理意義的填充函數(shù)。該方法從初始的局部平衡態(tài)出發(fā),首先進(jìn)行經(jīng)典分子動(dòng)力學(xué)模擬(如10000步MD),然后從哈密頓量H隨s的漲落評(píng)估感興趣CV的概率分布函數(shù)

(5)

(6)

此種形式的驅(qū)動(dòng)勢(shì)已經(jīng)考慮了原始自由能面的拓?fù)浣Y(jié)構(gòu)。所以,如果勢(shì)能盆地比較深,那么該填充驅(qū)動(dòng)勢(shì)就會(huì)比較大,從而使效率大大提高。另一方面,如果抽樣的勢(shì)能盆地比較淺,那么填充函數(shù)就會(huì)比較小,從而可以避免過填充,起到與適宜緩和元?jiǎng)恿W(xué)相似的作用。實(shí)踐證明,對(duì)于位錯(cuò)或擴(kuò)散塑性,個(gè)位數(shù)或幾十步的驅(qū)動(dòng)勢(shì)填充后,就可以觀測(cè)到構(gòu)型空間的實(shí)質(zhì)變化,從而使MD具備模擬稀有時(shí)間的功能。

圖10為使用自適應(yīng)加速方法,計(jì)算從銅晶體雙晶Σ9〈110〉{221}晶界上形核位錯(cuò)的算例[52]。圖10(a)顯示從10-4s-1到10-10s-1大應(yīng)變率范圍內(nèi)位錯(cuò)形核臨界應(yīng)力的變化趨勢(shì)。在不同的溫度下,均存在臨界應(yīng)力隨應(yīng)變率變化的不連續(xù)現(xiàn)象,對(duì)應(yīng)于形核機(jī)制的轉(zhuǎn)變。虛線對(duì)應(yīng)于經(jīng)典MD加載條件下的應(yīng)力響應(yīng),其微觀機(jī)制為多位錯(cuò)同時(shí)從晶界發(fā)射,而晶界并沒有發(fā)生重構(gòu),如圖10(b)所示?？墒?當(dāng)應(yīng)變率為典型的準(zhǔn)靜態(tài)實(shí)驗(yàn)條件10-3/s時(shí),微觀機(jī)制完全不同，單根位錯(cuò)通過晶界類擴(kuò)散重排獨(dú)立發(fā)射。這對(duì)經(jīng)典分子動(dòng)力學(xué)圖像提出了巨大挑戰(zhàn)。當(dāng)使用分子動(dòng)力學(xué)解析固體塑性機(jī)制時(shí),一定要留意時(shí)間尺度效應(yīng),不然可能誤導(dǎo)對(duì)微觀機(jī)制的認(rèn)知,在圖5所示的變形圖譜上錯(cuò)誤定位,極大誤導(dǎo)理性力學(xué)分析。

圖10 晶界位錯(cuò)形核機(jī)制的時(shí)間尺度效應(yīng)[52]

3.4 擴(kuò)展系綜分子動(dòng)力學(xué)

元?jiǎng)恿W(xué)和自適應(yīng)加速分子動(dòng)力學(xué)本質(zhì)上是借助統(tǒng)計(jì)力學(xué)的基本原理,在時(shí)間閾上的粗粒化操作,強(qiáng)烈依賴于過渡態(tài)理論的適用性和稀有事件假設(shè)。與之平行的最后一類方法為擴(kuò)展系綜分子動(dòng)力學(xué)模擬,包括副本交換分子動(dòng)力學(xué)(Parallel tempering or Multi-replica/replica exchange)[53]和溫度加速分子動(dòng)力學(xué)(Temperature -accelerated)[54-58]等。這類方法利用聯(lián)合的多個(gè)子系綜,每個(gè)系綜的溫度或某個(gè)物理量不同,而其他物理?xiàng)l件等價(jià),無論哪個(gè)子系綜有事件發(fā)生,都可以記入總事件數(shù)量并累積時(shí)間,從而使串行的時(shí)間序列事件實(shí)現(xiàn)并行化,達(dá)到直接加速分子動(dòng)力學(xué)模擬的目的。但該類方法強(qiáng)烈依賴于計(jì)算機(jī)并行規(guī)模,加速效率不高,如1000個(gè)處理器并行,原則上可以加速3個(gè)數(shù)量級(jí)的時(shí)間尺度,從而將典型的MD納秒尺度延長到微秒,但是與實(shí)驗(yàn)室時(shí)間尺度仍有不小差距。

4 挑戰(zhàn)與展望

分子動(dòng)力學(xué)的終極目標(biāo)為實(shí)現(xiàn)量子力學(xué)精度-大規(guī)模-實(shí)驗(yàn)室時(shí)間尺度的原子模擬。原則上,該目標(biāo)可以通過如量子霸權(quán)基于從頭算(ab initio)實(shí)現(xiàn)[59,60]。但是,目前最先進(jìn)的超級(jí)計(jì)算機(jī),通過大規(guī)模并行才可以實(shí)現(xiàn)十億量級(jí)、微秒尺度的直接分子動(dòng)力學(xué)模擬[16-18]。這樣的空間規(guī)模以及合適的邊界條件處理,基本上可以囊括材料中的各種局域或擴(kuò)展缺陷,包括位錯(cuò)、界面和析出相等[61],從而使模擬的空間尺度限制大大降低。所以,分子動(dòng)力學(xué)面臨的最艱巨的挑戰(zhàn)為勢(shì)函數(shù)的準(zhǔn)確性與時(shí)間尺度問題。

近些年,人工智能和機(jī)器學(xué)習(xí)結(jié)合材料基因工程思想和數(shù)據(jù)科學(xué),為材料科學(xué)提供了新的研究手段,甚至已經(jīng)并列為實(shí)驗(yàn)、理論和計(jì)算外的第四大研究范式。目前,分子動(dòng)力學(xué)可以模擬的體系受限于勢(shì)函數(shù)本身的準(zhǔn)確性,對(duì)于高溫合金和高熵合金等多組元體系的模擬受到很大的限制。為解決多元復(fù)雜體系勢(shì)函數(shù)構(gòu)建的普遍問題,Behler等[62]創(chuàng)新性地提出了基于高通量計(jì)算和機(jī)器學(xué)習(xí)構(gòu)建具有量子力學(xué)精度的神經(jīng)網(wǎng)絡(luò)勢(shì)函數(shù)的普適方案,從而有望擺脫多元?jiǎng)莺瘮?shù)的限制。鄂維南等[63,64]已經(jīng)基于此策略,實(shí)現(xiàn)了億量級(jí)和量子力學(xué)精度的直接分子動(dòng)力學(xué)模擬金屬材料力學(xué)行為,為此獲得了2020年度的戈登貝爾獎(jiǎng)。所以,基于機(jī)器學(xué)習(xí)與高通量數(shù)據(jù)庫構(gòu)建,建立任意元素組合的神經(jīng)網(wǎng)絡(luò)勢(shì)函數(shù),采用深度學(xué)習(xí)模型預(yù)測(cè)構(gòu)型演化,有望實(shí)現(xiàn)化學(xué)組分自由和量子精度的分子動(dòng)力學(xué)模擬,從而可以極大地拓展分子動(dòng)力學(xué)在力學(xué)、材料、物理、化學(xué)和生物等多學(xué)科中的應(yīng)用空間。

分子動(dòng)力學(xué)仍然受限于模擬的時(shí)間尺度問題,這可能是原子模擬最后的難點(diǎn)所在。雖然元?jiǎng)恿W(xué)、勢(shì)能面搜索與優(yōu)化方法在近20年取得了重要進(jìn)展,但是仍有諸多前沿問題需要解決。如集體變量(CV)的普適性問題,模擬加速的效率深度取決于選擇的CV能否代表最慢的物理過程。而對(duì)于復(fù)雜體系,選取合適的CV本身就構(gòu)成了不小的挑戰(zhàn)。所以,不依賴于CV的加速方法應(yīng)運(yùn)而生。Noé等[65]基于深度學(xué)習(xí)和統(tǒng)計(jì)力學(xué)提出了CV自由的玻爾茲曼生成器(Boltzmann generators)方法,可以跳過分子動(dòng)力學(xué)在平衡態(tài)構(gòu)型附近的長時(shí)間取樣,而直接從固有結(jié)構(gòu)出發(fā),生成一系列代表平衡態(tài)性質(zhì)的龐大等價(jià)構(gòu)型,通過深度學(xué)習(xí)驅(qū)使構(gòu)型能量符合正則(玻爾茲曼)分布?；趦?yōu)化的配分函數(shù)構(gòu)建自由能形式,并發(fā)現(xiàn)新的構(gòu)型,從而實(shí)現(xiàn)復(fù)雜結(jié)構(gòu)含時(shí)演化,而不需要任何先驗(yàn)的反應(yīng)坐標(biāo)知識(shí)。

最后,Invernizzi等[66]提出了增強(qiáng)抽樣的統(tǒng)一方案,將元?jiǎng)恿W(xué)與擴(kuò)展系綜模擬有機(jī)結(jié)合,通過副本交換等技術(shù),將不同系綜中不同集體變量相關(guān)的驅(qū)動(dòng)勢(shì)和不同的抽樣方法統(tǒng)一到整個(gè)擴(kuò)展系綜中,實(shí)現(xiàn)全局有效加速。這些進(jìn)展都為最終解決分子動(dòng)力學(xué)的時(shí)間尺度問題奠定了堅(jiān)實(shí)的基礎(chǔ)。