張寶花, 徐 順

(中國科學(xué)院計算機(jī)網(wǎng)絡(luò)信息中心 超級計算中心, 北京 100190)

GROMACS軟件并行計算性能分析①

張寶花, 徐 順

(中國科學(xué)院計算機(jī)網(wǎng)絡(luò)信息中心 超級計算中心, 北京 100190)

分子動力學(xué)模擬是對微觀分子原子體系在時間與空間上的運(yùn)動模擬, 是從微觀本質(zhì)上認(rèn)識體系宏觀性質(zhì)的有力方法. 針對如何提升分子動力學(xué)并行模擬性能的問題, 本文以著名軟件GROMACS為例, 分析其在分子動力學(xué)模擬并行計算方面的實(shí)現(xiàn)策略, 結(jié)合分子動力學(xué)模擬關(guān)鍵原理與測試實(shí)例, 提出MPI+OpenMP并行環(huán)境下計算性能的優(yōu)化策略, 為并行計算環(huán)境下實(shí)現(xiàn)分子動力學(xué)模擬的最優(yōu)化計算性能提供理論和實(shí)踐參考. 對GPU異構(gòu)并行環(huán)境下如何進(jìn)行MPI、OpenMP、GPU搭配選擇以達(dá)到性能最優(yōu), 本文亦給出了一定的理論和實(shí)例參考.

分子動力學(xué)模擬; GROMACS; 并行計算

2013年的諾貝爾化學(xué)獎授予了Martin Karplus、Michael Levitt和Arieh Warshel三位理論化學(xué)家, 以肯定他們在 “開發(fā)多尺度復(fù)雜化學(xué)系統(tǒng)模型”方面所做的貢獻(xiàn). 他們發(fā)展了大量分子動力學(xué)(MD: Molecular Dynamic)模擬基本理論和計算方法, 實(shí)現(xiàn)了更好地通過微觀相互作用研究體系的宏觀性質(zhì). 這一類強(qiáng)有力的計算化學(xué)方法滲透于迅速發(fā)展的分子生物學(xué)、生物信息學(xué)和生物物理學(xué)等前沿交叉學(xué)科, 同時逐漸形成一系列典型計算與模擬軟件, 包括: GROMACS、NAMD、LAMMPS、AMBER、CHARMM等[1-5]. 高性能并行計算技術(shù)應(yīng)用于分子動力學(xué)模擬, 使得體系模擬能夠在更大空間尺度和時間尺度上進(jìn)行. 比如GROMACS軟件的基礎(chǔ)框架設(shè)計就立足在高性能計算環(huán)境.

GROMACS是開源分子動力學(xué)通用軟件包, 面向百萬級甚至億萬級粒子體系的模擬[1]. 它具有多種并行模式: MPI、MPI+OpenMP、Thread-MPI, 還支持GPU等加速器的應(yīng)用. GROMACS 4.6以后針對GPU kernel進(jìn)行了重新設(shè)計, CPU和GPU的協(xié)調(diào)工作方式對計算性能變得更為敏感[6]. 針對MD模擬中的力場計算密集型任務(wù), 它采用區(qū)域分解技術(shù)和動態(tài)負(fù)載均衡技術(shù),并且計算長程作用力的任務(wù)(PME算法)能分類出獨(dú)立進(jìn)程負(fù)責(zé), 來加速優(yōu)化計算; 對于大規(guī)模粒子體系在并行計算時實(shí)行主區(qū)域分解技術(shù)(DD: Domain Decomposition)[7], Domain區(qū)域內(nèi)的計算任務(wù)會映射到一個PP rank (Particle-Particle (PP) interaction MPI rank)上面. 在并行規(guī)模較大時, 還會將PME計算分配獨(dú)立的ranks以減少rank間的通訊量來提高并行計算效率[8]. GROMACS雖然能自動調(diào)節(jié)rank間的計算負(fù)載, 但要達(dá)到最佳計算性能還很困難.

GROMACS的計算性能除了硬件影響外, 還有安裝編譯方式和并行計算策略選擇方式的影響, 在硬性約束固定和軟件安裝編譯優(yōu)化之后, 此時并行計算策略的選擇對運(yùn)行速度影響最為顯著. 并行策略方面可調(diào)節(jié)的方式有: 進(jìn)程并行數(shù)、節(jié)點(diǎn)劃分、區(qū)域分解、作用力計算相關(guān)的PME/PP任務(wù)分配和GPU/CPU協(xié)調(diào)工作等. 這些是提升分子動力學(xué)模擬并行計算性能的關(guān)鍵因素, 對于GROMACS程序的終端用戶意義重大.準(zhǔn)確高效地使用軟件, 不僅能夠節(jié)省時間, 更快獲取模擬計算結(jié)果加快科研進(jìn)度, 而且由于硬件設(shè)備的高效使用能夠大幅節(jié)省計算成本. 但現(xiàn)實(shí)是: 軟件程序?qū)τ诖蠖嘤脩魜碚f仍然是個黑盒子, 不可避免地存在使用效率上的陷阱. 比如選擇何種并行模式? 如何合理分配MPI+OpenMP? 如何有效使用GPU加速卡能夠達(dá)到計算絕對性能最優(yōu), 同時不會造成硬件資源的隱形浪費(fèi)? 這些都是用戶非常關(guān)注且迷惑的問題. 因此有必要對GROMACS的并行計算性能進(jìn)行分析研究.

本文分析GROMACS軟件在分子動力學(xué)模擬并行計算方面的實(shí)現(xiàn)策略, 先闡述MD模擬中涉及計算方法的并行化設(shè)計與實(shí)現(xiàn)關(guān)鍵問題, 在原理和算法分析之后, 以最優(yōu)并行計算為目標(biāo), 提出CPU并行中如何調(diào)節(jié)并行計算優(yōu)化策略方案. 文中還討論了異構(gòu)并行環(huán)境下GROMACS并行計算的實(shí)現(xiàn)原理, 并對GPU異構(gòu)加速時計算性能調(diào)節(jié)給出了有益參考.

1 GROMACS并行計算優(yōu)化策略

GROMACS采用了MD模擬中的區(qū)域分解并行設(shè)計方式, 以MPI為主要的并行環(huán)境, 加入OpenMP和CUDA-GPU實(shí)現(xiàn)熱點(diǎn)并行加速計算. 為獲得最佳計算性能, 常需要綜合考慮不同的并行組合方式. 根據(jù)GROMACS軟件并行設(shè)計方式以及影響并行計算性能的分子作用力計算的問題, 針對 MPI搭配OpenMP線程并行方式我們提出相應(yīng)的優(yōu)化策略方案, 并對存在GPU異構(gòu)加速情況下并行優(yōu)化作探索.

1.1 非成鍵力計算與PME計算性能優(yōu)化

優(yōu)化GROMACS MD模擬關(guān)鍵在于優(yōu)化分子作用力的計算. 分子作用力大致分為: 成鍵力和非成鍵力,非成鍵力有范德華力和靜電場力等形式. 直接計算非成鍵力的復(fù)雜度為O(N2), N為體系的粒子數(shù), 這也使得它成為并行優(yōu)化的關(guān)鍵問題. GROMACS采用PME(Particle Mesh Ewald, 包括LJ-PME)[9]算法來加速非成鍵力的計算, Ewald算法是PME的基礎(chǔ), 在粒子數(shù)為N的周期性邊界體系中, 靜電勢能U可以表達(dá)為

其中靜電常量k是關(guān)于介電常數(shù)ε的量k=1/4πε, 而rij,n=|ri-rj+n|表示電荷粒子i和j在跨越鏡像格矢n間的真實(shí)距離, 鏡像格矢n=(nx,nv,nz), 添加*號表示當(dāng)在原像中即n=(0,0,0)時排除i=j的情況. Ewald方法將慢收斂的式(1)轉(zhuǎn)換為兩個快收斂項和一個常數(shù)項

其中

實(shí)空間(direct space)表達(dá)式Udir中erfc為互補(bǔ)誤差函數(shù), 倒易空間(reciprocal space)表達(dá)式Urec中V為模擬盒子的體積, 系數(shù)β同時出現(xiàn)在三項中用于調(diào)節(jié)它們計算的相對權(quán)重. 比如將β系數(shù)取較大值,Udir會相對于r收斂很快, 因此可以采用截斷r來減少其計算量, 而不影響計算精度. 實(shí)際上Ewald方法的計算復(fù)雜度仍然為O(N2), 而PME在計算Urec時引入了快速傅里葉變換操作加快了計算, 同時采用大β系數(shù)加快Udir收斂便于計算中使用截斷操作, 最終PME總體上計算復(fù)雜度達(dá)到了O(Nlog(N)), 在大體系模擬中加速性能顯著.

非成鍵力按照作用范圍分為短程和長程非成鍵力.短程非成鍵力(如靜電場計算Udir)在實(shí)空間計算, 長程非成鍵力(如靜電場計算Urec)常變換到倒易空間計算. PME算法是用于加速長程非成鍵力計算的一種常用方法, 其它類似算法還有Reaction Field方法等.

圖1 GROMACS中分子間作用力的分類

GROMACS中計算作用力主要分為兩大類(見圖1): PP rank, 包括成鍵相互作用和短程非鍵相互作用(對應(yīng)于公式(2)中的Udir和U0), PP作用力計算的距離范圍閾值為rc(=max(rvdw,rcoulomb,rlist))[10]; PME rank, 包括長程非鍵相互作用(PME算法, 包括LJ-PME, 對應(yīng)于公式2中的Urec). PME rank通過從PP實(shí)空間處理器獲取電荷、坐標(biāo)等信息, 在倒易空間中進(jìn)行柵格化和3D FFT及反演計算, 最后將計算得到的力、能量及維里等返回實(shí)空間處理[11]. PME和PP rank的任務(wù)多寡與負(fù)載均衡在并行模擬優(yōu)化時需著重考慮.

GROMACS使用純CPU計算時, PP和PME的計算可以在同一個rank上, 也可以將PME計算部分開辟一些獨(dú)立的ranks處理. 前一種方式由于PME計算需要大量的全局通訊, 隨著并行規(guī)模的增加, 通訊復(fù)雜度會大大增加, 因此只針對小并行規(guī)模模擬. 后一種方式可以降低通訊復(fù)雜度, 故適合大規(guī)模并行模擬.例如使用64個ranks時, 其中分配16個PME ranks, 其余48個為PP ranks, PP和PME ranks僅在需要的時候通訊, 大大降低了通訊的復(fù)雜度(PME ranks減少4倍,通訊量就會減少16倍), 尤其是在大規(guī)模并行時, 這種方法往往能大幅提高計算效率. GROMACS在12 ranks以上并行時, 會自動進(jìn)行PME 獨(dú)立進(jìn)程的劃分. 在實(shí)際模擬計算中, 還可以使用tune_pme (GROMACS 5.0版本以上)小工具來測算最佳的PME ranks和3D-FFT倒易空間三維網(wǎng)格格點(diǎn)數(shù)[12].

GROMACS程序?qū)崿F(xiàn)了動態(tài)負(fù)載均衡, 通過在一定范圍內(nèi)調(diào)節(jié)截斷半徑以及PME格點(diǎn)數(shù)及格點(diǎn)間距離來調(diào)節(jié)成鍵和非成鍵相互作用的計算任務(wù)來實(shí)現(xiàn).但這種調(diào)節(jié)必須建立在對體系充分分析并設(shè)定合理參數(shù)的前提下和對硬件合理分配的基礎(chǔ)上進(jìn)行, 否則仍然會有大量負(fù)載不均衡的現(xiàn)象發(fā)生. GROMACS程序計算時, 由負(fù)載均衡導(dǎo)致的性能損失達(dá)到5%或以上時, 默認(rèn)會自動開啟動態(tài)負(fù)載均衡.

1.2 OpenMP線程計算加速

在CPU節(jié)點(diǎn)內(nèi)并行時, GROMACS可以使用兩種并行方式: 標(biāo)準(zhǔn)MPI方式和Thread-MPI并行方式. Thread-MPI并行方式是通過多線程機(jī)制模擬標(biāo)準(zhǔn)MPI方式實(shí)現(xiàn)的, 只限于節(jié)點(diǎn)內(nèi)應(yīng)用, 節(jié)點(diǎn)間并行模擬仍以標(biāo)準(zhǔn)MPI并行為主, 通過MPI通訊完成. 每一個MPI rank中, 可以開啟多個OpenMP線程加速. 需要注意的是: 在一定的可用核數(shù)時, 開啟OpenMP多線程并行, 由于硬件型號配置等原因, 計算速度不一定會提高[13].

1.3 GPU異構(gòu)計算加速

GROMACS 使用GPU加速時(目前只支持Verlet方法計算)[14], 如圖2[6]所示, 在每一步的計算中, 短程非鍵相互作用計算會加載到一個或多個GPU卡上. CPU同時在處理成鍵相互作用和長程非鍵相互作用(PME, 包括LJ-PME), 這部分計算可以使用OpenMP線程加速計算[15]. 在每步模擬中, GPU端接收CPU端傳輸?shù)泥徑恿斜怼⒆鴺?biāo)和電荷等信息. GPU在計算短程非鍵相互作用時, 當(dāng)有多個domain會分為local和nonlocal兩個stream處理, 得到短程非鍵計算力后再返回CPU端進(jìn)行其他計算處理, 每隔10-50步由CPU端進(jìn)行鄰接列表的更新. 因此, CPU和GPU之間的負(fù)載均衡對程序性能有很大的影響. 如果CPU和GPU任務(wù)分配不均, 就會造成一方超負(fù)荷運(yùn)轉(zhuǎn)狀態(tài), 而另一方一直在空閑等待狀態(tài), 造成硬件資源的隱形消耗.而影響這種負(fù)載均衡的因素有三: 一是CPU和GPU硬件的絕對計算能力; 二是CPU和GPU的使用搭配,即CPU/GPU 比例; 三是計算體系的特征.

在實(shí)際計算中, GPU異構(gòu)最佳并行模式與不同的硬件特征和運(yùn)行體系相關(guān). 例如對于一個非成鍵相互作用力計算量很大的體系, 假如配備的GPU性能非常強(qiáng)的話, 可以使用較少的CPU配備較多的GPU來達(dá)到這一目的, 但如果其他的CPU核心不能被其它任務(wù)所合理利用的話, 這樣勢必會造成硬件資源的閑置與浪費(fèi)[16]. 因此, 在運(yùn)行一個較大較長時間任務(wù)的時候, 應(yīng)該先做一個小的測試來確定最佳的計算配置.

GROMACS的非成鍵相互作用力有三種并行計算方式, 如表1分別可以用-nb gpu/cpu/gpu_cpu選項指定,它們在計算非鍵相互作用力的方面各有側(cè)重. 使用哪一種方式計算取決于硬件特征和體系的計算特性.

圖2 GROMACS中GPU與CPU協(xié)調(diào)工作機(jī)制

表1 GROMACS計算相互作用力的三種并行計算類型

本文提出了CPU和GPU異構(gòu)環(huán)境下的并行計算優(yōu)化策略:

1) 使用CPU并行環(huán)境模擬時, 我們將考慮影響性能的兩類主要因素: PP rank間任務(wù)負(fù)載不均衡對性能的影響; PME獨(dú)立rank對性能的影響. 并針對具體體系測試結(jié)果, 總結(jié)性給出計算性能高低判據(jù)和并行優(yōu)化策略方案.

2) 使用GPU異構(gòu)并行環(huán)境模擬時, 我們將考慮影響性能的兩類主要因素: CPU/GPU配比, 通過測試找到最佳搭配比, 使資源利用和計算性能最大化; CPU端使用MPI進(jìn)程和OpenMP線程加速, 通過測試找出合適的進(jìn)程與線程分配, 為GPU異構(gòu)并行環(huán)境下計算性能最大化提供實(shí)踐經(jīng)驗(yàn)參考.

2 實(shí)驗(yàn)設(shè)計

2.1 運(yùn)行環(huán)境

在中國科學(xué)院“元”超級計算機(jī)上安裝編譯GROMACS 5.1單精度浮點(diǎn)MPI版本, 并支持GPU加速. 采用了Intel C/C++ 2013_sp1.0.080版本編譯器, Intel MPI 4.1.3并行庫(支持了OpenMP線程并行), 鏈接了MKL高性能并行函數(shù)庫和fftw3.3.4函數(shù)庫, 采用了CUDA 6.0.37編譯器, 并針對AVX_256指令集進(jìn)行了優(yōu)化.

2.2 模擬體系

我們選擇了兩個具有代表性的生物體系: 一個為均相蛋白體系3MGO(抗原表位多肽與人白細(xì)胞抗原復(fù)合物)[17], 另一個為膜蛋白體系DPPC(膽固醇-磷脂膜蛋白)[18], 兩個體系結(jié)構(gòu)見圖3所示, 為突出計算核心部分已隱藏水分子.

圖3 3MGO(上)和DPPC(下)體系結(jié)構(gòu)圖

表2顯示了3MGO和DPPC體系的模擬參數(shù). 為了消除初始態(tài)對計算速度的影響, 我們先對體系進(jìn)行一段MD平衡模擬, 選擇分子動力學(xué)統(tǒng)計平衡后的體系態(tài)作為我們算例測試的初始態(tài).

表2 3MGO和DPPC體系模擬參數(shù)表

當(dāng)使用Verlet截斷方法計算時, 可以充分利用SSE、AVX 或CUDA, 且PP和PME rank都可以使用OpenMP線程加速. 當(dāng)使用group截斷方法時, 只有PME rank可以使用OpenMP加速, 相對Verlet截斷方式并行計算程度低, 且在非鍵力計算時會出現(xiàn)能量漂移, 故這種情況本文暫不考慮[10].

2.3 并行計算方式

我們首先考慮純CPU并行計算, 先考慮節(jié)點(diǎn)內(nèi)MPI+OpenMP如何搭配絕對性能最好, 之后對跨節(jié)點(diǎn)并行進(jìn)行性能分析, 并給出性能影響因素和調(diào)優(yōu)方向.加入GPU加速卡后, 首先測試驗(yàn)證節(jié)點(diǎn)內(nèi)CPU/GPU最佳配比, 之后再多節(jié)點(diǎn)時考慮MPI+OpenMP+GPU異構(gòu)的并行性能情況.

在大規(guī)模并行計算時, 我們首先采用程序自動評估PME獨(dú)立進(jìn)程, 根據(jù)性能表現(xiàn)決定是否有必要手動調(diào)節(jié)PME獨(dú)立進(jìn)程. 在運(yùn)行時, 動態(tài)負(fù)載均衡默認(rèn)打開, 可以調(diào)節(jié)PP任務(wù)間和PP/PME任務(wù)負(fù)載均衡.

“元”超級計算機(jī)CPU節(jié)點(diǎn)為雙路10核CPU, 共20個邏輯核心, 本文不考慮超線程情況, 因?yàn)檫@意味著兩個線程要共享一個單浮點(diǎn)數(shù)單元, 較大規(guī)模并行時性能并沒有帶來提升. 為方便表述, 在下文中, 我們均用Nrank表示MPI rank數(shù)量, Nth表示開啟的OpenMP線程數(shù), Npme和Npp表示使用的PME和PP的rank數(shù).

在本測試中, 由于長程非鍵力計算量較大, 而2塊GPU的計算性能相對20 CPU強(qiáng), 當(dāng)并行選項為“gpu_cpu”的情況時, 一部分CPU會模擬GPU的計算行為, 計算結(jié)果不甚理想. 因此我們僅考慮并行選項“-nb”為CPU和GPU的情況.

我們關(guān)注計算性能ns/day的變化, 每組數(shù)據(jù)都平行進(jìn)行了至少3次測試取平均值, 以消除系統(tǒng)運(yùn)行中網(wǎng)絡(luò)、內(nèi)存及I/O等繁忙造成的不穩(wěn)定因素影響.

3 結(jié)果與討論

3.1 僅CPU并行方式

3MGO和DPPC體系節(jié)點(diǎn)內(nèi)并行性能測試表明,當(dāng)用滿一個節(jié)點(diǎn)的所有邏輯核心數(shù)20(即Nth×Nrank=20)時, 獲得的ns/day數(shù)值最高, 性能最好.因此, 在接下來的測試中選擇總的模擬核心數(shù)均為20n(n為節(jié)點(diǎn)數(shù)).

我們將從PP任務(wù)負(fù)載不均衡性和PP/PME任務(wù)比兩個方面來衡量計算性能是否優(yōu)良. PP任務(wù)負(fù)載不均衡率包括平均負(fù)載不均衡率(Avg. imb%)與每步力計算的不均衡率(Load imb%), 在模擬中會通過改變DD cell的尺寸來調(diào)節(jié)PP rank的負(fù)載, 還會受到體系各項異性的影響. 一般認(rèn)為這兩個值越小越好, 在本測試中認(rèn)為小于5%即為相對均衡(小規(guī)模并行時小于3%). PME/PP比值反映了PME和PP rank的任務(wù)分配是否均衡, 理想值為1, 即PME和PP rank的任務(wù)分配是完全均衡的. 在實(shí)際中, 由于網(wǎng)絡(luò)通訊等原因,實(shí)現(xiàn)PP-PME負(fù)載完全平衡幾乎是不可能的. 在本測試中認(rèn)為小于此值在0.8-1.2之間相對均衡.

在較小并行規(guī)模時, PP rank間任務(wù)負(fù)載不均衡是導(dǎo)致性能損失的主要原因.

如圖4, 當(dāng)計算節(jié)點(diǎn)為1時, 我們分別取Nrank為1、2、10、20, Nth為20/Nrank. 測試表明, 對于3MGO和DPPC體系, Avg. imb%與Load imb%均為1.2%以下,即每個rank的PP計算任務(wù)和每計算步所有力的計算都是相對均衡的. 由ns/day數(shù)值來看, 選取較大的Nrank數(shù), 更有利于絕對性能提升. 如3MGO體系: 當(dāng)Nrank為10和20時, ns/day分別為14.8和14.5. 而Nrank為1和2時, ns/day分別為13.8和13.5. DPPC體系結(jié)果相似.

當(dāng)并行規(guī)模擴(kuò)展到2個節(jié)點(diǎn)時, Nrank分別取2、4、20、40, 相應(yīng)的Nth為 20、10、2、1. 如圖5所示, 2個節(jié)點(diǎn)計算結(jié)果與1個節(jié)點(diǎn)相一致, 當(dāng)Nrank＞Nth時, 絕對性能更優(yōu).

圖4 3MGO和DPPC體系在1個節(jié)點(diǎn)時的計算結(jié)果

值得注意的是, 當(dāng)Nrank=4時, PP任務(wù)負(fù)載不均衡率導(dǎo)致的Avg. imb%較高, 由此造成了部分性能損失(DPPC體系為4, 3MGO為2.9), 這時如果我們調(diào)整Nrank與Nth的相對數(shù)值, PP任務(wù)負(fù)載不均衡現(xiàn)象就會得到改善, 計算性能將會得到提升. 如DPPC體系, 我們?nèi)rank為8, Nth為5時, Avg. imb%降低到了1.9%, 計算絕對性能提升了12%. 對于3MGO體系, 調(diào)節(jié)Nrank數(shù)為8后, 計算性能提升了6.4%.

圖5 3MGO和DPPC體系在2個節(jié)點(diǎn)時的計算結(jié)果

較大并行規(guī)模下, 除PP任務(wù)負(fù)載不均衡外, PP/PME的任務(wù)分配不均是導(dǎo)致性能損失的主要原因.

圖6 3MGO體系和DPPC體系在5個節(jié)點(diǎn)內(nèi)計算結(jié)果

如圖6所示, 當(dāng)并行規(guī)模達(dá)到5個節(jié)點(diǎn)時, Nrank分別取5、10、50和100, Nth分別為20、10、2和1, Nrank＞=20時有獨(dú)立PME rank分配. 在較大并行規(guī)模下, PME/PP任務(wù)分配不均是導(dǎo)致并行性能損失的主要原因. 例如3MGO體系, 當(dāng)Nrank=50時, GROMCAS核心程序mdrun自動估算的Npme=15, PME/PP比為1.37, 這說明PME的任務(wù)負(fù)載比PP rank多出37%, 絕對性能為48.3 ns/day. 當(dāng)調(diào)整Npme=20時, PME/PP比為1.06, PP和PME任務(wù)分配均衡, 且PP rank內(nèi)任務(wù)負(fù)載也相對均衡, 由此帶來絕對性能提升了19.3%. 當(dāng)Nrank=100時, mdrun估算的Npme=28, PME/PP為1.67, 據(jù)對性能為47.73 ns/day. 當(dāng)調(diào)整Npme=40時, PME/PP比為1.11, 由此帶來絕對性能提升了22.8%.

合適的PME rank數(shù)可以根據(jù)tune_pme小工具估算, 亦可根據(jù)前面小規(guī)模測試數(shù)據(jù)估算. 前一種方法掃描計算所有可能情況, 計算比較耗時. 本文使用后一種方法. 例如對于3MGO體系, 在單節(jié)點(diǎn)Nrank=20時, Npme=8, PME/PP=1.03, 任務(wù)分配比較均衡. 由此推算, Npme≈0.4Nrank, 當(dāng)Nrank為50時, PME獨(dú)立rank數(shù)選擇20較為合適, 當(dāng)Nrank為100時, PME獨(dú)立rank數(shù)選擇40較為合適, 在實(shí)際計算時, 還需要通過PME/PP比例確定調(diào)整的方向.

對于DPPC體系, 從圖6可以看出, 當(dāng)Nrank(Npme)=100 (20)時, PME/PP任務(wù)分配不均(PME/PP=0.73)且PP任務(wù)負(fù)載不均衡率很高(Avg. imb%達(dá)到了32.5%), 當(dāng)調(diào)整Npme=16時, PME/PP比值為1.03, PME和PP的任務(wù)分配已達(dá)到比較均衡, PP任務(wù)負(fù)載不均衡性亦得到了較大改觀(Avg. imb降低到了3.4%), 性能提升了9.8%. PP任務(wù)負(fù)載不均衡性高時,與小并行規(guī)模類似, 還可以調(diào)節(jié)Nrank和Nth的相對數(shù)值, 進(jìn)行優(yōu)化, 例如當(dāng)Nrank(Npme)=20(4)時, Avg. imb%降低到了3.8%, 絕對性能較Nrank(Npme)=100 (20)時候提升了20%.

另外, 負(fù)載不均衡情況都得到改善后, 與小規(guī)模并行情況相同, Nrank＞Nth時, 計算的絕對性能更優(yōu). 因此為獲取最佳絕對計算性能, 在實(shí)際計算中盡量選擇較少的線程和較多的進(jìn)程, 減少多線程開銷, 使Domain數(shù)據(jù)本地化, 減少Cache一致性開銷和socket間的通訊開銷.

綜上分析, 我們總結(jié)出了CPU并行時性能調(diào)優(yōu)的方向或判據(jù), 如圖7所示.

圖7 CPU并行時性能調(diào)優(yōu)判據(jù)/方向圖

①當(dāng)存在PME獨(dú)立rank, 且PME/PP比小于0.8或者大于1.2時, PME和PP任務(wù)分配不均衡較為明顯,需要調(diào)整Npme來調(diào)整PME和PP的任務(wù)分配 (Npme可通過小規(guī)模數(shù)據(jù)測算).

②當(dāng)沒有Npme, 或PME和PP任務(wù)分配達(dá)到均衡時, 若PP rank內(nèi)任務(wù)負(fù)載不均衡時(Avg. imb%＜5%)通過調(diào)節(jié)增大Nrank和減少Nth的相對數(shù)值來減弱或消除不均衡性現(xiàn)象.

3.2 CPU與GPU異構(gòu)并行方式

加入GPU卡后, 我們首先在單節(jié)點(diǎn)中測試CPU/GPU最佳配比(此處設(shè)定Nth=1, 先考察總核心數(shù)與GPU卡的最佳配比)

圖8展示了3MGO和DPPC單節(jié)點(diǎn)內(nèi) CPU/GPU配比測試結(jié)果, 其中紅色框使用了2個GPU, 綠色框使用了1個GPU, 結(jié)果表明在本測試硬件平臺上, 3MGO體系采用18 CPU+1 GPU計算的絕對性能優(yōu)于20 CPU+1 GPU, 我們推論本測試算例CPU/GPU的理論最佳配比約為18(對于DPPC體系, 最佳配比約為16), 限于節(jié)點(diǎn)CPU最大僅為20, 測試顯示用滿節(jié)點(diǎn)內(nèi)所有CPU和GPU卡獲得的絕對性能最佳(3MGO為29.46 ns/day, DPPC為28.88 ns/day), 從使用者角度來講, 可以在同樣時間內(nèi)模擬更多步數(shù). 從硬件有效利用來看, 添加第二個GPU, 對于DPPC體系更加有利(DPPC絕對性能提升了75.6%, 3MGO性能提升了28.3%), 可以綜合考慮性能提升與硬件耗電造成的損失比重來決定全部用滿CPU和GPU或者選擇18CPU+ 1 GPU配比. 在本文中, 我們從使用者角度出發(fā), 選擇用滿節(jié)點(diǎn)內(nèi)所有CPU和GPU卡. 在以下的計算中, 我們選擇的CPU/GPU均為20/2.

圖8 3MGO和DPPC單節(jié)點(diǎn)內(nèi) CPU/GPU配比測試

上面測試CPU端不考慮線程, 接下來我們加入OpenMP線程的影響, 加入GPU后, 由于異構(gòu)并行的影響, PP rank任務(wù)負(fù)載很容易出現(xiàn)較大不均衡現(xiàn)象,因此我們只關(guān)注每步力的計算的負(fù)載均衡和計算的絕對性能. 在下面的測試中Load imb%均在2.5%以下,因此不做單獨(dú)記錄分析.

圖9 3MGO和DPPC體系多節(jié)點(diǎn)多線程測試數(shù)據(jù)圖

圖9 展示了兩個體系多節(jié)點(diǎn)多線程測試結(jié)果, 圖中隱藏了Nrank值, Nrank=20n/nth, n為節(jié)點(diǎn)數(shù). 由圖知,隨著并行規(guī)模增大, 使用線程更有利于絕對性能的提升. 這是因?yàn)? 在較大并行度時, 使用進(jìn)程和線程結(jié)合方式可以減少區(qū)域分解的數(shù)量和MPI ranks, 即間接減少了CPU-GPU-CPU之間的通訊開銷, 因此計算性能優(yōu)越性會逐漸顯現(xiàn). 但開啟多線程會有一定開銷,特別使對cache一致性開銷和socket間的通訊開銷, 線程并非開的越多越好. 根據(jù)圖9測試結(jié)果, 大規(guī)模并行時選擇線程數(shù)為2或5較為合適.

圖10 3MGO和DPPC體系CPU和GPU計算性能對比圖

值得注意的是, 當(dāng)Nnode＞1時, 由于每一個主機(jī)CPU計算時間比GPU計算時間慢(GPU的負(fù)載較少),增加MPI計算節(jié)點(diǎn)時, 使用GPU卡計算時性能會快速提高, 但在計算節(jié)點(diǎn)持續(xù)增加后, CPU-GPU-CPU之間的通訊開銷也會增加, 因此如果計算節(jié)點(diǎn)繼續(xù)增加, 使用GPU的加速效果就會越來越不明顯(如3MGO體系從1節(jié)點(diǎn)到2節(jié)點(diǎn)并行效率近90%, 到10節(jié)點(diǎn)時效率下降到不足40%), 直至與CPU計算性能相近. 如圖10,加入GPU后, 相較于純CPU計算, 性能提升最高為2.4倍, 且隨著并行規(guī)模的增加, GPU的加速效果下降.

在本測試中, 我們僅考慮了‘-nb’為cpu/gpu兩種情況, 當(dāng)短程非鍵相互作用相對長程非鍵相互作用力很多時(例如coulombtype=Reaction-field), 選擇‘-nb’為gpu_cpu模式可能會更好的改進(jìn)性能, 這種情況我們會在以后的測試中進(jìn)行討論.

4 總結(jié)

本文對GROMACS軟件在單純CPU和CPU/GPU異構(gòu)并行計算環(huán)境下, 如何提升計算性能進(jìn)行了較為深入的分析與研究. 對于單純的CPU并行方式, 負(fù)載均衡水平是一個重要的考察指標(biāo). 在低并行度時, 這種負(fù)載不均衡主要是PP rank任務(wù)負(fù)載的不均衡; 高并行度時, PME獨(dú)立rank與PP rank一起作用, PP/PME的任務(wù)分配不均是導(dǎo)致性能損失的主要原因. 測試結(jié)果顯示, 在總并行數(shù)N一定情況下, 采用較多的進(jìn)程,較少的線程更有利于計算性能的提高. 使用CPU/GPU異構(gòu)并行計算時, 著重考慮CPU和GPU的相對運(yùn)算性能和配比, 需測試決定最佳的CPU/GPU配比. 在多個節(jié)點(diǎn)并行時候, 需要考慮CPU-GPU-CPU通訊開銷的影響以及線程開銷. 通過實(shí)例測試可得每個rank使用2-5個線程會提高絕對運(yùn)算性能.

我們從原理闡述到實(shí)例討論, 為GROMACS的終端用戶提供了性能優(yōu)化的方向和依據(jù), 幫助其更容易理解這個軟件黑盒子并行優(yōu)化機(jī)理, 乃至其它分子動力學(xué)模擬軟件.

1 Abraham MJ, Murtola T, Schulz R, Páll S, Smith JC, Hess B, Lindahl E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. Softwarex, 2015, 1-2: 19–25.

2 Brooks BR, Bruccoleri RE, et al. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. Journal of Computational Chemistry, 1983, 4(2): 187–217.

3 Nelson MT, Humphrey W, Gursoy A, Dalke A, Kalé LV, Skeel RD, Schulten K. NAMD: A parallel, object oriented molecular dynamics program. International Journal of High Performance Computing Applications, 1996, 10(4): 251–268.

4 Pearlman DA, Case DA, Caldwell JW, Ross WS, Iii TEC, Debolt S, Ferguson D, Seibel G, Kollman P. AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Computer Physics Communications, 1995, 91(1): 1–41.

5 Plimpton S. Fast Parallel algorithms for short-range molecular dynamics. J Comp Phys, 1995, 117, 1–19.

6 http://www.gromacs.org/.

7 Hess B, Kutzner C, Van DSD, Lindahl E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. Journal of Chemical Theory and Computation, 2008, 4(3): 435–47.

8 Gruber CC, Pleiss J. Systematic benchmarking of large molecular dynamics simulations employing GROMACS on massive multiprocessing facilities. Journal of Computational Chemistry, 2011, 32(4): 600–606

9 Darden T, York D, Pedersen L. Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems. The Journal of Chemical Physics, 1993, 98(12): 10089.

10 ftp://ftp.gromacs.org/pub/manual/manual-5.0.7.pdf.

11 Essmann U, Perera L, Berkowitz M L, Darden T, Lee H, Pedersen LG. A smooth particle mesh Ewald method. The Journal of Chemical Physics, 1995, 103(19): 8577–8593.

12 Abraham MJ, Gready JE. Optimization of parameters for molecular dynamics simulation using smooth particlemesh Ewald in GROMACS 4.5. Journal of Computational Chemistry, 2011, 32(9): 2031–2040.

13 http://www.gromacs.org/documentation/acceleration_and_ parallelization.

14 Grubmüller H, Heller H, Windemuth A, Schulten K. Generalized Verlet algorithm for efficient molecular dynamics simulations with long-range interactions. Molecular Simulation, 1991, 6(1-3): 121–142.

15 Pall S, Hess B. A flexible algorithm for calculating pair interactions on SIMD architectures. Computer Physics Communications, 2013, 184(12): 2641–2650.

16 Kutzner C, Páll S, Fechner M, Esztermann A, deGroot BL, Grubmüller H. Best bang for your buck: GPU nodes for GROMACS biomolecular simulations. Journal of Computational Chemistry, 2015, 36(26): 1990–2008.

17 Sun Y, Liu Q. Differential structural dynamics and antigenicity of two similar influenza H5N1 virus HA-specific HLA-A* 0201-restricted CLT epitopes. RSC Adv, 2014, 5(3): 2318–2327.

18 Chiu SW, Pandit SA, Scott H, Jakobsson E. An improved united atom force field for simulation of mixed lipid bilayers. The Journal of Physical Chemistry B, 2009, 113(9): 2748–2763.

Parallel Computing Performance Analysis on GROMACS Software

ZHANG Bao-Hua, XU Shun
(Supercomputing Center, Computer Network Information Center, CAS, Beijing 100190, China)

Molecule dynamic (MD) simulation is used for simulation of atomic systems. It is a powerful tool for understanding relationships between microcosmic nature and macroscopic properties. In view of the problem of how to improve the performance of parallel MD simulations, in this paper taking the famous GROMACS software as an example, we study its implementation strategy on parallel computing performance, combining with the key principles of MD and the results of test cases, and we propose optimization strategies of performance under MPI and OpenMP parallel environment. It can provide theoretical and practical reference for optimizing performance of MD software. In addition, we also study that in the GPU heterogeneous parallel environment how to tune MPI/OpenMP and GPU to obtain the best performance.

molecule dynamic simulation; GROMACS; parallel computing

中國科學(xué)院青年創(chuàng)新促進(jìn)會(2016156);“一三五”規(guī)劃重點(diǎn)培育方向?qū)ｍ楉椖?CNIC_PY_1404)

2016-03-24;收到修改稿時間:2016-04-27

10.15888/j.cnki.csa.005471