基于FPGA的分子動力學模擬匹配單元設計

王 鑫,劉 炫

(1. 江南大學 教育部輕工過程先進控制重點實驗室，江蘇 無錫 214122; 2. 江蘇省未來網絡創(chuàng)新研究院，南京 211111)

0 引言

分子動力學模擬利用計算機求解生物分子體系內原子和分子的運動方程來模擬這些粒子的運動軌跡，從而獲得系統(tǒng)的溫度、體積、壓力等宏觀和微觀過程量。分子動力學模擬在計算物理、計算化學、材料科學、生命科學、生物醫(yī)學等多個領域有著廣泛的應用，已成為與實驗同等重要的科學研究方法。

分子動力學模擬過程中，需要花費大量時間用于計算范德華相互作用和靜電相互作用，而模擬具有現實意義的生物大分子體系動力學過程往往需要幾個月甚至幾年時間。因此，針對生物大分子進行微秒甚至毫秒級別的模擬，對范德華相互作用和靜電相互作用的優(yōu)化加速將是首要任務，其中進行硬件加速是提高計算速度的重要方法。近年來，隨著集成電路技術的快速發(fā)展，分子動力學模擬在專用集成電路(ASIC, application specific integrated circuit)以及現場可編程門陣列(FPGA，field programmable gate array)上都得到了很好的實現，并在并行化、算法進展和硬件專業(yè)化方面進行了不同程度的改進。這些改進大大提升了系統(tǒng)的性能。需要指出的是，相比于ASIC，FPGA設計靈活，功能強大，風險低。因此，FPGA軟硬件協同設計環(huán)境下的硬件加速技術具有重要的研究意義。

在分子動力學模擬中，粒子間需要計算的力包括鍵合力和非鍵合力[15]，其中非鍵合力包括范德華力和靜電力。為了降低非鍵合力的計算復雜度，將非鍵合力分解為短程力與長程力，并對兩者分別設計專用的低復雜度計算方法。具體來說，對固定的一個粒子，可以定義一個截止半徑。這個粒子與范圍內的其它粒子之間的非鍵合力就稱為短程力。如何有效地從大量粒子中匹配出滿足短程力計算要求的粒子對，是本文要解決的問題。

為了提升分子動力學模擬中粒子間非鍵合力的計算效率，本文提出了匹配單元的設計方案，用以高效地篩選出滿足短程力計算條件的粒子對。

1 系統(tǒng)結構及原理

圖1是本文研究的短程力計算系統(tǒng)的總體框圖。首先，上位機通過PCIe(peripheral component interconnect express)總線將二進制的粒子數據傳遞到FPGA加速卡上的BRAM(Block RAM)。當數據預處理單元接收到開始讀取數據的指令時，就從BRAM中讀取數據。讀取完成后，數據預處理單元將讀取的數據對應到各個粒子，再將粒子數據送入匹配單元進行匹配，匹配成功的粒子對進入到短程力計算單元。最后，計算得到的力將在力累加器中進行累加，并通過力回寫單元將累加結果寫回BRAM，供上位機取回。

圖1 短程力計算系統(tǒng)總體框圖

綜上所述，短程力計算單元是短程力計算系統(tǒng)的核心部分。為了保證短程力計算單元的利用率，就需要從預處理后的粒子數據中準確快速地匹配出符合短程力計算條件的粒子對。針對該問題，本文設計了匹配單元。

本文設計的匹配單元是根據截止半徑來篩選粒子對的。在半徑為的圓球內，以球心處粒子為參考粒子，匹配單元將濾除圓球外的粒子。與此同時，將圓內的粒子依次和參考粒子配對送入短程力計算單元。如果匹配單元每次只對一個參考粒子進行匹配，匹配效率不高。若多個參考粒子和周圍粒子同時進行匹配，就能大幅度提升匹配效率。此外，分子動力學模擬是通過離散時間步長來模擬化學系統(tǒng)中單個粒子的運動[21]，而粒子在化學系統(tǒng)中的空間分布具有不確定性，因此，參考粒子的選取方式極大程度上影響著匹配單元的工作效率。

為保證匹配單元的工作效率，避免參考粒子的重復選取，需要對分子系統(tǒng)進行區(qū)域劃分。因為化學系統(tǒng)是在一個(非物理)周期性邊界條件盒子中的，所以，實際模擬的系統(tǒng)相當于將一個無限的空間在三維坐標軸上平鋪為大小相同的盒子。粒子就分布在這些盒子內。盒子的大小和粒子大小有關，每個盒子都有對應的盒子編號。除此之外，每個盒子內的粒子還有相應的盒內編號，并且每個粒子都有對應的位置坐標。為了便于遍歷計算粒子對之間的相互作用，將其中的某個粒子作為參考粒子，稱為靜態(tài)粒子。靜態(tài)粒子周圍的粒子則稱為動態(tài)粒子。由于化學系統(tǒng)內的粒子是隨機分布的，故盒內粒子數具有不確定性。

在設計中，本文每個盒子最多可容納64個粒子。為了補足盒內粒子數，使動態(tài)粒子必有靜態(tài)粒子與其進行相互作用，本文又將靜態(tài)粒子分為有效靜態(tài)粒子和無效靜態(tài)粒子(即空粒子)，空粒子是用來補足盒內粒子數。與此同時，為了方便判斷靜態(tài)粒子是不是空粒子，本文引入了空粒子標志位。

當對化學系統(tǒng)進行盒子劃分后，結合FPGA的DSP(digital signal processor)資源的限制，匹配單元的匹配原則就不單單是在截止半徑范圍內，就需要更細致的匹配原則。匹配原則的內容主要有以下兩個方面：

1)偏序法：偏序法主要是為了防止同一盒子內的兩個粒子重復計算短程力。具體而言，如果兩個粒子在同一個盒子內，若將其中一個粒子作為參考粒子，參考粒子同周圍粒子進行短程力的計算，同時盒內每個粒子都將輪流作為參考粒子，此時兩個粒子之間就進行了兩次短程力的計算。根據牛頓第三運動定律，兩次計算的短程力大小相等，完全不需要計算兩次。因此，參考粒子只需和盒內編號比自己大的粒子進行短程力的計算，盒內編號比自己小的已經計算過了，不需要再次計算。

|X1-X2|

(1)

|Y1-Y2|

(2)

|Z1-Z2|

(3)

(4)

(5)

(6)

(7)

2 硬件設計

匹配單元的設計如圖2所示。圖中的匹配單元有8個距離計算器，可實現1個動態(tài)粒子分別和8個靜態(tài)粒子進行配對。首先將盒子中的靜態(tài)粒子和動態(tài)粒子送入FIFO(first input first output)進行暫存，當匹配單元接收到數據預處理完成的信號后，FIFO中靜態(tài)粒子和動態(tài)粒子進入流寄存器中，以粒子對形式暫存，共計8組粒子對。一旦粒子對發(fā)送到距離計算器，就將粒子之間的距離與截止半徑進行比較，最后輸出粒子對匹配結果。距離計算器的匹配結果存放在向量寄存器(向量寄存器本質上是多個寄存器進行捆綁)中。由于可能存在空粒子，所以，需要排除空粒子與動態(tài)粒子生成的粒子對。于是，通過將距離計算器的匹配結果和流寄存器中的空粒子標志位進行與運算，經過與運算后的8組粒子對流向仲裁器，最后輸出滿足短程力計算條件的粒子對。

圖2 匹配單元設計框圖

值得注意的是，數據預處理單元的工作速度也限制匹配單元的工作效率。由于數據預處理的速度是未知的，因此需要一個FIFO，通過緩存數據來控制進入匹配單元的粒子。除此之外，粒子對在進行計算的時候需要的時間比較久，為了節(jié)約時間，將短程力計算單元設計成流水的，也就是說第一組粒子對進入計算單元計算后，第二組能馬上跟上。在本文系統(tǒng)中，短程力計算單元只有一個輸入口。由于匹配單元匹配成功的粒子對個數是大于或等于1的，所以，需要一個仲裁器對匹配單元輸出粒子對進行仲裁，使匹配成功的粒子對流水輸出到短程力計算單元。

匹配單元的設計，重點在于粒子對的生成、過濾和仲裁。

2.1 粒子對生成

為了使進入匹配單元的粒子能夠高效有序配對，本文系統(tǒng)的匹配單元設計了粒子對生成，如圖2所示粒子對生成主要包括FIFO和流寄存器部分。匹配單元每輪將有8組粒子對進行匹配。當本輪匹配結束后，就需要新一輪的粒子進入到匹配單元，此時，需要通過FIFO使暫存的粒子進入到匹配單元。此外，進入匹配單元的粒子將開始進行匹配，但數據預處理可能依舊有粒子輸出，因此，需要FIFO將它們暫存，等下一輪再輸出。FIFO中輸出的粒子在進入到流寄存器中時，為使靜態(tài)粒子和動態(tài)粒子進行高效合理地配對，同時滿足時序要求，本文使用流寄存器，將靜態(tài)粒子和動態(tài)粒子以粒子對的形式暫存起來。

2.2 粒子對過濾

根據本文第1節(jié)，實際模擬的系統(tǒng)相當于將一個無限的空間在三維坐標軸上平鋪為大小相同的盒子，圖3是粒子P的模擬空間，其中以粒子P所在的盒子為主盒，主盒的鄰域以非灰色表示，盒子邊緣大小為截止半徑rc，則三維坐標系中可能與粒子P進行匹配的其他粒子在相鄰的27個盒子內。

圖3 粒子P的模擬空間

平均配對成功率見式(8):

(8)

從式(8)可以看出，平均配對成功率不高，只有15.5%。為了盡可能提高短程力計算的利用率(盡量達到不停計算，而不產生空轉)，使每次短程力計算至少有一個有效輸入就需要7個以上距離計算器，理論上這就可以達到使短程力計算流水線充分運轉。所以，本文系統(tǒng)的匹配單元使用了8個距離計算器。

如圖2所示，粒子對過濾主要包括距離計算器和向量寄存器部分。當流寄存器中的粒子對在進入到8個距離計算器后，將根據第1節(jié)的匹配原則對粒子對進行匹配。

圖4的狀態(tài)轉移圖詳細介紹了距離計算器的工作過程。當距離計算器收到粒子對信息后，首先對盒子號和盒內編號進行比較，其次計算粒子間距離，最后輸出匹配成功的粒子對。如果靜態(tài)粒子和動態(tài)粒子不在同一個盒子內，就根據式(1)～(7)對粒子間距離進行計算。如果靜態(tài)粒子和動態(tài)粒子在同一個盒子內，比較靜態(tài)粒子和動態(tài)粒子的盒內編號，如果動態(tài)粒子盒內編號大于靜態(tài)粒子盒內編號，則進行粒子間距離的計算，反之則不需要進行粒子間距離的計算。如果粒子間距離都滿足式(1)～(7)，距離計算器的匹配結果為1，反之為0。

圖4 距離計算器狀態(tài)圖

2.3 粒子對仲裁

由于本文系統(tǒng)的短程力計算單元只有一個輸入端口，并且為了排除空粒子，本文系統(tǒng)的匹配單元設計了粒子對仲裁，粒子對仲裁主要包括與運算和仲裁器部分，如圖2所示。當距離計算器的匹配結果輸出時，由于存在空粒子，故還需要對圖2中暫存在向量寄存器中的8組距離計算器匹配結果和流寄存器中的空粒子標志位進行與運算，避免空粒子進入短程力計算單元。

由于8組與運算結果同時輸出，并且8組結果中含有無效粒子對，此時，需要流水輸出有效粒子對，故本文系統(tǒng)的匹配單元使用了仲裁器，即可以判斷粒子對的有效性又可以流水輸出。同時考慮到仲裁器有多個輸入端口和一個輸出端口，并且每個有效的粒子對請求輸出的機會是平等的，在這種情況下，存在粒子對輸出先后問題。為了在保證公平的基礎上解決這個問題，采用了仲裁器的輪詢仲裁策略。通過輪詢仲裁依次輸出多個匹配成功的粒子對，為短程力計算單元實現流水供數。下面對輪詢仲裁進行詳細介紹。

輪詢仲裁策略[22]的基本思路是，當請求者發(fā)出一個請求并獲得許可后，它的優(yōu)先級在接下來的仲裁中就變成了最低。簡而言之，就是每個請求者的優(yōu)先級不是固定不變的，而是在優(yōu)先級達到最高之后(即獲得許可后)變成最低，并且還會根據其他請求者獲得許可的情況進行相應的調整。這樣當有多個請求者時，每個請求者將依次獲得許可，即便之前高優(yōu)先級的請求者重新發(fā)出新的請求，也會等前面的請求者都獲得許可后再輪到它。關于輪詢仲裁策略的硬件實現，就必須提到固定優(yōu)先級仲裁策略，固定優(yōu)先級，顧名思義，就是每個模塊的優(yōu)先級是固定的，提前分配好的，如果有兩個模塊同時發(fā)出請求，則優(yōu)先級高的模塊可以獲得許可。固定優(yōu)先級的硬件設計是通過輸入一個多位的請求信號，其中每一位代表一個模塊的請求，輸出也是多位的許可信號，每一位代表對應模塊的許可信號。默認最低位優(yōu)先級最高，最高位優(yōu)先級最低。從而實現固定優(yōu)先級仲裁。輪詢仲裁的硬件設計是在固定優(yōu)先級仲裁的基礎上進行的，輪詢仲裁硬件設計需要兩個并行的固定優(yōu)先級仲裁，還需要一個屏蔽信號。屏蔽信號原理是將一個多位的請求信號之前獲得許可的位以及之前的位都屏蔽掉，允許通過的是之前優(yōu)先級更低的位。如果那些位上之前有發(fā)出請求但沒獲得許可，在經過屏蔽信號后就可以輪到它們。每次獲得許可之后，屏蔽信號里0的位數會變多，從而屏蔽掉更多位，直到所有低優(yōu)先級發(fā)出的請求均都獲得過一次許可后，請求信號和屏蔽信號邏輯與的結果變成全0，此時表明輪詢結束，需要進行新一輪的輪詢。圖5是輪詢仲裁的硬件設計圖，屏蔽信號和請求信號在邏輯與的結果不全為0時，即存在沒被屏蔽掉的請求，此時優(yōu)先獲得許可。當屏蔽信號和請求信號邏輯與結果全為0時，即沒有請求需要屏蔽。此時需要將邏輯與的結果和請求信號進行邏輯與。最后將邏輯與的結果和之前的結果進行邏輯或，最后輸出的即獲得許可。

圖5 輪詢仲裁硬件設計圖

3 軟件設計

系統(tǒng)實現過程中，一方面需要不斷將硬件實現的結果和軟件的結果進行對比，以確保硬件設計實現的正確性。另一方面，需要不斷調整優(yōu)化算法，實現對硬件資源的高效利用。系統(tǒng)的軟件設計是基于C語言實現，系統(tǒng)的開發(fā)流程包括以下幾點：

1)獲取軟件代碼和輸入數據，生成最終的輸出數據(標準結果，作為參照)；

2)根據硬件頂層設計，重新劃分模塊，每個模塊都有軟硬件兩種實現；

3)運行軟件模塊獲得軟件結果文件，軟件結果文件經過比對后，可以作為回歸測試的標準結果文件使用；

4)運行硬件模塊獲得硬件結果文件；

5)對軟硬件結果文件進行比對。

軟件主要使用浮點數據作為輸入輸出，同時軟硬件輸入數據，均由軟件統(tǒng)一提供。匹配單元軟件設計的主要思路是判斷兩個粒子之間的距離是否在截止半徑范圍內。根據第1節(jié)引入的偏序法和平面法，實現對粒子對的篩選，將不滿足短程力計算的粒子對過濾，其軟件流程如圖6所示。

圖6 軟件流程圖

如圖6所示，首先通過偏序法判斷靜態(tài)粒子和動態(tài)粒子的盒子號和盒內編號來排除粒子對重復計算的可能性。如果靜態(tài)粒子和動態(tài)粒子在同一盒子內，并且靜態(tài)粒子的盒內編號大于等于動態(tài)粒子的盒內編號，則直接排除掉該粒子對。當且僅當靜態(tài)粒子和動態(tài)粒子不在同一盒子內或者靜態(tài)粒子的盒內編號小于動態(tài)粒子的盒內編號，才可以進一步進行篩選。下一步通過平面法將滿足偏序法的粒子對進行距離比較，根據線性放縮法，將兩點距離公式從高階降為低階，使粒子到球面的距離簡化為到平面的距離。在對粒子進行距離比較時，根據式(1)～(7)，將公式左邊數據進行運算處理，此過程中不涉及高階運算，全部運算都停留在低階，通過使用c語言庫中自帶函數，使整個運算過程簡單易處理，相比于硬件這部分運算內容更易實現，運算結束再進行距離比較得到同時滿足偏序法和平面法的粒子對組合。最后，匹配成功的粒子對即滿足短程力計算的粒子對。據此，實現對粒子對的過濾。

4 實驗結果與分析

4.1 軟件測試

測試使用的粒子數據集(7 051個有效粒子數據)是來自蛋白質晶體結構資料數據庫(PDB, protein data bank)的多聚泛素酶基因(UBQ)數據集(https://www.rcsb.org/structure/1UBQ)。由于數據過多只說明部分數據，表1列出了該數據集中的部分粒子數據。輸入的數據主要是0號盒子和5號盒子內的粒子。前期工作準備完畢后，就根據第3節(jié)的匹配原則對粒子進行匹配。其中截止半徑為12?。通過軟件編譯得到的匹配結果如表2所示。

表1 粒子數據

表2 匹配結果

4.2 硬件測試

本文使用SpinalHDL實現前述匹配單元的硬件設計，SpinalHDL是一種開源的高級硬件描述語言，是基于Scala的全新硬件設計語言，更準確的說是一個基于scala的HDL開發(fā)框架，不同于以前的HDL，它解決了Verilog的例化不方便，大量重復聲明，函數不能帶參，參數化笨拙，錯誤檢測弱等問題。并且和傳統(tǒng)的集成電路流和諧共存，適合大規(guī)模的系統(tǒng)級芯片(Soc，system on chip)的開發(fā)。在系統(tǒng)通過SpinalHDL實現硬件設計之后，還需對其進行板級測試驗證，驗證所設計的邏輯代碼是否符合預期的要求，測試結果是否與軟件仿真結果一致，是否滿足時序要求。在現在的集成電路的驗證中，主要有兩種方法，一種是隨機驗證，另外一種是形式驗證(formal verification)。隨機驗證就是通過大量的隨機向量來驗證測試是否符合要求(通常是將硬件的結果和軟件的結果進行比較，或者通過斷言比較)。形式驗證的主要思想是通過使用形式證明的方式來驗證一個設計的功能是否正確。形式驗證可以分為3大類：等價性檢查(equivalence checking)、模型檢查(model checking)和定理證明(theory prover)。本文所使用的方法是基于FPGA開發(fā)平臺的隨機驗證?；痉椒ㄊ牵涸贔PGA中調用生成的設計模塊(以IP形式)，添加一些外圍邏輯，來給此設計模塊提供測試輸入；通過主機發(fā)送控制命令，對設計模塊進行測試；回收測試結果，將測試結果與預期的軟件仿真結果比較。使用的FPGA開發(fā)平臺是Xilinx公司的VCU128開發(fā)套件。VCU128 開發(fā)板采用全新 Xilinx VU37P HBM FPGA，利用堆疊芯片互連將 8 GB HBM DRAM裸片添加到封裝基板上的 FPGA 裸片旁邊。支持高帶寬存儲器(HBM) 的 Xilinx FPGA 是計算帶寬問題(與在 PCB 上使用 DDR4 等并行內存相關)的最優(yōu)解決方案。VCU128提供PCIe Gen3 x16接口，插到一臺服務器主機上。主機平臺安裝Ubuntu18.04 Linux操作系統(tǒng)，在其上運行驅動程序和測試程序。測試框架如圖7所示。

圖7 驗證平臺框架示意圖

圖7中，測試腳本程序調用測試程序，測試程序使用驅動程序中的讀寫函數，跟VCU128開發(fā)平臺上的FPGA芯片通訊。FPGA芯片接收命令，輸入控制邏輯從存儲器指定位置讀取測試輸入數據，提供給待測模塊進行測試，測試結果通過測試輸出邏輯寫到指定存儲器位置。測試程序從指定存儲器讀取結果數據，跟期望值比較，輸出比較結果。

因數據量較大，本節(jié)主要說明基于部分數據的硬件測試結果。為了便于計算，將表1中的浮點數據都轉換成了16進制的定點數(22位整數，42位小數)。圖8是0號盒子內粒子的匹配結果，從圖的最左邊可以看到匹配單元的相關信息。其中重要信息有static_vec(靜態(tài)粒子向量)、is_real(空粒子標志)、dynamic(動態(tài)粒子)、number(盒子號)、idx(盒內編號)、loc_x(粒子x軸坐標)、loc_y(粒子y軸坐標)、loc_z(粒子z軸坐標)、match_cal_io_count(距離計算器匹配結果)、pairs_out(最后輸出粒子對)、pairs_out_valid(輸出粒子對有效標志)。其中靜態(tài)粒子和動態(tài)粒子都在0號盒子中，且動態(tài)粒子盒內編號是1。因為靜態(tài)粒子個數比較多，所以只展示盒內編號為0的粒子。從圖中可以看到距離計算器只有一對粒子匹配成功(0表示匹配不成功，1表示匹配成功)，最后經過仲裁器輸出的粒子對也僅有一對。

圖8 0號盒子粒子匹配仿真波形

圖9是8組粒子對都匹配成功情況。可以看出靜態(tài)粒子在0號盒子中，動態(tài)粒子在5號盒子中，由于共有8個靜態(tài)粒子，導致粒子數據信息比較多，故只展示了第1個靜態(tài)粒子的數據信息，經過距離計算器進行匹配，發(fā)現距離計算器的結果均為1，即所有粒子對都匹配成功，最后匹配成功的粒子對全部通過仲裁器一一輸出。

圖9 不同盒子粒子匹配仿真波形

結果表明，圖8、圖9硬件測試結果與表2的理論匹配結果完全一致。

與此同時，為了說明偏序法和平面法可以節(jié)省系統(tǒng)的資源消耗，對使用偏序法和平面法的情況以及只使用直接計算方法的情況分別進行了測試驗證。其中使用直接計算方法和使用偏序法、平面法的資源對比如表3?？梢钥闯?，使用偏序法、平面法比使用直接計算方法要多消耗7.6%的LUT(Look-Up-Table)，但卻節(jié)省了70%的DSP資源。相比于LUT資源，DSP資源對短程力計算系統(tǒng)更為珍貴。

表3 不同情況下粒子匹配資源消耗

5 結束語

本文針對提升分子動力學模擬中粒子間非鍵合力的計算效率的問題，提出了盒子劃分、靜態(tài)粒子、動態(tài)粒子和空粒子標志位等概念。通過將盒子內的靜態(tài)粒子和動態(tài)粒子進行配對生成粒子對，有效地提高了匹配效率。此外，還提出了粒子對篩選的兩種方法：偏序法和平面法，使用偏序法和平面法節(jié)省了系統(tǒng)70%的DSP資源。在分子動力學模擬中，短程力計算單元是計算量最大的部分，而匹配單元過濾掉不需要計算的粒子對減少了整個系統(tǒng)的計算量。從測試結果可知，匹配單元有效地過濾掉了粒子間距離較大的粒子對，輸出了與理論分析一致的結果。

目前本文系統(tǒng)的匹配單元最高運行頻率達到115 MHz，下一步可繼續(xù)優(yōu)化匹配單元，如減少組合邏輯內容，使資源消耗中LUT數量減少。LUT數量的減少有助于提升工作頻率，這對于進一步提升短程力計算系統(tǒng)整體性能具有重要意義。