基于gpu加速分子動力學模擬的長程受力計算部分

基于gpu加速分子動力學模擬的長程受力計算部分

0gpu長件模型的改進分子模擬是對微觀分子體系的有效研究工具，在化學、生物和物理方面得到了廣泛的應用。例如，生物分子的模擬已經(jīng)從不到500個原子系統(tǒng)的皮秒模擬（ps）發(fā)展到了數(shù)千萬原子系統(tǒng)的微秒模擬。然而，復雜生物分子系統(tǒng)中的大多數(shù)過程的典型時空-微秒（甚至秒）規(guī)?？梢赃_到非常小的微秒（甚至秒）規(guī)模。在這方面，當前的計算還不夠充分。通常所說的計算能力來自于中央處理器(CPU),而近年來圖形處理器(GPU(GraphicProcessingUnit圖形處理單元))在某種意義上顯示了更強大的計算能力.最新的單GPU峰值已超過1TFlops(單精度).而且隨著GPU程序開發(fā)平臺的日益進步,GPU的編程難度也大為降低.對于適合GPU硬件體系結構的算法,用GPU實現(xiàn)可以達到幾十倍甚至一百多倍的加速.分子動力學模擬中諸如蛙跳迭代、鄰居列表生成,近程作用(范德華力與靜電力的近程部分)計算等部分能在GPU上獲得很好的加速.同時,在納秒級的生物體系模擬中發(fā)現(xiàn),如果以截斷法計算靜電作用,即使采用全原子力場和顯式溶劑模型,計算與實驗結果間也有很大差別.但如果其他條件不變而嚴格計算靜電作用,就能吻合較好.可見靜電作用的準確計算對模擬結果的真實性非常關鍵.生物體系的分子動力學模擬一般采用“周期性邊界條件”,靜電作用的計算本質上是“周期性邊界條件”下的N-body問題,其直接計算的復雜性為O(N2).對于較大規(guī)模的體系,這種方法的計算量非常大,因此,雖然已有GPU加速的N-body問題程序,但尚難以滿足生物體系的模擬需求.而對靜電作用的長程部分的優(yōu)化算法包括Ewald,PPPM(Particle-Particle-ParticleMesh),PME(Particle-MeshEwald),MSM(multilevelsummationmethod)等.Ewald方法最初用于周期性的晶體體系中,仍是O(N2)或O(N2/3)的算法,不適合大體系的計算.PPPM算法在某些情況下計算速度比PME快,但是計算精度較PME差,而且在計算長程力對體系維里的貢獻時只能夠得到維里張量的跡,不能準確得到維里張量的各個分量,不適合有壓力控制的體系的計算.MSN算法是一種比較新的算法,其GPU程序的效率較高,但還沒有長程力對體系維里等物理量貢獻的公式,應用上有困難.PME是目前較常用的一種O(NlogN)的算法,精度也較高,但某些環(huán)節(jié)不很適合GPU計算.為此,本文在CUDA(ComputeUnifiedDeviceArchitecture計算統(tǒng)一設備架構)GPU編程平臺上針對Gromacs-3.3.3(GMX)軟件包中的PME算法提出了一套優(yōu)化實現(xiàn),以盡量避免效率損失.1靜電作用的估計考慮一個帶有點電荷的N個原子的體系,各原子的位置為r1,r2,…,rN,電荷為q1,q2,···,qN,且q1+q2+···+qN=0.模擬區(qū)域用三個向量aα(α=1,2,3)表示(圖1).在周期性邊界條件下,所有原子產(chǎn)生的靜電勢能為Ewald將式(1)分解為兩個無條件收斂的級數(shù):近程部分(Edir),長程部分(Erec),再加上修正項(Ecorr).于是式(1)可表示為其中Edir與Ecorr的計算量較小,本文不予討論.是靜電作用中計算量最大的部分,也是PME算法的核心部分.將模擬空間劃分為傅立葉變換的網(wǎng)格,再引入B樣條插值,式(3)可以近似為式(3),(4),(5)中各變量的物理意義及數(shù)學表達式請參閱文.PME的精度與傅立葉變換網(wǎng)格的大小及B樣條插值的階數(shù)有關.在一定的傅立葉變換網(wǎng)格大小及4階或6階B樣條插值的情況下,PME能量計算的誤差為10-5~10-6量級,作用力計算的誤差為10-4~10-5,GPU單精度計算也能夠達到這樣的精度.2fft網(wǎng)格的離散GPU適于單指令多數(shù)據(jù)(SIMD)的算法,這類算法應滿足以下特點:(1)GPU內核函數(shù)各線程處理的數(shù)據(jù)間相關性小;(2)CPU主存與GPU顯存間數(shù)據(jù)交換少.更詳盡的介紹請參閱文.本章分析PME算法各步驟在GMX中的實現(xiàn),并設計其在GPU中的實現(xiàn).PME算法的實現(xiàn)即求解公式(4)與(5),其效率提高的一個重要因素為應用FFT代替對波矢的直接求和,即用B樣條插值將點電荷分布到FFT的網(wǎng)格中,然后在該網(wǎng)格上作三維正向傅立葉變換,再對變換后的網(wǎng)格直接加和以得到靜電作用長程部分的能量.如需得到靜電作用長程部分的受力,對正向FFT變換的結果處理后作三維逆向傅立葉變換,再對變換的結果插值可得到靜電力長程部分.在GMX中,可分為以下幾步:1)PME算法參數(shù)確定.首先根據(jù)模擬體系的性質與模擬精度的要求選取FFT網(wǎng)格的大小與PME插值的階數(shù).在采用GMX力場、4階插值及FFT網(wǎng)格為0.12nm時,PME算法得到的靜電勢能的精度為5×10-3,這個精度比范德華力度精度高,在某些模擬中可以適當增大FFT網(wǎng)格,從而減少計算量.這步是確定PME算法的參數(shù),計算量可忽略.2)點電荷在FFT網(wǎng)格的離散,即得到式(4)的Q數(shù)組.在生物大分子體系中,原子間的平均距離要比0.12nm大,因此在FFT網(wǎng)格為0.12nm的情況下,FFT網(wǎng)格的數(shù)目要比體系中的原子數(shù)N大,即一個粒子所帶的點電荷可離散到幾個FFT網(wǎng)格.需先根據(jù)各原子的實際坐標求得其比例坐標及各原子所對應的FFT網(wǎng)格,然后建立B樣條插值函數(shù)系數(shù)(M)矩陣,如需計算各原子受力,還需建立B樣條插值函數(shù)導數(shù)系數(shù)矩陣.此步是整個PME算法中花費時間最長的,可分為三步,分析如下:(1)原子的比例坐標與FFT網(wǎng)格索引可以看作原子自身的屬性,模擬體系中共有N個原子,CPU程序中,通過對N個原子的循環(huán)操作求得各個原子的比例坐標,再根據(jù)各個原子的比例坐標與模擬體系的FFT網(wǎng)格求得各原子的FFT網(wǎng)格索引,每個原子的循環(huán)中沒有與其他原子相關的量.這是典型的SIMD算法,適合GPU實現(xiàn).(2)B樣條插值函數(shù)的系數(shù)與B樣條插值函數(shù)導數(shù)的系數(shù)是與體系的原子數(shù)、各原子的比例坐標、B樣條插值階數(shù)相關的量.每個原子在x,y,z方向上各有插值階數(shù)個系數(shù),這樣在每個方向共有(原子個數(shù)×插值階數(shù))個系數(shù).且每個原子的系數(shù)與其他原子無關.與上一環(huán)節(jié)一樣,這步適合GPU計算.(3)將原子的電荷離散到FFT網(wǎng)格,這個環(huán)節(jié)可以抽象成N個原子對應Ngrid(FFT網(wǎng)格數(shù))個網(wǎng)格的問題.由于一個點電荷可離散到幾個網(wǎng)格,網(wǎng)格的值對應原子點電荷離散值的累加值.這種問題在CPU程序中很容易實現(xiàn),只需在原子的每步循環(huán)中用網(wǎng)格收集原子對網(wǎng)格的累加值即可,在GPU中卻困難得多.如果在該環(huán)節(jié)的GPU程序內核函數(shù)中,用一個線程處理一個原子,每個原子根據(jù)自己的FFT網(wǎng)格索引將電荷離散到相應網(wǎng)格,內核函數(shù)各個線程是并發(fā)的,但會出現(xiàn)多個原子對應一個網(wǎng)格的情況,即當一個線程從這個網(wǎng)格讀取數(shù)值,進行累加,再寫回網(wǎng)格時會發(fā)生“同時寫內存”的沖突,很大程度上限制了GPU效率的發(fā)揮.為了避免寫內存的沖突,同時還能用上GPU多線程的特性,將這個環(huán)節(jié)分解為兩步:ue591收集對FFT網(wǎng)格有貢獻的原子信息,這步需要申請額外的內存存儲對網(wǎng)格有貢獻的原子信息,該內存中有一個二維數(shù)組,其列數(shù)為FFT網(wǎng)格個數(shù),行數(shù)為對網(wǎng)格有貢獻的最大原子數(shù).這步內核函數(shù)發(fā)起的有效線程數(shù)為原子的個數(shù),每個線程處理一個原子,計算這個原子對哪些網(wǎng)格有貢獻,將原子信息存儲到二維數(shù)組的相應位置;ue592根據(jù)上一步結果與B樣條插值系數(shù)將各原子電荷離散到網(wǎng)格,這步內核函數(shù)發(fā)起的有效線程數(shù)為FFT網(wǎng)格數(shù),每個線程操作一個網(wǎng)格.經(jīng)過兩步分解,雖然用了更多的GPU內存,但可以應用GPU多線程的特性,計算速度有較大的提高.3)對FFT網(wǎng)格進行三維FFT變換.在GMX中此步采用FFTW-3.1.2軟件包.NVIDIA的CUDA編程平臺也提供了FFT函數(shù)CUFFT,本文的GPU程序直接采用之.4)靜電勢能的求解.在CPU程序中,由變換后的Q數(shù)組,直接在傅氏空間求和得到靜電勢能,并將變換后Q的數(shù)組與式(4)中的B數(shù)組、C數(shù)組相乘,為下一步三維FFT逆變換做準備.在GPU程序中,求Q數(shù)組與B數(shù)組及C數(shù)組的乘積很容易用上GPU多線程的特性,而能量求和卻難以在同一個內核函數(shù)內實現(xiàn).為此,需要申請一個一維數(shù)組,用來存放每個FFT網(wǎng)格對能量的貢獻,然后通過兩個內核函數(shù),先分塊歸約該數(shù)組,再得到整個數(shù)組的和.5)長程靜電力的求解.這步可分解為兩步:(1)CPU程序中,首先對上一步中得到的Q與B,C數(shù)組的乘積進行三維FFT逆變換,得到卷積θrec*Q.GPU程序中,與步驟3)一樣,先調用CUFFT中的傅立葉逆變換函數(shù)得到θrec*Q.(2)CPU程序中,對原子循環(huán),將各原子對應的θrec*Q與ue014Q/ue014rαi(在第2)步求B樣條函數(shù)導數(shù)的系數(shù)時構建)相乘得到各原子所受長程靜電力ue014珘Erec/ue014rαi.GPU程序中,用一個線程處理一個原子所對應的θrec*Q與ue014Q/ue014rαi,得到各原子的受力.該步GPU內核函數(shù)發(fā)起的有效線程數(shù)為原子的總數(shù),效率較高.3運行gpu和cpu的性能3.1gpu加速案例為了比較GPU程序與CPU程序的性能,選擇了7個生物分子模擬的體系,見表1.測試平臺為惠普工作站,Intel(R)Xeon(R)E54302.66GHzCPU,8條667MHz(1.5ns)2GBDDR2內存.GPU為GeForceGTX2951.3版,主頻1.242GHz.Linux版本為2.6.18-128.el5,CPU程序以gcc-4.1.2編譯.GPU程序開發(fā)編譯環(huán)境為NVIDIA-Linux-x86_64-180.29-pkg2.run驅動,cuda-sdk-linux-2.10.1215.2015-3233425.runSDK,cudatoolkit_2.1_linux64_rhel5.2.run的CUDAToolkit.分別對PME的CPU與GPU程序進行測試,各步驟的運行時間如表2與3所示.GPU程序相對CPU程序各步驟的加速比如表4.表2與3中的每個體系運行10次,所列為平均時間.表中各步驟的時間為實際消耗值,而總時間還包括申請與清空內存等的時間,故各步驟時間的和與總時間有偏差.3.2gpucs-gs軟件加速比分析根據(jù)PMEGPU程序各步驟的加速比可以將各步驟分為三種類型,如圖2(a)~(c)所示.1)適合GPU計算模式的算法.這類算法的GPU/CPU加速比隨計算量的增大而提高,當計算量達到GPU硬件的飽和值時(每個線程處理一個原子),加速比會達到平臺,如圖2(a)中描述的步驟2)(2)與4).而圖2(a)中描述的步驟2)(1)與5)(2)則沒有.2)不很適合GPU計算的算法.主要有兩類原因抑制了GPU性能的發(fā)揮:(1)線程間數(shù)據(jù)的相關性,可能會導致線程間的同步;(2)訪問GPU顯存的量,GPU各線程訪問其全局內存需要的時間比浮點計算的時間長一到兩個量級,如果算法需要訪問內存的操作比浮點計算的操作多,會影響GPU性能的發(fā)揮.圖2(b)中的步驟屬于前者,共有8次對GPU全局內存的讀寫,而浮點計算只有18個操作,如果按讀寫一次內存為100個時鐘周期,一次浮點操作為1個時鐘周期計算,讀寫內存的時間是浮點計算的44倍多.3)CUDA編程環(huán)境自帶函數(shù)的應用.目前它們包括FFT與BLAS數(shù)學庫,分別稱為CUFFT,CUBLAS.編程人員可以很方便地應用這兩個庫中的函數(shù)接口開發(fā)自己的程序.但從圖2(c)看,CUFFT對正向與反向FFT的加速比在1~7倍之間,在PMEGPU程序各個步驟中的加速比是比較低的,而且CUFFT函數(shù)的性能不太穩(wěn)定.但由于沒有CUFFT庫的源碼,無法分析導致波動的原因.從表2,3中可以看出,圖2(b)中的步驟占用總時間的50%多,故此,描述整個程序的加速比圖2(d)與圖2(b)的趨勢很相近.而其中的步驟2)(3),計算量大且不太適合GPU計算,是程序整體性能的主要限制.在Nvidia1.3版的GPU硬件體系結構中,可以應用原子函數(shù)(AtomicFunction),但目前只支持int型與unsignedint型數(shù)據(jù)的算術與邏輯運算(數(shù)據(jù)的交換可以應用于浮點型數(shù)據(jù)).如果原子函數(shù)能夠支持浮點數(shù)據(jù)類型,步驟2)(3)可以與步驟2)(2)合為一個內核函數(shù),而只消耗步驟2)(2)的時間,這樣整個程序的性能將進一步提升2倍左右.從圖2中可看到,加速比并不是隨著粒子數(shù)的增大而線性提高的,出現(xiàn)了跳躍和波動,這和CPU程序與GPU程序不同的計算模式有關.CPU程序是單指令單數(shù)據(jù)(SISD)的模式,用循環(huán)來實現(xiàn)對所有原子的操作.PME算法中的某些步驟不是對所有原子進行操作,而是只對帶有電荷的原子有操作,如步驟2)(1),2)(3),5)(2).在CPU程序中,先判斷當前循環(huán)的原子是否帶有電荷,如果不帶有電荷則直接跳到下一次循環(huán),對CPU程序這些步驟的計算量隨帶電荷原子數(shù)線性增加,而不是總原子數(shù).而GPU程序是SIMD模式的,內核函數(shù)同時發(fā)起多個線程,盡管每個線程也是先判斷其操作的原子是否帶