基于FPGA的分子動力學模擬動態數據流分析

吳子剛 柴志雷

江南大學人工智能與計算機學院 江蘇 無錫 214000

引言

分子動力學（MolecularDynamics，MD）模擬是過去幾十年科研人員從微觀入手研究宏觀性質的重要手段之一。然而，由于MD算法的計算復雜度，在傳統處理器上的一個簡單生物實體模擬通常需要數周甚至數月的時間。因此對于MD模擬的加速就顯得很是重要。GROMACS 能很好地利用從 SSE到 AVX512 的 X86 架構的 SIMD 指令集，從而獲得較大的性能提升[1]; AMBER 則在 NVIDIA 的 GPU 上進行了大量的優化工作，成為用 GPU 模擬生物體系的代表作之一[2]。

MD模擬的關鍵路徑是粒子間的作用力計算，具有很高的并行性。在眾多的加速平臺中，FPGA以低功耗和高能效而成為一種具有吸引力的解決方案。因此很多研究者也都嘗試將FPGA的優勢引入分子動力學模擬中。例如，S.Kasap等[3]人首次嘗試將LAMMPS移植到基于FPGA的并行計算機上，取得了顯著的性能提升。但是傳統RTL開發的較高編程難度和設計時間，目前主流的MD模擬軟件包還沒有采用FPGA加速。因此，本文使用HLS探尋加速分子動力學模擬中關鍵路徑的一些優化手段，并分析得出目前HLS還不支持動態數據流行為的描述。

1 分子動力學模擬

隨著計算機的發展和算力的提高，科學研究的手段不再局限于宏觀實驗。當宏觀實驗難以進行或者很難通過宏觀實驗來得到滿意結果時，微觀模擬為科學研究提供了一個嶄新的角度和方向。蒙特卡羅計算法（MC計算法）是最早對龐大系統采用非量子計算的微觀模擬方法。但MC的缺點也很明顯，即只能計算統計的平均值，無法得到系統的動態信息。因此目前更多地采用的是分子動力學模擬。

1.1 MD關鍵路徑

MD模擬是在牛頓力學下分子或原子集合中的迭代應用，每次迭代中分為兩個過程：①計算粒子間的受力；②更新粒子之間的速度和位置。計算的力跟模擬系統的粒子構成有關：

其中，公式中的前三項為鍵成力，最后一項為非鍵成力。非鍵成力因為引入了截斷半徑又分為長程非鍵成力和短程非鍵成力（RL），而RL占了整個MD的90%FLOPS，是MD中的關鍵路徑。算法1說明了RL的計算流程。

算法1RL計算偽代碼

foreach cell do

foreach neighbor cell do

foreach atom in home cell do

force = 0;

foreach atom in neighbor cell doif(distance() ＜ cutoff radius)

compute force();

end foreach

end foreach

end foreach

end foreach

RL又包含了LJ勢能作用力和短程庫侖力。兩者都是粒子對之間的作用力。為了簡化問題突出本質，以液氬（不帶電）為研究對象。因此只需要考慮LJ勢能作用力。

LJ勢能作用力可以寫成公式（2）形式：

1.2 晶胞法

在分子動力學模擬中通常采用截斷半徑來將降低非鍵成力的計算復雜度，即只計算中心粒子和截斷半徑以內粒子之間的力。在引入截斷半徑后，有兩種方法來進行鄰域索引：鄰接列表法和晶胞法。

鄰接列表法的做法是對每個粒子都將其截斷半徑內的所有粒子保存在這個粒子的鄰接表中。這樣做雖然可以獲得100%的效率，但是鄰接列表的構造會消耗大量的存儲和時間。尤其是截斷半徑很大的情況下，每個粒子的鄰接表中會有上千個粒子。

晶胞法是將模擬空間劃分成邊長一樣的三維晶胞。一般截斷半徑等于晶胞邊長，使用晶胞法的效率為15.5%。這樣粒子只需要跟所在晶胞的其他粒子以及相近晶胞中的粒子進行力的計算。同時晶胞的構造只需要遍歷一遍模擬空間的位置信息。

因為晶胞法所體現出的空間局部性，結合FPGA的有限片上資源，晶胞法成為FPGA加速MD的首選方法。

2 HLS介紹

FPGA是一個可編程器件。它由查找表（LUT）、寄存器（FF）、連線資源以及I/O等組成。傳統的RTL級別開發使用的是HDL硬件描述語言，主要分為VHDL、verilog以及system Verilog。HLS以FPGA為執行結構，使軟件工程師能夠優化代碼的整體性能、功耗和延遲，而不需要解決單個內存空間和有限計算資源的性能瓶頸。

3 RL計算的HLS加速

目前HLS在縮短加速器設計時間和工程成本方面取得了巨大的成功。然而，在HLS中將C/C++程序映射為想要的硬件架構仍然需要開發人員具有足夠的FPGA架構知識和優化技術。本文針對MD關鍵路徑RL的計算，討論了HLS下的實現手段，分析了RL計算中的動態數據流行為。

3.1 HLS加速

對于RL的計算可以簡單分為三個模塊：①粒子數據預取模塊；②粒子間距離計算模塊；③計算RL力的模塊。為生成HLS加速器，首先考慮對圖一中算法的最內層循環進行流水線設計，包括粒子間距離的計算和力的計算兩個模塊。此時的計算結構如圖1（a）。

圖1 RL計算架構

圖1 （a）中距離計算模塊和力計算模塊雖然都經過了流水線設計，但都達不到很高的效率。因為采用的是晶胞法，粒子距離計算模塊的效率僅為15.5%。粒子距離計算模塊每個時鐘周期都會從粒子數據緩存模塊中讀取一個粒子對，但不會每個周期都輸出粒子對，從而使力計算模塊的效率也僅為15.5%。

3.2 力計算模塊的優化

粒子i和粒子j的相互作用力可以由公式（2）簡化成公式（3）的形式。

可以看到需要計算 r-8和 r-14兩個高次冪。為了減少DSP的消耗，使FPGA芯片上運行更多的流水線，本文采用插值的方式來得到這兩個高次冪的值。先沿著X軸把函數分成幾段，并使后一段為前一段的2倍，然后在每段中等間距地取得一些離散點，得到相應的區間。這個過程被稱為插值，在一段中插入的離散點越多，采用插值方式獲得的高次冪就越精確。

為了避免開方操作，使用r2作為X軸。計算中如果r2的值落在由相鄰離散點A和B構成的區間（A，B）之中，在使用一階線性插值時r-k的值可由公式（4）得到：

公式中的x就等于 r2，a就等于A點的橫坐標。參數C0和 C1預先計算好并分別存放在兩個查找表中。計算時使用 r2的值來生成查找表的地址。因為使用的是一階插值，所以在每段中需要插入足夠多的離散點，得到足夠多的區間。本文在每段中會有256個區間。

3.3 動態數據流

為此通常會為力計算模塊配備多個距離計算模塊，以使力計算模塊的效率達到接近100%，從而提升整體性能，如圖1（b）。在距離計算模塊和力計算模塊之間加入了一個仲裁模塊。仲裁模塊使用FIFO存儲距離計算模塊輸出的有效粒子對，并決定傳輸哪個距離計算模塊輸出的有效粒子對到力計算模塊。

由于距離計算模塊不是每個時鐘都會有輸出且輸出的時間是沒有規律的，圖1（b）中距離計算模塊、仲裁器和力計算模塊三者構成的數據路徑就是動態數據路徑。

雖然現在HLS工具在綜合靜態數據流和規則算法時會取得很好的效果，但是在綜合描述動態數據流行為的算法時通常都會是很差的效果。這是因為現在絕大多數HLS工具都假定了模塊每接收一個輸入就會產生一個輸出。在這種假定下就不會產生動態數據流行為。因此無法直接使用HLS很好地實現圖1（b）中描述的結構。

值得注意的是，除了上述原因外使用HLS還會出現全局狀態機（FSM）的問題。在圖1（a）中只生成了一個全局FSM，這導致了三個模塊之間的強制順序執行。圖1（b）如果只使用一個全局FSM，那么實現的系統顯然是低效的。理想的情況應該是各個模塊都有自己獨立的FSM來控制處理邏輯，而不是只使用一個HLS生成的全局FSM。這樣各個模塊都是并行運行的。

4 結束語

本文討論了在FPGA上利用HLS工具對于MD算法關鍵路徑的加速，發現關鍵路徑中存在動態數據流行為，分析得到了這限制了加速器的性能。