苗樹偉,解德欣,劉建鵬,葉中豪,郭 峰,*
(1.聊城大學(xué)物理科學(xué)與信息工程學(xué)院,聊城 252000;2.山東省光通信科學(xué)與技術(shù)重點(diǎn)實(shí)驗(yàn)室,聊城 252000)
使用計(jì)算機(jī)模擬,來實(shí)現(xiàn)一些實(shí)驗(yàn)中不容易監(jiān)測或一些試驗(yàn)成本過高的狀態(tài)或過程(例如爆炸和高溫燃燒過程,亦或是一些反應(yīng)非常迅速缺乏測定技術(shù)的實(shí)驗(yàn))??蒲腥藛T通過對微觀世界的原子與分子的動(dòng)態(tài)模擬,來探索新材料的物理特性及其性能,為物理與材料科學(xué)的發(fā)展奠定基礎(chǔ)。目前使用較多的研究方法為分子力學(xué)模擬與量子力學(xué)計(jì)算,量子力學(xué)的計(jì)算是可以不依靠那些基于實(shí)驗(yàn)所得的數(shù)據(jù),而直接進(jìn)行相關(guān)的運(yùn)算并獲得物質(zhì)的幾何結(jié)構(gòu)以及電子結(jié)構(gòu)等相關(guān)的物質(zhì)結(jié)構(gòu),由此也能夠獲得物質(zhì)所具有的能量及其物理性質(zhì)[1]。分子力學(xué)則是通過對體系進(jìn)行相關(guān)的數(shù)學(xué)建模,并對模型中的粒子相互作用抽象化來實(shí)現(xiàn),模型中包含了物理、化學(xué)、力學(xué)等特性。最后使用分子動(dòng)力學(xué)模擬(molecule dynamics)等一系列計(jì)算機(jī)運(yùn)算手段,分析微觀體系的運(yùn)動(dòng)軌跡以及結(jié)構(gòu),得出宏觀狀態(tài)下該模型所擁有的物理化學(xué)特性。
分子動(dòng)力學(xué)模擬被設(shè)定為遵守經(jīng)典牛頓定律,體系中的所有粒子都在力場中遵循經(jīng)典牛頓定律進(jìn)行運(yùn)動(dòng)。且體系中的粒子都遵守疊加原理[2?3]。通過計(jì)算系統(tǒng)中粒子在多個(gè)離散的時(shí)刻下的牛頓運(yùn)動(dòng)方程就可以得到粒子在這些時(shí)刻下的運(yùn)動(dòng)信息。
雖然量子力學(xué)(QM)方法能夠?yàn)檠芯刻峁┖芏嘤杏玫闹笇?dǎo),但也存在部分的缺陷,一旦系統(tǒng)想要考慮整體的運(yùn)動(dòng)演化模擬,就需要相當(dāng)大的硬件系統(tǒng)資源[4]。其他的對于經(jīng)驗(yàn)原子間勢能雖然可以模擬動(dòng)態(tài)的演化,但也需要初始化原子間的連接,對于反應(yīng)的模擬也存在著一定的缺陷。因此ReaxFF 就被開發(fā),走入人們的視野中[5]。
ReaxFF 是以鍵級(jí)為基礎(chǔ),鍵級(jí)這一概念是由分子軌道理論提出的,能夠描述原子間的相互作用以及原子間鍵的生成與斷裂,原子與原子的相互作用不再是只取決于自身的所處位置,而是取決于和該原子相互影響的所有原子[6]。ReaxFF 做到了在保證精度不下降的同時(shí)降低了計(jì)算的復(fù)雜度,更好地加快了計(jì)算速度。當(dāng)我們把ReaxFF 與MD 相結(jié)合后,便可應(yīng)用于更加龐大的模擬體系,并且不需要預(yù)設(shè)初始的反應(yīng)狀態(tài)等條件。ReaxFF反應(yīng)力場已經(jīng)在高能材料、有機(jī)小分子體系中廣泛地使用。ReaxFF 反應(yīng)力場能夠在不多的計(jì)算資源下很好地表現(xiàn)出物理與化學(xué)的運(yùn)動(dòng)及反應(yīng)過程。
ReaxFF反應(yīng)力場與經(jīng)驗(yàn)非反作用力場類似,反作用力場將系統(tǒng)能量劃分為不同的局部能量貢獻(xiàn)[7]。公式為
反應(yīng)力場的核心是鍵級(jí)Ebond,將原子間的關(guān)系建模為一個(gè)關(guān)于鍵級(jí)的函數(shù),準(zhǔn)確描述原子間鍵所蘊(yùn)含的能量。函數(shù)中包含鍵能、范德華勢能等物理化學(xué)能量。Eangle和Eetor為三體價(jià)角應(yīng)變與四體鈕角應(yīng)變所蘊(yùn)含的相關(guān)能量。Eover是一種能量懲罰,它阻止原子的超配位,這是基于原子的價(jià)態(tài)規(guī)則(例如,如果一個(gè)碳原子形成超過4 個(gè)鍵,則施加一個(gè)嚴(yán)格的能量懲罰)[8]。Ecol和Evdw是計(jì)算所有原子之間的靜電和色散貢獻(xiàn),而不考慮連通性和鍵序。反應(yīng)力場中最小的單位是原子,鍵級(jí)包含了分子內(nèi)各原子間、各部分的能量。鍵級(jí)的計(jì)算函數(shù)為
式中:rij為原子間距離;Pbo是經(jīng)驗(yàn)參數(shù)。式中BO'σij、BO'πij、BO'ππij則依次表示單鍵σ、雙鍵π、三鍵ππ 的鍵級(jí),σ、π、ππ 鍵之間的躍遷不包含不連續(xù)。r0為平衡鍵長。由于鍵級(jí)還未經(jīng)過矯正,加上各元素的成鍵也不盡相同,因此會(huì)產(chǎn)生很多不合理或者較弱的成鍵,為了解決這一問題則需要引入修正項(xiàng)。經(jīng)過一系列的運(yùn)算,最終得到的總鍵能公式為
其中:P和De是能量系。
本研究在ReaxFF已有的基礎(chǔ)上對ReaxFF反應(yīng)力場進(jìn)行了擴(kuò)展,編寫了一套名為Irff 的代碼。Irff 的代碼使用了Facebook 人工智能研究院(FAIR)基于Torch推出的PyTorch 庫。PyTorch 提供了更靈活的API,可以定義動(dòng)態(tài)計(jì)算圖,且為張量上的所有操作提供了自動(dòng)微分(也稱為反向傳播)[9]。因此我們可以使用Adam 優(yōu)化器,Adam 算法可以為不同的權(quán)重參數(shù)分配不同的自適應(yīng)學(xué)習(xí)率[10],如此進(jìn)一步提升了代碼的優(yōu)化能力。并且基于PyTorch的tensor為矩陣的特性,我們對ReaxFF 原有的函數(shù)進(jìn)行了改進(jìn),將其改為矩陣運(yùn)算。經(jīng)過更改矩陣運(yùn)算后的公式如下,其中BO'是未經(jīng)修正的鍵級(jí)。
Etor是表示為二面角勢能的鍵級(jí)函數(shù),與鍵級(jí)呈正相關(guān)變化,最小值為0。sin Θijk以及sin Θijk均為設(shè)置的控制函數(shù)在價(jià)角180°時(shí)保證二面角的能量是0。Etor具有多極值點(diǎn)。
Econj為描述四體作用共軛情況的四體共軛項(xiàng):
Eval為附加能量,只有當(dāng)價(jià)角處在唯一的最優(yōu)角度時(shí)Eval才為0。Eval與鍵級(jí)呈現(xiàn)正相關(guān)的關(guān)系,并且是連續(xù)的。σ、π 鍵形成的價(jià)角能量各不相同[11]。
我們從LAMMPS程序包中選取了一部分相應(yīng)的力場文件,并重建一個(gè)與其研究過程中建立的模型相近似的模型。使用編寫的軟件進(jìn)行計(jì)算且將計(jì)算結(jié)果與通過GULP(general utility lattice pro?gram)軟件計(jì)算所得的結(jié)果進(jìn)行比對,以此來證明并保證本軟件的計(jì)算結(jié)果是正確可靠的。
我們參考了文獻(xiàn)[12]構(gòu)建了C2H4 模型,參考了文獻(xiàn)[13]構(gòu)建了C5H8O4N12模型,參考了文獻(xiàn)[14]構(gòu)建了C8H10 模型,以及文獻(xiàn)[15]構(gòu)建了H6CF2 模型。Cell 的晶格向量分別設(shè)置為10 × 10 × 10 以及20 × 20 × 20,且具有周期性邊界,詳情見圖1。所有結(jié)構(gòu)建模都采用球棍模型,小球中灰色、紅色、藍(lán)色、綠色分別對應(yīng)的元素為碳、氧、氮、氟。后續(xù)小節(jié)我們將會(huì)采用上述文獻(xiàn)中所用的力場文件以及模型來驗(yàn)證我們的程序。
圖1 測試本程序設(shè)計(jì)所用的模型
分子動(dòng)力學(xué)模擬(MD)是根據(jù)所有分子的當(dāng)前坐標(biāo)計(jì)算分子的受力,根據(jù)受力更新分子的坐標(biāo),在此過程中收集用于計(jì)算宏觀性質(zhì)的有關(guān)信息。本文中MD 模擬Time Step 設(shè)置為0.1,分子動(dòng)力學(xué)描述的運(yùn)動(dòng)是不連續(xù)的以時(shí)間原子當(dāng)時(shí)所有的各項(xiàng)物理參數(shù)。步數(shù)設(shè)置為100。我們使用了NVE 模擬來生成訓(xùn)練所需的數(shù)據(jù),而Velocity Verlet 是實(shí)現(xiàn)NVE 集成的唯一動(dòng)力。步長對于動(dòng)力學(xué)來說是必須的,太小會(huì)造成系統(tǒng)資源的浪費(fèi),太大會(huì)出現(xiàn)“爆炸”。牛頓第二定律的直接積分保留了系統(tǒng)總能量等的參數(shù),因此Velocity Verlet算法是最合適的,在步長大的時(shí)候它依然可以使總能量有不俗的穩(wěn)定表現(xiàn)。需要注意的是NVE模擬中溫度一般保持動(dòng)態(tài)的恒定,如果結(jié)構(gòu)產(chǎn)生顯著的變化那么溫度有可能會(huì)發(fā)生急劇變化。相較于Runge?Kutta 等的算法,它們確實(shí)能在短期能量保存效果上獲得不俗的表現(xiàn),但長時(shí)間的模擬則會(huì)產(chǎn)生能量緩慢的漂移。
在分子動(dòng)力學(xué)模擬中力場(FFs)的精度是核心要素,力場文件中包含了分子內(nèi)、分子間的勢能面(PES)等各種參數(shù)。因此我們從LAMMPS中選取了文獻(xiàn)[12]中使用的力場參數(shù)文件,然后構(gòu)建乙烯(PE)結(jié)構(gòu)。我們對此PE 結(jié)構(gòu)進(jìn)行分子動(dòng)力學(xué)模擬并產(chǎn)生一個(gè)有一百幀的運(yùn)動(dòng)軌跡文件md.traj。然后直接讀取該軌跡文件,利用IRFF 函數(shù)計(jì)算出相應(yīng)的力場參數(shù)。我們分別對比了Ebond鍵能項(xiàng)、Eover過配位修正項(xiàng)、Eunder配位能校正項(xiàng)、Eang鍵角彎曲能、Ecoul靜電能、Evdw范德華能,以及總能量Total?Energy。
其中:θ為鍵角;θ0為平衡鍵角;kb為鍵角彎折力常數(shù)。二面角扭轉(zhuǎn)能。
其中:Vn為勢壘高度;N為多重度;ω為扭轉(zhuǎn)角度;γ為相因子。由于不同的方法計(jì)算得出的能量不盡相同,因此只有在同體系的情況下其他構(gòu)象計(jì)算得到的能量才有比較意義。
另外,我們還從LAMMPS 中選取了文獻(xiàn)[13?15]中優(yōu)化使用的力場文件,進(jìn)行了與上述相同的操作與計(jì)算。得到了四組IRFF 與GULP的參數(shù)對比,詳見圖2~圖5。
圖2 結(jié)構(gòu)為C2H4時(shí)相關(guān)參數(shù)對比
圖3 結(jié)構(gòu)為CHON時(shí)對比
圖4 結(jié)構(gòu)為CH4O2時(shí)對比
圖5 結(jié)構(gòu)為H6C3F2時(shí)對比
上述對比可以看出,經(jīng)過我們矩陣化運(yùn)算的轉(zhuǎn)換,原子間的各項(xiàng)能量均保持在正常水平范圍內(nèi)。保證準(zhǔn)確度不下降的情況下,我們的軟件IRFF可以詳細(xì)地展示出分子運(yùn)動(dòng)模擬中各參數(shù)的具體數(shù)值以及變化情況。將一個(gè)反應(yīng)力場運(yùn)動(dòng)模擬的細(xì)節(jié)直觀地展示在了研究人員的面前。
本文通過對ReaxFF 反應(yīng)力場代碼的改進(jìn)構(gòu)建,改變原有的計(jì)算方式為矩陣運(yùn)算。能夠直接使用分子動(dòng)力學(xué)的軌跡文件作為輸入數(shù)據(jù)集,以及引入PyTorch 模塊實(shí)現(xiàn)了能夠自動(dòng)微分的反應(yīng)動(dòng)力學(xué)模型。直白地將動(dòng)力學(xué)運(yùn)動(dòng)模擬過程中的各種勢能函數(shù)展現(xiàn)在科研人員的面前,為了解運(yùn)動(dòng)模擬過程中的詳細(xì)情況提供了理論支持。且經(jīng)過對比驗(yàn)證,證明了改進(jìn)后的代碼的精確性與可靠性。