基于FPGA的分子動力學(xué)模擬動態(tài)數(shù)據(jù)流分析

吳子剛 柴志雷

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

引言

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

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

1 分子動力學(xué)模擬

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

1.1 MD關(guān)鍵路徑

MD模擬是在牛頓力學(xué)下分子或原子集合中的迭代應(yīng)用，每次迭代中分為兩個過程：①計算粒子間的受力；②更新粒子之間的速度和位置。計算的力跟模擬系統(tǒng)的粒子構(gòu)成有關(guān)：

其中，公式中的前三項為鍵成力，最后一項為非鍵成力。非鍵成力因為引入了截斷半徑又分為長程非鍵成力和短程非鍵成力（RL），而RL占了整個MD的90%FLOPS，是MD中的關(guān)鍵路徑。算法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勢能作用力和短程庫侖力。兩者都是粒子對之間的作用力。為了簡化問題突出本質(zhì)，以液氬（不帶電）為研究對象。因此只需要考慮LJ勢能作用力。

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

1.2 晶胞法

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

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

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

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

2 HLS介紹

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

3 RL計算的HLS加速

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

3.1 HLS加速

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

圖1 RL計算架構(gòu)

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

3.2 力計算模塊的優(yōu)化

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

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

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

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

3.3 動態(tài)數(shù)據(jù)流

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

由于距離計算模塊不是每個時鐘都會有輸出且輸出的時間是沒有規(guī)律的，圖1（b）中距離計算模塊、仲裁器和力計算模塊三者構(gòu)成的數(shù)據(jù)路徑就是動態(tài)數(shù)據(jù)路徑。

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

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

4 結(jié)束語

本文討論了在FPGA上利用HLS工具對于MD算法關(guān)鍵路徑的加速，發(fā)現(xiàn)關(guān)鍵路徑中存在動態(tài)數(shù)據(jù)流行為，分析得到了這限制了加速器的性能。