基于自動(dòng)微分的反應(yīng)動(dòng)力學(xué)程序設(shè)計(jì)

苗樹偉，解德欣，劉建鵬，葉中豪，郭 峰,*

（1.聊城大學(xué)物理科學(xué)與信息工程學(xué)院，聊城 252000；2.山東省光通信科學(xué)與技術(shù)重點(diǎn)實(shí)驗(yàn)室，聊城 252000）

0 引言

使用計(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)過程。

1 模擬方法

1.1 ReaxFF反應(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是能量系。

1.2 軟件介紹

本研究在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］。

2 結(jié)果與討論

我們從LAMMPS程序包中選取了一部分相應(yīng)的力場文件，并重建一個(gè)與其研究過程中建立的模型相近似的模型。使用編寫的軟件進(jìn)行計(jì)算且將計(jì)算結(jié)果與通過GULP（general utility lattice pro?gram）軟件計(jì)算所得的結(jié)果進(jìn)行比對，以此來證明并保證本軟件的計(jì)算結(jié)果是正確可靠的。

2.1 模型構(gòu)建

我們參考了文獻(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ì)所用的模型

2.2 MD模擬

分子動(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)生能量緩慢的漂移。

2.3 實(shí)驗(yàn)具體細(xì)節(jié)

在分子動(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é)直觀地展示在了研究人員的面前。

3 結(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)后的代碼的精確性與可靠性。