基于自動微分的反應動力學程序設計

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

（1.聊城大學物理科學與信息工程學院，聊城 252000；2.山東省光通信科學與技術重點實驗室，聊城 252000）

0 引言

使用計算機模擬，來實現一些實驗中不容易監測或一些試驗成本過高的狀態或過程（例如爆炸和高溫燃燒過程，亦或是一些反應非常迅速缺乏測定技術的實驗）。科研人員通過對微觀世界的原子與分子的動態模擬，來探索新材料的物理特性及其性能，為物理與材料科學的發展奠定基礎。目前使用較多的研究方法為分子力學模擬與量子力學計算，量子力學的計算是可以不依靠那些基于實驗所得的數據，而直接進行相關的運算并獲得物質的幾何結構以及電子結構等相關的物質結構，由此也能夠獲得物質所具有的能量及其物理性質［1］。分子力學則是通過對體系進行相關的數學建模，并對模型中的粒子相互作用抽象化來實現，模型中包含了物理、化學、力學等特性。最后使用分子動力學模擬（molecule dynamics）等一系列計算機運算手段，分析微觀體系的運動軌跡以及結構，得出宏觀狀態下該模型所擁有的物理化學特性。

分子動力學模擬被設定為遵守經典牛頓定律，體系中的所有粒子都在力場中遵循經典牛頓定律進行運動。且體系中的粒子都遵守疊加原理［2?3］。通過計算系統中粒子在多個離散的時刻下的牛頓運動方程就可以得到粒子在這些時刻下的運動信息。

雖然量子力學（QM）方法能夠為研究提供很多有用的指導，但也存在部分的缺陷，一旦系統想要考慮整體的運動演化模擬，就需要相當大的硬件系統資源［4］。其他的對于經驗原子間勢能雖然可以模擬動態的演化，但也需要初始化原子間的連接，對于反應的模擬也存在著一定的缺陷。因此ReaxFF 就被開發，走入人們的視野中［5］。

ReaxFF 是以鍵級為基礎，鍵級這一概念是由分子軌道理論提出的，能夠描述原子間的相互作用以及原子間鍵的生成與斷裂，原子與原子的相互作用不再是只取決于自身的所處位置，而是取決于和該原子相互影響的所有原子［6］。ReaxFF 做到了在保證精度不下降的同時降低了計算的復雜度，更好地加快了計算速度。當我們把ReaxFF 與MD 相結合后，便可應用于更加龐大的模擬體系，并且不需要預設初始的反應狀態等條件。ReaxFF反應力場已經在高能材料、有機小分子體系中廣泛地使用。ReaxFF 反應力場能夠在不多的計算資源下很好地表現出物理與化學的運動及反應過程。

1 模擬方法

1.1 ReaxFF反應力場

ReaxFF反應力場與經驗非反作用力場類似，反作用力場將系統能量劃分為不同的局部能量貢獻［7］。公式為

反應力場的核心是鍵級Ebond，將原子間的關系建模為一個關于鍵級的函數，準確描述原子間鍵所蘊含的能量。函數中包含鍵能、范德華勢能等物理化學能量。Eangle和Eetor為三體價角應變與四體鈕角應變所蘊含的相關能量。Eover是一種能量懲罰，它阻止原子的超配位，這是基于原子的價態規則（例如，如果一個碳原子形成超過4 個鍵，則施加一個嚴格的能量懲罰）［8］。Ecol和Evdw是計算所有原子之間的靜電和色散貢獻，而不考慮連通性和鍵序。反應力場中最小的單位是原子，鍵級包含了分子內各原子間、各部分的能量。鍵級的計算函數為

式中：rij為原子間距離；Pbo是經驗參數。式中BO'σij、BO'πij、BO'ππij則依次表示單鍵σ、雙鍵π、三鍵ππ 的鍵級，σ、π、ππ 鍵之間的躍遷不包含不連續。r0為平衡鍵長。由于鍵級還未經過矯正，加上各元素的成鍵也不盡相同，因此會產生很多不合理或者較弱的成鍵，為了解決這一問題則需要引入修正項。經過一系列的運算，最終得到的總鍵能公式為

其中：P和De是能量系。

1.2 軟件介紹

本研究在ReaxFF已有的基礎上對ReaxFF反應力場進行了擴展，編寫了一套名為Irff 的代碼。Irff 的代碼使用了Facebook 人工智能研究院（FAIR）基于Torch推出的PyTorch 庫。PyTorch 提供了更靈活的API，可以定義動態計算圖，且為張量上的所有操作提供了自動微分（也稱為反向傳播）［9］。因此我們可以使用Adam 優化器，Adam 算法可以為不同的權重參數分配不同的自適應學習率［10］，如此進一步提升了代碼的優化能力。并且基于PyTorch的tensor為矩陣的特性，我們對ReaxFF 原有的函數進行了改進，將其改為矩陣運算。經過更改矩陣運算后的公式如下，其中BO'是未經修正的鍵級。

Etor是表示為二面角勢能的鍵級函數，與鍵級呈正相關變化，最小值為0。sin Θijk以及sin Θijk均為設置的控制函數在價角180°時保證二面角的能量是0。Etor具有多極值點。

Econj為描述四體作用共軛情況的四體共軛項：

Eval為附加能量，只有當價角處在唯一的最優角度時Eval才為0。Eval與鍵級呈現正相關的關系，并且是連續的。σ、π 鍵形成的價角能量各不相同［11］。

2 結果與討論

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

2.1 模型構建

我們參考了文獻［12］構建了C2H4 模型，參考了文獻［13］構建了C5H8O4N12模型，參考了文獻［14］構建了C8H10 模型，以及文獻［15］構建了H6CF2 模型。Cell 的晶格向量分別設置為10 × 10 × 10 以及20 × 20 × 20，且具有周期性邊界，詳情見圖1。所有結構建模都采用球棍模型，小球中灰色、紅色、藍色、綠色分別對應的元素為碳、氧、氮、氟。后續小節我們將會采用上述文獻中所用的力場文件以及模型來驗證我們的程序。

圖1 測試本程序設計所用的模型

2.2 MD模擬

分子動力學模擬（MD）是根據所有分子的當前坐標計算分子的受力，根據受力更新分子的坐標，在此過程中收集用于計算宏觀性質的有關信息。本文中MD 模擬Time Step 設置為0.1，分子動力學描述的運動是不連續的以時間原子當時所有的各項物理參數。步數設置為100。我們使用了NVE 模擬來生成訓練所需的數據，而Velocity Verlet 是實現NVE 集成的唯一動力。步長對于動力學來說是必須的，太小會造成系統資源的浪費，太大會出現“爆炸”。牛頓第二定律的直接積分保留了系統總能量等的參數，因此Velocity Verlet算法是最合適的，在步長大的時候它依然可以使總能量有不俗的穩定表現。需要注意的是NVE模擬中溫度一般保持動態的恒定，如果結構產生顯著的變化那么溫度有可能會發生急劇變化。相較于Runge?Kutta 等的算法，它們確實能在短期能量保存效果上獲得不俗的表現，但長時間的模擬則會產生能量緩慢的漂移。

2.3 實驗具體細節

在分子動力學模擬中力場（FFs）的精度是核心要素，力場文件中包含了分子內、分子間的勢能面（PES）等各種參數。因此我們從LAMMPS中選取了文獻［12］中使用的力場參數文件，然后構建乙烯（PE）結構。我們對此PE 結構進行分子動力學模擬并產生一個有一百幀的運動軌跡文件md.traj。然后直接讀取該軌跡文件，利用IRFF 函數計算出相應的力場參數。我們分別對比了Ebond鍵能項、Eover過配位修正項、Eunder配位能校正項、Eang鍵角彎曲能、Ecoul靜電能、Evdw范德華能，以及總能量Total?Energy。

其中：θ為鍵角；θ0為平衡鍵角；kb為鍵角彎折力常數。二面角扭轉能。

其中：Vn為勢壘高度；N為多重度；ω為扭轉角度；γ為相因子。由于不同的方法計算得出的能量不盡相同，因此只有在同體系的情況下其他構象計算得到的能量才有比較意義。

另外，我們還從LAMMPS 中選取了文獻［13?15］中優化使用的力場文件，進行了與上述相同的操作與計算。得到了四組IRFF 與GULP的參數對比，詳見圖2～圖5。

圖2 結構為C2H4時相關參數對比

圖3 結構為CHON時對比

圖4 結構為CH4O2時對比

圖5 結構為H6C3F2時對比

上述對比可以看出，經過我們矩陣化運算的轉換，原子間的各項能量均保持在正常水平范圍內。保證準確度不下降的情況下，我們的軟件IRFF可以詳細地展示出分子運動模擬中各參數的具體數值以及變化情況。將一個反應力場運動模擬的細節直觀地展示在了研究人員的面前。

3 結語

本文通過對ReaxFF 反應力場代碼的改進構建，改變原有的計算方式為矩陣運算。能夠直接使用分子動力學的軌跡文件作為輸入數據集，以及引入PyTorch 模塊實現了能夠自動微分的反應動力學模型。直白地將動力學運動模擬過程中的各種勢能函數展現在科研人員的面前，為了解運動模擬過程中的詳細情況提供了理論支持。且經過對比驗證，證明了改進后的代碼的精確性與可靠性。