張寶花, 徐 順

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

GROMACS軟件并行計(jì)算性能分析①

張寶花, 徐 順

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

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

分子動(dòng)力學(xué)模擬; GROMACS; 并行計(jì)算

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

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

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

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

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

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

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

優(yōu)化GROMACS MD模擬關(guān)鍵在于優(yōu)化分子作用力的計(jì)算. 分子作用力大致分為: 成鍵力和非成鍵力,非成鍵力有范德華力和靜電場力等形式. 直接計(jì)算非成鍵力的復(fù)雜度為O(N2), N為體系的粒子數(shù), 這也使得它成為并行優(yōu)化的關(guān)鍵問題. GROMACS采用PME(Particle Mesh Ewald, 包括LJ-PME)[9]算法來加速非成鍵力的計(jì)算, 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), 添加*號(hào)表示當(dāng)在原像中即n=(0,0,0)時(shí)排除i=j的情況. Ewald方法將慢收斂的式(1)轉(zhuǎn)換為兩個(gè)快收斂項(xiàng)和一個(gè)常數(shù)項(xiàng)

其中

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

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

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

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

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

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

1.2 OpenMP線程計(jì)算加速

在CPU節(jié)點(diǎn)內(nèi)并行時(shí), 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通訊完成. 每一個(gè)MPI rank中, 可以開啟多個(gè)OpenMP線程加速. 需要注意的是: 在一定的可用核數(shù)時(shí), 開啟OpenMP多線程并行, 由于硬件型號(hào)配置等原因, 計(jì)算速度不一定會(huì)提高[13].

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

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

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

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

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

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

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

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

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

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

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

在中國科學(xué)院“元”超級(jí)計(jì)算機(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編譯器, 并針對(duì)AVX_256指令集進(jìn)行了優(yōu)化.

2.2 模擬體系

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

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

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

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

當(dāng)使用Verlet截?cái)喾椒ㄓ?jì)算時(shí), 可以充分利用SSE、AVX 或CUDA, 且PP和PME rank都可以使用OpenMP線程加速. 當(dāng)使用group截?cái)喾椒〞r(shí), 只有PME rank可以使用OpenMP加速, 相對(duì)Verlet截?cái)喾绞讲⑿杏?jì)算程度低, 且在非鍵力計(jì)算時(shí)會(huì)出現(xiàn)能量漂移, 故這種情況本文暫不考慮[10].

2.3 并行計(jì)算方式

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

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

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

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

我們關(guān)注計(jì)算性能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)用滿一個(gè)節(jié)點(diǎn)的所有邏輯核心數(shù)20(即Nth×Nrank=20)時(shí), 獲得的ns/day數(shù)值最高, 性能最好.因此, 在接下來的測試中選擇總的模擬核心數(shù)均為20n(n為節(jié)點(diǎn)數(shù)).

我們將從PP任務(wù)負(fù)載不均衡性和PP/PME任務(wù)比兩個(gè)方面來衡量計(jì)算性能是否優(yōu)良. PP任務(wù)負(fù)載不均衡率包括平均負(fù)載不均衡率(Avg. imb%)與每步力計(jì)算的不均衡率(Load imb%), 在模擬中會(huì)通過改變DD cell的尺寸來調(diào)節(jié)PP rank的負(fù)載, 還會(huì)受到體系各項(xiàng)異性的影響. 一般認(rèn)為這兩個(gè)值越小越好, 在本測試中認(rèn)為小于5%即為相對(duì)均衡(小規(guī)模并行時(shí)小于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之間相對(duì)均衡.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

②當(dāng)沒有Npme, 或PME和PP任務(wù)分配達(dá)到均衡時(shí), 若PP rank內(nèi)任務(wù)負(fù)載不均衡時(shí)(Avg. imb%＜5%)通過調(diào)節(jié)增大Nrank和減少Nth的相對(duì)數(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個(gè)GPU, 綠色框使用了1個(gè)GPU, 結(jié)果表明在本測試硬件平臺(tái)上, 3MGO體系采用18 CPU+1 GPU計(jì)算的絕對(duì)性能優(yōu)于20 CPU+1 GPU, 我們推論本測試算例CPU/GPU的理論最佳配比約為18(對(duì)于DPPC體系, 最佳配比約為16), 限于節(jié)點(diǎn)CPU最大僅為20, 測試顯示用滿節(jié)點(diǎn)內(nèi)所有CPU和GPU卡獲得的絕對(duì)性能最佳(3MGO為29.46 ns/day, DPPC為28.88 ns/day), 從使用者角度來講, 可以在同樣時(shí)間內(nèi)模擬更多步數(shù). 從硬件有效利用來看, 添加第二個(gè)GPU, 對(duì)于DPPC體系更加有利(DPPC絕對(duì)性能提升了75.6%, 3MGO性能提升了28.3%), 可以綜合考慮性能提升與硬件耗電造成的損失比重來決定全部用滿CPU和GPU或者選擇18CPU+ 1 GPU配比. 在本文中, 我們從使用者角度出發(fā), 選擇用滿節(jié)點(diǎn)內(nèi)所有CPU和GPU卡. 在以下的計(jì)算中, 我們選擇的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)注每步力的計(jì)算的負(fù)載均衡和計(jì)算的絕對(duì)性能. 在下面的測試中Load imb%均在2.5%以下,因此不做單獨(dú)記錄分析.

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

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

圖10 3MGO和DPPC體系CPU和GPU計(jì)算性能對(duì)比圖

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

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

4 總結(jié)

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

我們從原理闡述到實(shí)例討論, 為GROMACS的終端用戶提供了性能優(yōu)化的方向和依據(jù), 幫助其更容易理解這個(gè)軟件黑盒子并行優(yōu)化機(jī)理, 乃至其它分子動(dòng)力學(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)會(huì)(2016156);“一三五”規(guī)劃重點(diǎn)培育方向?qū)ｍ?xiàng)項(xiàng)目(CNIC_PY_1404)

2016-03-24;收到修改稿時(shí)間:2016-04-27

10.15888/j.cnki.csa.005471