王云江
(1.中國科學院力學研究所 非線性力學國家重點實驗室,北京 100190;2.中國科學院大學 工程科學學院,北京 100049)
隨著超級計算機硬件的飛速發(fā)展,以及大規(guī)模并行技術的進步,計算力學已經(jīng)深入到固體力學、流體力學以及工程力學研究的各方面,成為與實驗和理論并駕齊驅(qū)的三大研究范式之一。如圖1所示的多尺度計算力學框架,計算力學觸及的范疇包括從電子和原子尺度的基礎研究到大型工程構件的整體優(yōu)化,涉及量子力學到連續(xù)介質(zhì)力學的多層次跨越與耦合,以及與各種物理、化學和生物等外場的相互作用。
圖1 多尺度計算力學示意圖
多尺度計算力學包含多個層次和不同范疇的計算方法,如微觀原子尺度的密度泛函理論和分子動力學、介觀尺度的相場、離散位錯動力學以及宏觀尺度的彈塑性有限元模擬。這些計算方法跨越從埃米(10-10m)到米、從阿秒(10-18s)到秒等10個數(shù)量級以上的空間和時間尺度,且理論基礎迥異。目前,不同尺度的力學問題通常需要考量計算的精度與效率,從而選擇合適的單尺度研究方法,而統(tǒng)一的跨尺度方案仍是一個棘手的問題。
其中,分子動力學MD(Molecular dynamics)作為一種典型的原子尺度模擬技術,半個多世紀以來,在物理、化學和生物等多個學科取得了重要的應用,為原子層次觀察和解析宏觀現(xiàn)象背后所蘊涵的微觀機理提供了一個重要的研究手段。分子動力學模擬起始于Alder等[1,2]使用IBM 704計算機模擬了硬球模型的彈性碰撞問題。而第一個真實體系的MD模擬,為Gibson等[3]使用Born-Mayer排斥勢模型研究了固體銅的輻射損傷演化過程。之后,Rahman[4]使用Lennard-Jones經(jīng)驗勢函數(shù)研究了液體氬的原子結構特征和擴散行為,并與實驗數(shù)據(jù)對照。這些先驅(qū)工作,拉開了分子動力學計算的序幕。
如圖2所示,從1960年開始,分子動力學的研究論文從每年的數(shù)篇發(fā)展到2020年大約5萬篇之巨,已經(jīng)成為一個重要的研究領域,可見其為科學界帶來的影響。為此,對分子動力學的發(fā)展作出重要貢獻的三位科學家Karplus,Levitt和Warshel在2013年獲得了諾貝爾化學獎,其頒獎詞為“……為發(fā)展復雜化學系統(tǒng)的多尺度方法做出了卓越貢獻”。該結論也同樣適用于力學領域,尤其是固體力學與物理學和材料學相結合,使分子動力學為理解固體宏觀塑性的微觀機理,以及液體和氣體的動力學行為提供了原子尺度的微觀視角和模擬工具。
圖2 以molecular dynamics或atomistic simulation為關鍵詞在Web of Science數(shù)據(jù)庫中檢索到的論文隨年代的變化
分子動力學在理性理解固體力學行為隱藏的微觀塑性機制方面具有絕對優(yōu)勢,其中的一個杰出代表為李曉雁等[5]在揭示納米孿晶銅強化-軟化尺寸效應方面的工作。如圖3(a)所示,基于大規(guī)模分子動力學模擬,重現(xiàn)了實驗中發(fā)現(xiàn)的隨孿晶間距的減小,納米孿晶銅屈服應力先增大后減小的趨勢,即發(fā)現(xiàn)從Hall-Petch效應到逆Hall-Petch效應的轉變。圖3(b,c)的原子圖像解釋了該轉變的微觀機制為位錯塑性行為的轉變,即發(fā)生了與孿晶面具有一定傾角的位錯形核(強化)到平行于孿晶面位錯形核(弱化)的轉變。從而基于分子動力學模擬提出了新的納米金屬的弱化機制,破解了逆Hall-Petch效應的難題,體現(xiàn)了分子動力學的魅力與重要性。
圖3 分子動力學揭示納米孿晶銅中Hall-Petch到逆Hall-Petch效應轉變的位錯機制[5]
分子動力學的思想為牛頓運動方程的數(shù)值積分。首先,根據(jù)經(jīng)驗力場和系統(tǒng)中原子的位置確定每個原子的勢能;然后,將勢能對位置求導得出每個原子的受力;根據(jù)原子的初始速度、位置和加速度,基于牛頓運動方程,在比較短的時間(~10-15s)內(nèi),預測下一時刻原子的速度和位置,并且更新受力信息。如此往復,直至達到模擬要求的運行時間或其他收斂判據(jù)。由于數(shù)值積分的時間步長過短,通常百萬步的分子動力學模擬,只能達到納秒量級(10-9s)的時間尺度[6]。而該時間尺度與實驗室時間尺度(100s)存在數(shù)個數(shù)量級的區(qū)別,限制了分子動力學在更大時間范圍的應用。如分子動力學在模擬諸如擴散[7-9]、位錯攀移[10]和蠕變[11]等常見的金屬材料塑性變形方面存在著顯著缺陷。
時間尺度問題為分子動力學的可信性帶來了困擾。一個典型的問題為固體塑性變形的應變率效應[12-15]。如目前最先進的分子動力學模擬應變率可以低至106s-1的量級(相當于微秒級模擬)[16-18]。但是,該應變率與典型的準靜態(tài)實驗條件(10-4s-1)仍存在10個數(shù)量級的區(qū)別。如圖4所示,實際金屬玻璃材料的強度隨應變率會發(fā)生顯著的變化。注意,為了公平比較,此處的強度已用各自材料的模量(楊氏模擬E或剪切模量G)約化。圖4的曲線為應變率相關的非線性擬合??梢?與實驗值相比[19,20],分子動力學遠高估了材料的強度。原因為材料塑性的熱激活事件發(fā)生速率具有強烈的應力依賴性,高應變率弛豫時間短,微觀結構來不及充分優(yōu)化,導致應力無法完成松弛[21-23]。
圖4 金屬玻璃的強度(模量約化)隨應變率的變化
圖5 三維Ashby變形圖譜
為了在實驗室時間尺度模擬固體材料的塑性變形行為,相應的加速分子動力學技術應運而生[14,26-28]。如果可以將分子動力學的模擬時間擴展到秒量級,那么模擬的結果就可以同步捕捉到實驗背后的微觀物理圖像(本文假設勢函數(shù)準確)。本文提到的加速方法并不是指簡單增加數(shù)值積分時間步長或者大規(guī)模并行實現(xiàn)。實際上,并行技術對于解決空間尺度問題比較有效,因為可以根據(jù)超級計算機構架采用分區(qū)并行策略。但是,時間序列上依次發(fā)生的動力學過程,通常具有串行排列的特征,簡單并行無法直接實現(xiàn)。為此,本文將討論幾種基于物理的加速技術,在保持分子動力學原子級精度的同時,實現(xiàn)模擬時間的物理跨越。
目前流行的加速分子動力學方法主要包括間接法與直接法。間接法是通過搜索最小能量(勢能)路徑確定激活能,基于簡諧過渡態(tài)理論預測事件發(fā)生的頻率或時間尺度。此外,也可以通過繁瑣的自由能面抽樣過程,如路徑積分或基于傘形抽樣(umbrella sampling)的平均力勢(potential of mean-force)等方法計算自由能面,更精確地預測在有限溫度條件下事件發(fā)生的時間尺度。另一類加速分子動力學方法為直接法,如元動力學法(metadynamics)及其衍生方法和擴展系綜模擬方法(replica exchange/parallel tempering)等。本文將簡要介紹加速分子動力學相關的統(tǒng)計力學基礎,并通過具體的算例詳細闡述其基本內(nèi)核與實踐。
該方法基于經(jīng)典的簡諧過渡態(tài)理論[29],如圖6所示。該理論描述某一個物理、化學反應或塑性事件,對應于從一個構型穩(wěn)定的能量極小出發(fā),沿某一廣義自由度(或序參量)描述的最小勢能路徑,翻越能壘,發(fā)生躍遷,并弛豫到臨近能量極小的全過程。注意,該路徑從其他任何維度觀察均為能量最小值。所以,能量最高點對應于鞍點,其勢能與初態(tài)能量極小的差為激活能ΔE。那么,基于簡諧近似的過渡態(tài)理論預測此躍遷發(fā)生的頻率為
圖6 過渡態(tài)理論和元動力學算法
(1)
目前,有很多有效的方法可以搜索復雜的高維勢能面,最流行的方法為微動彈性帶法NEB(nudged elastic band)[30,31]。NEB根據(jù)路徑的復雜程度,使用插值等技術在初態(tài)與末態(tài)間插入數(shù)量不等的構型鏡像。并將鏡像之間用彈性作用固定,其垂直分切面的位置通過勢能曲面的梯度進行優(yōu)化,直至達到非受限方向的力收斂條件。但是,NEB方法需要事先明確初態(tài)與末態(tài)構型,使用受限。針對無序體系中原子環(huán)境復雜和大激活能分布的特征[32-35],激活-弛豫方法ART(Activation-relaxation technique)可以實現(xiàn)大規(guī)模的高效抽樣[36,37]。并且,ART不需要指定末態(tài),從而為處理復雜體系的動力學問題提供了更多選項。
圖7以純金屬晶界位錯塑性為例,說明NEB方法如何實現(xiàn)跨時間尺度模擬晶體塑性。該算例基于銅Σ5{110}{411}傾斜晶界,建立雙晶體模型,在垂直于晶界方向上施加小于0 K屈服強度的拉伸應力,即σ<σath。圖7(b)的0 K應力-應變曲線顯示該模型的強度上限為σath=4.34 GPa。圖7(c)的最小能量路徑曲線顯示,在σ=3.2 GPa的應力條件下,位錯形核塑性事件需要克服的激活焓為ΔH(σ)=1.82 eV。根據(jù)式(1)預測,如此大的激活焓,使得該塑性事件在室溫(300 K)下發(fā)生的頻率為10-20Hz,即該事件幾乎不會發(fā)生。圖7(c)展示了激活焓隨應力增大而減小的趨勢。只有當施加的應力大于3.47 eV時,其對應的激活焓才可能小于0.71 eV。那么,300 K發(fā)生該塑性事件的時間尺度大概為秒的量級。但是,對于分子動力學來說,秒量級的過程仍然屬于稀有時間,遠超過了經(jīng)典分子動力學可以達到的時間尺度上限(微秒)。而通過NEB方法搜索動力學過程的最小能量路徑,確定相應的激活能,可以通過過渡態(tài)理論間接預測事件對應的時間尺度,并與實驗數(shù)據(jù)直接比較,佐證跨時間尺度計算的準確性。
OC4 DeepCwind半潛式浮式風機平臺有三組雙浮筒結構、中柱、橫撐和斜撐組成。如圖10所示,即為其中的一個雙浮筒結構。該浮筒分為上浮筒和下浮筒兩部分,中間由隔板分割。上、下浮筒為半潛式浮式風力機系統(tǒng)提供浮力,浮筒中間為壓艙水,用來維持系統(tǒng)的穩(wěn)定性。當浮筒發(fā)生破損時,平臺會發(fā)生傾斜,影響風機系統(tǒng)的穩(wěn)定運行。
圖7 搜索勢能平面預測時間尺度
此外,動力學蒙特卡洛KMC(kinetic Monte Carlo)算法也可以克服時間尺度限制,但其強烈依賴于對某動力學過程的物理認知,即激活能為最關鍵的輸入?yún)?shù)。KMC是一種專門模擬體系長時間構型演化的一類時間相關蒙特卡洛算法。首先,確定一類微觀過程,如空位擴散和位錯運動等,使用NEB等方法計算這些過程的激活能。然后,利用偽隨機數(shù)在能量極小態(tài)產(chǎn)生擾動,根據(jù)激活能大小,利用式(1)確定動力學過程的反應速率。之后,根據(jù)反應速率隨機選擇反應路徑,更新構型和勢壘,累計總時間。以上過程不斷重復,直至設定的時間要求,從而實現(xiàn)模擬的時間尺度跨越。Mousseau[38-40]提出的動力學激活-弛豫方法(kinetic ART)可以根據(jù)原子環(huán)境變化適時搜索激活能,實現(xiàn)無平移和旋轉對稱性的復雜體系動力學蒙特卡洛模擬。但是,KMC并沒有實現(xiàn)真正的熱擾動,也忽略了激活熵等不確定因素的作用,原則上只是對蒙特卡洛在時間閾上的粗粒化,其精度強烈相關于構型空間的復雜程度。
另一類跨時間尺度算法為直接法。其中,Parrinello等[41,42]提出的Metadynamics方法,為這一領域的開創(chuàng)性研究工作。鑒于國內(nèi)還沒有該方法準確的中文翻譯,本文權且稱之為元動力學方法。元動力學完全基于經(jīng)典統(tǒng)計力學,可以在統(tǒng)一理論框架內(nèi)加速復雜原子體系中的稀有時間,適時評估自由能曲面,達到延長原子模擬時間的目的。簡言之,元動力學使用一系列高斯函數(shù)形式的驅(qū)動勢累加填充原始的勢能面盆地。這些高斯函數(shù)可以基于構型的原子坐標在運行中適時調(diào)整位置和高度,基于抽象的廣義反應坐標,即集體變量CV(Collective variable)定義。當勢能盆地填平時,系統(tǒng)可以在整個構型空間等概率運動,從而實現(xiàn)構型空間的躍遷。此外,文獻[43-46]提出的ABC(autonomous basin climbing)算法,廣義上也屬于此類元動力學跨時間尺度策略。
元動力學的具體抽樣過程如圖6所示。首先,需要基于所有原子的坐標矢量x定義集體變量s=s(x)。當初始構型處于熱力學平衡態(tài)時,系統(tǒng)構型圍繞平衡位置作熱振動。元動力學過程以s為序參量,根據(jù)即時位置添加高度為W,寬度為σ的驅(qū)動勢
(2)
式中τ為每次抽樣,即添加驅(qū)動勢的MD時間間隔,如1000步MD運行時間,一共進行k次添加高斯操作,從而運行了時長為t的經(jīng)典MD抽樣。當原始的勢能面填平時,系統(tǒng)將從一個能量極小躍遷到臨近的能量較小。該過程對應的自由能面可以通過回憶添加的驅(qū)動式重構,即
F(s)=-ΔV(s)+Const.
(3)
常數(shù)項代表驅(qū)動勢填滿整個相空間后,可以在不同的s處等概率巡游的結果。元動力學完成后,所經(jīng)歷的分子動力學抽樣時間記錄為tsampling。那么,考慮勢能面修正后的系統(tǒng)歷經(jīng)的等效物理時間為
(4)
圖8以非晶態(tài)合金Cu50Zr50的物質(zhì)擴散為例,闡述元動力學的計算過程[8]。因為擴散為單個原子或原子集團的位置遷移,本文選擇跳躍原子相對于原始位置的跳躍距離為CV。圖8(a)為300 K條件下20納秒分子動力學的抽樣過程,縱軸為CV隨時間的變化。可見,CV起初圍繞0位置附近漲落;大概在11納秒附近,CV發(fā)生了一次跳躍,對應于一個擴散事件的發(fā)生。該擴散事件的形貌為鏈狀形式,如圖8(c)所示。在整個20納秒的抽樣時間中,一個發(fā)生了幾次擴散事件。圖7(b)為重構后的原始自由能面。如第一個擴散事件對應的激活自由能為0.84 eV,相當于100 s發(fā)生一次,從而元動力學實現(xiàn)了在實驗室時間尺度模擬擴散過程,同時保留了原子細節(jié)。圖8(c)顯示擴散事件的時間隨溫度變化的趨勢。從300 K到700 K,擴散事件跨越了10個數(shù)量級以上的時間范圍,大大超出了經(jīng)典分子動力學的計算能力。此外,元動力學還捕捉到非晶態(tài)物質(zhì)的動力學不均勻性特征[47],即同一個溫度條件下不同原子環(huán)境的擴散過程對應于非常寬頻時間譜。
圖8 跨時間尺度模擬非晶態(tài)物質(zhì)擴散[8]
元動力學不僅可以處理原子擴散這種體積比較小的塑性機制,還可以處理涉及多原子集體運動的塑性形式。此處,本文以非晶態(tài)固體(Cu50Zr50)的剪切轉變塑性(shear transformation)[48]或β弛豫為例說明[49],具體結果如圖9所示。在該算例中,選取一團近鄰原子,排除集團的剛體轉動和平動后,原子相對于起始位置的均方位移RMSD作為CV。圖9(a)顯示的CV跳躍表示原子集團發(fā)生了整體的重新排列。圖9(a)的位置信息顯示,該塑性事件為涉及4個近鄰原子的閉合環(huán)狀運動。由于可以在不同的溫度下深度抽樣該動力學過程,故可得出各個溫度下的平均激活自由能與事件發(fā)生頻率,如圖9(b)所示??梢?元動力學模擬的結果趨勢與動態(tài)力學分析儀的測量一致。并且,由于激活自由能ΔF=ΔE-TΔS,元動力學還可以評估熵效應(ΔS),這在NEB方法中是忽略的因素[25,50-52]。如對于非晶態(tài)物質(zhì)的局部激活,圖9(b)的高溫極限頻率推演出ΔS=14.8kB,對事件的發(fā)生頻率可以貢獻多個數(shù)量級的前因子。如此顯著的熵效應,揭示了激活路徑的復雜性[49]。
圖9 非晶剪切轉變塑性的元動力學模擬[49]
雖然元動力學實現(xiàn)了實驗室時間尺度的分子動力學模擬,但是其高斯形式的填充驅(qū)動勢沒有明確的物理含義,而且高斯函數(shù)寬度和高度等參數(shù)的選擇具有一定的經(jīng)驗性。所以,演化出各種元動力學的變體。其中,以適宜緩和元動力學(well-tempered metadynamics)的應用最為廣泛[42]。其思想是將式(2)所示的填充高斯函數(shù)改造為隨時間衰減的形式,這樣可以避免勢能填充過飽和,限制構型在物理相關的相空間抽樣。實際上,圖8和圖9所示的模擬均采用了適宜緩和元動力學的形式。
除了對填充函數(shù)數(shù)學形式的改造,Ishii等[9]還從統(tǒng)計力學的角度出發(fā),發(fā)展了自適應加速分子動力學(Adaptive -boost molecular dynamics),提出了具有明確物理意義的填充函數(shù)。該方法從初始的局部平衡態(tài)出發(fā),首先進行經(jīng)典分子動力學模擬(如10000步MD),然后從哈密頓量H隨s的漲落評估感興趣CV的概率分布函數(shù)
(5)
(6)
此種形式的驅(qū)動勢已經(jīng)考慮了原始自由能面的拓撲結構。所以,如果勢能盆地比較深,那么該填充驅(qū)動勢就會比較大,從而使效率大大提高。另一方面,如果抽樣的勢能盆地比較淺,那么填充函數(shù)就會比較小,從而可以避免過填充,起到與適宜緩和元動力學相似的作用。實踐證明,對于位錯或擴散塑性,個位數(shù)或幾十步的驅(qū)動勢填充后,就可以觀測到構型空間的實質(zhì)變化,從而使MD具備模擬稀有時間的功能。
圖10為使用自適應加速方法,計算從銅晶體雙晶Σ9〈110〉{221}晶界上形核位錯的算例[52]。圖10(a)顯示從10-4s-1到10-10s-1大應變率范圍內(nèi)位錯形核臨界應力的變化趨勢。在不同的溫度下,均存在臨界應力隨應變率變化的不連續(xù)現(xiàn)象,對應于形核機制的轉變。虛線對應于經(jīng)典MD加載條件下的應力響應,其微觀機制為多位錯同時從晶界發(fā)射,而晶界并沒有發(fā)生重構,如圖10(b)所示??墒?當應變率為典型的準靜態(tài)實驗條件10-3/s時,微觀機制完全不同,單根位錯通過晶界類擴散重排獨立發(fā)射。這對經(jīng)典分子動力學圖像提出了巨大挑戰(zhàn)。當使用分子動力學解析固體塑性機制時,一定要留意時間尺度效應,不然可能誤導對微觀機制的認知,在圖5所示的變形圖譜上錯誤定位,極大誤導理性力學分析。
圖10 晶界位錯形核機制的時間尺度效應[52]
元動力學和自適應加速分子動力學本質(zhì)上是借助統(tǒng)計力學的基本原理,在時間閾上的粗?;僮?強烈依賴于過渡態(tài)理論的適用性和稀有事件假設。與之平行的最后一類方法為擴展系綜分子動力學模擬,包括副本交換分子動力學(Parallel tempering or Multi-replica/replica exchange)[53]和溫度加速分子動力學(Temperature -accelerated)[54-58]等。這類方法利用聯(lián)合的多個子系綜,每個系綜的溫度或某個物理量不同,而其他物理條件等價,無論哪個子系綜有事件發(fā)生,都可以記入總事件數(shù)量并累積時間,從而使串行的時間序列事件實現(xiàn)并行化,達到直接加速分子動力學模擬的目的。但該類方法強烈依賴于計算機并行規(guī)模,加速效率不高,如1000個處理器并行,原則上可以加速3個數(shù)量級的時間尺度,從而將典型的MD納秒尺度延長到微秒,但是與實驗室時間尺度仍有不小差距。
分子動力學的終極目標為實現(xiàn)量子力學精度-大規(guī)模-實驗室時間尺度的原子模擬。原則上,該目標可以通過如量子霸權基于從頭算(ab initio)實現(xiàn)[59,60]。但是,目前最先進的超級計算機,通過大規(guī)模并行才可以實現(xiàn)十億量級、微秒尺度的直接分子動力學模擬[16-18]。這樣的空間規(guī)模以及合適的邊界條件處理,基本上可以囊括材料中的各種局域或擴展缺陷,包括位錯、界面和析出相等[61],從而使模擬的空間尺度限制大大降低。所以,分子動力學面臨的最艱巨的挑戰(zhàn)為勢函數(shù)的準確性與時間尺度問題。
近些年,人工智能和機器學習結合材料基因工程思想和數(shù)據(jù)科學,為材料科學提供了新的研究手段,甚至已經(jīng)并列為實驗、理論和計算外的第四大研究范式。目前,分子動力學可以模擬的體系受限于勢函數(shù)本身的準確性,對于高溫合金和高熵合金等多組元體系的模擬受到很大的限制。為解決多元復雜體系勢函數(shù)構建的普遍問題,Behler等[62]創(chuàng)新性地提出了基于高通量計算和機器學習構建具有量子力學精度的神經(jīng)網(wǎng)絡勢函數(shù)的普適方案,從而有望擺脫多元勢函數(shù)的限制。鄂維南等[63,64]已經(jīng)基于此策略,實現(xiàn)了億量級和量子力學精度的直接分子動力學模擬金屬材料力學行為,為此獲得了2020年度的戈登貝爾獎。所以,基于機器學習與高通量數(shù)據(jù)庫構建,建立任意元素組合的神經(jīng)網(wǎng)絡勢函數(shù),采用深度學習模型預測構型演化,有望實現(xiàn)化學組分自由和量子精度的分子動力學模擬,從而可以極大地拓展分子動力學在力學、材料、物理、化學和生物等多學科中的應用空間。
分子動力學仍然受限于模擬的時間尺度問題,這可能是原子模擬最后的難點所在。雖然元動力學、勢能面搜索與優(yōu)化方法在近20年取得了重要進展,但是仍有諸多前沿問題需要解決。如集體變量(CV)的普適性問題,模擬加速的效率深度取決于選擇的CV能否代表最慢的物理過程。而對于復雜體系,選取合適的CV本身就構成了不小的挑戰(zhàn)。所以,不依賴于CV的加速方法應運而生。Noé等[65]基于深度學習和統(tǒng)計力學提出了CV自由的玻爾茲曼生成器(Boltzmann generators)方法,可以跳過分子動力學在平衡態(tài)構型附近的長時間取樣,而直接從固有結構出發(fā),生成一系列代表平衡態(tài)性質(zhì)的龐大等價構型,通過深度學習驅(qū)使構型能量符合正則(玻爾茲曼)分布?;趦?yōu)化的配分函數(shù)構建自由能形式,并發(fā)現(xiàn)新的構型,從而實現(xiàn)復雜結構含時演化,而不需要任何先驗的反應坐標知識。
最后,Invernizzi等[66]提出了增強抽樣的統(tǒng)一方案,將元動力學與擴展系綜模擬有機結合,通過副本交換等技術,將不同系綜中不同集體變量相關的驅(qū)動勢和不同的抽樣方法統(tǒng)一到整個擴展系綜中,實現(xiàn)全局有效加速。這些進展都為最終解決分子動力學的時間尺度問題奠定了堅實的基礎。