李 靖,祝愛琦,韓 林,侯超峰
(1.鄭州大學(xué) 信息工程學(xué)院,鄭州 450001;2.中國科學(xué)院過程工程研究所,北京 100190;3.鄭州大學(xué) 國家超級計算鄭州中心,鄭州 450001)
近年來,利用圖形處理器(Graphics Processing Unit,GPU)加速分子動力學(xué)(Molecular Dynamics,MD)模擬已成為計算化學(xué)、計算材料學(xué)等領(lǐng)域的研究熱點[1-2]。分子動力學(xué)模擬通常從原子、分子尺度上研究凝聚態(tài)材料的熱力學(xué)性質(zhì)和行為,被廣泛應(yīng)用于微觀計算模擬領(lǐng)域中[3]。晶體硅是信息技術(shù)產(chǎn)業(yè)中應(yīng)用最廣泛和最重要的材料之一,其納微結(jié)構(gòu)的熱傳遞性質(zhì)和導(dǎo)熱性能會顯著影響電子設(shè)備的穩(wěn)定性,導(dǎo)致器件性能降低。因此,分子動力學(xué)模擬方法對接近真實尺度硅納微結(jié)構(gòu)的直接模擬的實現(xiàn)與其導(dǎo)熱機制的研究具有重要意義。
模擬體系的計算量和時空分辨率的提高是分子動力學(xué)模擬領(lǐng)域的研究熱點。高性能加速芯片(如GPU)為處理模擬體系的巨量計算提供契機。GPU作為一種具有良好并行性能且適應(yīng)密集計算的新型加速計算工具,在通用科學(xué)計算領(lǐng)域具有較大的應(yīng)用潛力。分子動力學(xué)模擬優(yōu)化的研究大多集中在通用分子動力學(xué)模擬算法優(yōu)化和程序并行算法設(shè)計方面,很少專門針對固體材料并結(jié)合GPU 硬件結(jié)構(gòu)特性對MD 算法進行設(shè)計和優(yōu)化。文獻[4]提出一種高效的多粒子碰撞動力學(xué)算法,并在多個GPU 上實現(xiàn)。針對多體勢的高效分子動力學(xué)模擬,文獻[5]提出一種基于GPU 的力求解算法。文獻[6]實現(xiàn)分子動力學(xué)三體勢的GPU 加速計算。文獻[7]提出不同于對勢分子動力學(xué)模擬的方法,根據(jù)復(fù)雜的多體勢描述一種GPU 加速的固態(tài)共價晶體分子動力學(xué)模擬的高效和可擴展算法。
針對固態(tài)共價晶體硅的分子動力學(xué)模擬,本文結(jié)合原子間多體勢函數(shù)與CUDA 編程模型,通過對GPU 并行計算特征和模擬熱點進行分析,提出面向GPU 計算平臺的固定鄰居算法設(shè)計與優(yōu)化。根據(jù)GPU 硬件架構(gòu)特性,將并行線程處理的float4 數(shù)據(jù)結(jié)構(gòu)重新映射為連續(xù)float 布置的形式,采用三目運算符控制指令排隊,消除條件分支語句造成的流水線停頓,減少全局訪存消耗和分支結(jié)構(gòu)耗時,提升算法的計算效率。
為充分發(fā)揮GPU 的計算性能,從軟件與硬件相結(jié)合的角度分析固態(tài)材料分子動力學(xué)模擬算法。以Tesla V100 GPU 為實驗平臺描述GPU 硬件架構(gòu)[8-10]。Tesla V100 GPU 架構(gòu)如圖1 所示,由多個圖形處理簇(Graphics Processing Cluster,GPC)、紋理處理簇(Texture Processing Cluster,TPC)、流式多處理器(Streaming Multiprocessor,SM)以及內(nèi)存控制器組成。每個SM 內(nèi)有64 個流處理器(Streaming Processor,SP),其中,每個SP 具有自己的緩存空間和寄存器組,并通過線程束(warp)調(diào)度器的方式控制SP 的并行執(zhí)行,將L1 高速緩存和共享內(nèi)存組合到單個內(nèi)存塊中,內(nèi)存總?cè)萘繛?28 KB,在SM 上運行的線程對L1 高速緩存和共享內(nèi)存的訪問相對于全局內(nèi)存具有更快的訪問速度,在不同SM 之間的L1 高速數(shù)據(jù)緩存和共享內(nèi)存相互獨立。Tesla V100 GPU 總共擁有80 個SM、5 120 個32 位浮點處理單元、2 560 個64 位浮點處理單元、640 個Tensor 核心以及320 個紋理單元(Tex)[8]。
圖1 Tesla V100 GPU 架構(gòu)Fig.1 Architecture of Tesla V100 GPU
在硬件與執(zhí)行程序的對應(yīng)關(guān)系上,一個SP 控制一個線程的執(zhí)行。為提高并行效率,NVIDIA 把32 個線程組織成線程束的形式,并在warp 內(nèi)以單指令多線程(Single Instruction Multiple Threads,SIMT)的方式執(zhí)行。在線程塊中的多個warp 由線程束調(diào)度器調(diào)度進入SM 執(zhí)行[11]。在Tesla V100 GPU 中大量的SM 雖然提高GPU 硬件的并行性,但是其具有較高并行性的復(fù)雜架構(gòu)給MD 模擬算法的優(yōu)化帶來挑戰(zhàn)。在線程并行優(yōu)化方面,不僅減少線程之間的相互依賴關(guān)系,還需要充分利用其硬件特性以減少訪存開銷。
MD 模擬是通過刻畫微觀原子的運動過程,由相關(guān)的數(shù)據(jù)統(tǒng)計來研究模擬體系的結(jié)構(gòu)和熱力學(xué)性質(zhì)的方法。其中,在周圍鄰居原子勢場作用下每個原子根據(jù)牛頓定律運動[12]。原子間的相互作用是MD 模擬算法的主要組成部分。與原子間對勢相比,原子間多體勢具有更復(fù)雜的模型形式,但通常能更準(zhǔn)確、更有效地體現(xiàn)共價材料模擬體系的物理特性[13]。
本文針對固態(tài)共價晶體硅,基于固定鄰居算法,利用Tersoff 多體勢實現(xiàn)晶體硅的MD 模擬計算。Tersoff 勢函數(shù)的數(shù)學(xué)表達式如式(1)和式(2)所示:
其中:UR為排斥勢項;UA為吸引勢項;fC為截斷函數(shù);bij為鍵級項。UR、UA、fC、bij、ζij和g(θ)的數(shù)學(xué)表達式如式(3)~式(8)所示:
其中:ζij為原子i與原子j間的角勢能。上述公式的模型參數(shù)見文獻[14-15]。固態(tài)共價晶體硅的MD 模擬計算首先需要生成模擬體系晶體硅原子的初始位置,并在特定溫度下給出原子的初始速度,通過Tersoff 多體勢計算每個硅原子受其他硅原子的作用合力,然后選擇蛙跳算法積分牛頓運動方程,求解硅原子的新速度和新位置。
固態(tài)共價晶體硅具有規(guī)則的晶胞結(jié)構(gòu),在低溫狀態(tài)下每個原子周圍有4 個鄰居原子。在Tesla V100 GPU 平臺實現(xiàn)算法時,模擬體系中所有的硅原子按超胞[7]的順序在分配的設(shè)備存儲器上排列,并且在后續(xù)的迭代模擬中按超胞順序分配的內(nèi)存位置保持不變。超胞中原子的分組情況如圖2 所示。根據(jù)每個原子的鄰居原子在超胞內(nèi)的數(shù)目,把超胞內(nèi)的原子分為三組:第一組原子在超胞內(nèi)有一個鄰居原子,如圖2(a)中的灰色原子;第二組原子在超胞內(nèi)有兩個鄰居原子;第三組原子在超胞內(nèi)有四個鄰居原子,如圖2(a)中的黑色原子。算法執(zhí)行時將超胞內(nèi)上述三組原子按順序布置,如圖2(b)所示。此外,為提高緩存效率,采用靜態(tài)重排方法代替動態(tài)更新鄰居列表的過程,根據(jù)超胞內(nèi)原子分組情況,把原子在超胞內(nèi)和超胞外的鄰居原子信息分組聚合并進行連續(xù)排列。由于固定鄰居算法沒有鄰居列表動態(tài)更新的過程,因此原子間作用力的計算、原子位置和速度的更新占大量的模擬總時間。
圖2 每個超胞中計算原子的分組示意圖Fig.2 Schematic diagram of grouping of calculated atoms in each supercell
固定鄰居算法流程如圖3 所示,在CPU 端通過控制模擬的原子數(shù)目和迭代步數(shù),初始化原子的位置、速度等信息,將原子按其下標(biāo)映射到與超胞排序相應(yīng)的CPU 端內(nèi)存數(shù)組中,然后完成CPU 到GPU 的數(shù)據(jù)拷貝以及GPU 端核函數(shù)的生成和調(diào)度。在GPU 端,首先將原子按其下標(biāo)映射到超胞排序的相應(yīng)GPU 數(shù)組中,根據(jù)“一個線程模擬一個原子”的分配原則,使超胞中的原子數(shù)等于線程塊中的線程數(shù),并利用線程分支控制超胞內(nèi)的分組情況。在GPU設(shè)備端并行計算每個超胞中每個原子和其鄰居原子之間的距離,進一步計算每個原子受鄰居原子的全部作用力。最后,通過積分時間計算原子的新速度和新位置。當(dāng)算法沒有達到模擬設(shè)定的迭代步數(shù)時,GPU 端繼續(xù)執(zhí)行MD 模擬,直到達到設(shè)定的模擬迭代步數(shù),程序結(jié)束。
圖3 固定鄰居算法流程Fig.3 Procedure of fixed neighbor algorithm
雖然固定鄰居算法采用超胞分組的方式提高緩存命中率,但是在設(shè)備端并行執(zhí)行原子受力計算時也增加了if 分支,為GPU 程序并行執(zhí)行帶來了額外的分支消耗。Tersoff 多體勢比對勢具有更復(fù)雜的模型形式,當(dāng)計算原子間作用力時使用數(shù)據(jù)類型(float4)并未實現(xiàn)計算吞吐量和數(shù)據(jù)全局訪存的最佳匹配。因此,為提高程序的性能,本文結(jié)合GPU硬件架構(gòu)特性,對固定鄰居算法進行設(shè)計優(yōu)化。
在對Tesla V100 GPU 平臺的硬件特性進行分析之后,本文分別對固定鄰居算法的數(shù)據(jù)結(jié)構(gòu)和分支結(jié)構(gòu)進行優(yōu)化,充分利用GPU 計算資源來提高線程并行效率,進一步提升固定鄰居算法的性能。
早期的GPU 適用于圖形渲染設(shè)計,其運算對象是表示顏色的R、G、B、A 和表示坐標(biāo)的X、Y、Z、W,都是對一個四維向量進行操作[16]。隨著GPU 應(yīng)用越來越廣泛,數(shù)據(jù)結(jié)構(gòu)可能不僅具有表示坐標(biāo)的X、Y、Z、W和表示顏色的R、G、B、A,還有可能存在一些非向量的數(shù)據(jù)。因此,本文對不同的數(shù)據(jù)做不同的處理,而不同的數(shù)據(jù)結(jié)構(gòu)方式會影響數(shù)據(jù)的內(nèi)存局部性。為提高通用性,NVIDIA 引入標(biāo)量計算的概念,當(dāng)發(fā)展到Tesla V100 時,其內(nèi)部包含大量專門針對標(biāo)量的計算單元(CUDA CORE),大幅提高運算性能。
在固定鄰居算法中,數(shù)據(jù)結(jié)構(gòu)主要包括各個原子的位置、速度、力等屬性,在三維空間中,它們都是以四維向量方式定義的,四個分量x、y、z、w分別表示在三個方向上的方位坐標(biāo)和空值。線程對數(shù)據(jù)結(jié)構(gòu)的重映射示意圖如圖4 所示。float4 數(shù)據(jù)結(jié)構(gòu)如圖4(a)所示,原子屬性以四維向量x、y、z、w的形式連續(xù)存儲在全局內(nèi)存中,從三個維度分別計算原子間的距離和受力。當(dāng)計算x方向上兩個原子間的距離和受力時,在warp 中,如Thread 0 訪問x1,Thread 1訪問x2,x1與x2之間相隔y1、z1、w1,兩個線程對原子數(shù)據(jù)的訪問和計算會造成計算吞吐量的降低。在圖4(b)中,結(jié)合Tesla V100 GPU 的硬件特性,修改算法數(shù)據(jù)結(jié)構(gòu)布局,使原子屬性以單位標(biāo)量x、y、z、w的形式各自分開存儲在全局內(nèi)存中,在計算原子間x方向上的距離和受力時,在warp 中,Thread 0 訪問x1,Thread 1 訪問x2,類似操作實現(xiàn)了線程對原子數(shù)據(jù)的連續(xù)訪問。這種新的數(shù)據(jù)結(jié)構(gòu)布局方式有助于充分利用硬件特性,提升計算吞吐量和訪存帶寬性能,提高算法的運行效率。
圖4 線程對數(shù)據(jù)結(jié)構(gòu)的重映射示意圖Fig.4 Remapping schematic diagram of data structures by threads
GPU 分支結(jié)構(gòu)執(zhí)行方式不同于CPU 分支結(jié)構(gòu)執(zhí)行方式。在CPU 中設(shè)置分支預(yù)測機制,有效縮短執(zhí)行時間。在GPU 執(zhí)行中,如果warp 中任何一個線程進入到不同的分支路徑,則其余的線程都處于等待狀態(tài),直到該線程執(zhí)行完分支程序。warp 中線程按順序串行通過多條不同分支路徑,當(dāng)所有分支路徑都執(zhí)行后,warp 中的所有線程才會回到同一條執(zhí)行路徑上。
warp 分支結(jié)構(gòu)如圖5 所示。一個線程束warp 包含32 個線程,在warp 內(nèi)以SIMT 的方式執(zhí)行[17]。在執(zhí)行if 分支時,當(dāng)前18 個線程指向的數(shù)據(jù)執(zhí)行時,剩余后14 個線程指向的數(shù)據(jù)都需要等待,造成線程資源浪費。若程序中的分支路徑較多,且warp 中的線程都需要執(zhí)行時,會延長程序的執(zhí)行時間,造成性能的嚴(yán)重損失。
圖5 warp 分支結(jié)構(gòu)Fig.5 Structure of warp branch
大多數(shù)對于分支結(jié)構(gòu)的優(yōu)化保留在線程層,線程交換是通過減少同一個warp 中不同分支的方式來提高程序并行性,根據(jù)分支控制數(shù)據(jù)的依賴關(guān)系,重新排列數(shù)據(jù)數(shù)組,建立線程與數(shù)據(jù)的重映射關(guān)系,從而提升性能[18-20]。然而,一些特定算法的分支結(jié)構(gòu)無法通過上述方式進行優(yōu)化,本文引入三目運算符的方法對分支結(jié)構(gòu)進行優(yōu)化。分支結(jié)構(gòu)轉(zhuǎn)換為三目運算符的示意圖如圖6 所示。
圖6 分支結(jié)構(gòu)轉(zhuǎn)換為三目運算符的示意圖Fig.6 Schematic diagram of branch structure converted to the ternary operator
固定鄰居算法在超胞中的原子受力需要分組賦值的情況下,采用三目運算符代替分支結(jié)構(gòu)。在執(zhí)行分支結(jié)構(gòu)代碼時,if-else 語句會對分支指令進行排隊處理,而三目運算符對分支的兩種情況同時進行預(yù)處理,使三目運算符代碼的耗時周期相較于分支結(jié)構(gòu)代碼更短,有效減少分支耗時并提升算法的運行效率。
本文固態(tài)共價晶體硅的MD 模擬是在Tesla V100 GPU 單卡上實現(xiàn)并行執(zhí)行。固定鄰居算法的應(yīng)用環(huán)境如表1 所示。采用算法性能最優(yōu)的超胞尺寸,其包含256 個硅原子,在x、y、z方向上長度分別為2、4、4 個單位晶胞,模擬的初始溫度為300 K。
表1 固定鄰居算法的計算環(huán)境Table 1 Computational environment of fixed neighbor algorithm
4.2.1 正確性分析
本文通過算例驗證固定鄰居算法優(yōu)化前后的有效性和正確性,總的迭代步數(shù)設(shè)置為100 萬步,原子模擬數(shù)目為32 768 個,每間隔1 000 迭代步,輸出一次統(tǒng)計結(jié)果。GPU 固定鄰居算法優(yōu)化前后每個原子的平均勢能隨時間變化曲線如圖7 所示,SP 表示單精度,DP 表示雙精度。
圖7 GPU 固定鄰居算法優(yōu)化前后每個原子的平均勢能隨時間變化曲線Fig.7 Time varying curves of average potential energy of each atom before and after optimizing the fixed neighbor algorithm on GPU
在Tesla V100 GPU 上固定鄰居算法優(yōu)化前后原子的平均勢能變化穩(wěn)定,雙精度的原子勢能趨向于-4.610 24 eV,單精度的原子勢能趨向于-4.610 20 eV,模擬結(jié)果與文獻[7]一致,證明了固定鄰居算法優(yōu)化前后的計算正確可靠。
4.2.2 性能測試與分析
在超胞尺寸固定不變的情況下,本文考慮到Tesla V100 GPU 單卡硬件資源的限制,分別對不同原子數(shù)目進行模擬,實驗選擇原子數(shù)目為32 768、262 144、1 024 000、4 096 000、10 240 000,以模擬1 000 迭代步的運行時間作為評價算法性能的指標(biāo),通過對比算法優(yōu)化前后的耗時來驗證優(yōu)化效果。
將數(shù)據(jù)結(jié)構(gòu)與分支結(jié)構(gòu)優(yōu)化方法相結(jié)合,固定鄰居算法優(yōu)化前后的性能對比如表2 所示。t1表示雙精度固定鄰居算法優(yōu)化前的運行時間,t2表示雙精度固定鄰居算法優(yōu)化后的運行時間,t3表示單精度固定鄰居算法優(yōu)化前的運行時間,t4表示單精度固定鄰居算法優(yōu)化后的運行時間。當(dāng)原子數(shù)目為4 096 000時,優(yōu)化后雙精度固定鄰居算法的性能最優(yōu),與固定鄰居算法優(yōu)化前的運行時間相比,雙精度固定鄰居算法與單精度固定鄰居算法的加速效果分別提高了27.7%與20.7%。
表2 固定鄰居算法優(yōu)化前后性能對比Table 2 Performance comparison of fixed neighbor algorithms before and after optimization
固定鄰居算法優(yōu)化后的加速比如圖8 所示。當(dāng)原子數(shù)目為32 768 時,優(yōu)化后雙精度固定鄰居算法加速比是優(yōu)化后單精度固定鄰居算法加速比的1 倍多。其原因為當(dāng)原子規(guī)模較小時,優(yōu)化后雙精度固定鄰居算法計算1 000 迭代步的運行時間較短。此時,單精度固定鄰居算法計算力的核函數(shù)耗時占GPU 端總耗時的比例較小,雙精度固定鄰居算法計算力的核函數(shù)耗時占GPU 端總耗時比例較大,因此,優(yōu)化后雙精度固定鄰居算法的加速比較高。固定鄰居算法從原子數(shù)目為262 144 開始,優(yōu)化后單精度和雙精度固定鄰居算法的加速比隨著原子數(shù)目的增大而增大,并在硬件資源限制的條件下達到一個最高值,然后,加速比開始下降。
圖8 固定鄰居算法優(yōu)化后的加速比Fig.8 Acceleration ratio of fixed neighbor algorithm after optimization
在Tesla V100 GPU 處理器資源可承受范圍內(nèi),固定鄰居算法的優(yōu)化效果依賴于計算的浮點精度和計算規(guī)模。在GPU 計算資源和存儲資源相同的條件下,原子數(shù)目多的體系優(yōu)化效果更好。其原因為當(dāng)原子數(shù)目較小時,浮點計算不能充分利用Tesla V100 GPU 處理器并行資源,計算性能不能達到最優(yōu),當(dāng)原子數(shù)目較大時,Tesla V100 GPU 處理器才能充分發(fā)揮計算優(yōu)勢,加速效果明顯。
LAMMPS 是一種支持多種原子間勢函數(shù)并被廣泛應(yīng)用的開源MD 軟件,利用空間分解技術(shù)將三維模擬區(qū)域劃分為更小的三維子域,利用MPI 并行化提高計算規(guī)模和性能。KOKKOS 加速包是基于LAMMPS 的C++加速庫,可以有效地在GPU 硬件上運行[21]。目前,KOKKOS 加速包沒有精度選項,所有編譯和計算都必須以雙精度執(zhí)行[22-23]。HOOMDblue 是專門支持GPU 加速的分子動力學(xué)模擬軟件,能夠在多種力場下執(zhí)行MD 模擬計算[24-25]。本文使用LAMMPS 軟件KOKKOS 加速包和支持GPU 加速的HOOMD-blue 軟件在Tesla V100 GPU 上對Tersoff多體勢進行單雙精度計算,測試環(huán)境與4.1 節(jié)相同,選擇模擬1 000 迭代步的運行時間作為性能對比的基準(zhǔn)。
LAMMPS 和HOOMD-blue 與優(yōu)化后的固定鄰居算法性能對比如表3 所示,表中T1表示雙精度固定鄰居算法的執(zhí)行時間,T2表示單精度固定鄰居算法的執(zhí)行時間,T3表示LAMMPS 雙精度固態(tài)晶體硅分子動力學(xué)模擬的執(zhí)行時間,T4表示HOOMD-blue 雙精度固態(tài)晶體硅分子動力學(xué)模擬的執(zhí)行時間,T5表示HOOMDblue 單精度固態(tài)晶體硅分子動力學(xué)模擬的執(zhí)行時間,加速比1 表示LAMMPS 雙精度固態(tài)晶體硅分子動力學(xué)模擬與雙精度固定鄰居算法(LAMMPS 雙精度/雙精度算法)的執(zhí)行時間比值,加速比2 表示HOOMD-blue雙精度固態(tài)晶體硅分子動力學(xué)模擬與雙精度固定鄰居算法(HOOMD-blue 雙精度/雙精度算法)的執(zhí)行時間比值,加速比3 表示HOOMD-blue 單精度固態(tài)晶體硅分子動力學(xué)模擬與單精度固定鄰居算法(HOOMD-blue單精度/單精度算法)的執(zhí)行時間比值。
表3 LAMMPS 和HOOMD-blue 與優(yōu)化后的固定鄰居算法性能對比Table 3 Performance comparison of LAMMPS and HOOMD-blue with the optimized fixed neighbor algorithm
相比兩個MD 開源模擬軟件,固定鄰居算法具有較優(yōu)的加速效果。當(dāng)原子數(shù)目為10 240 000 時,雙精度固定鄰居算法加速效果最優(yōu),LAMMPS 和 HOOMD-blue 軟件與雙精度固定鄰居算法的執(zhí)行時間的比值分別為 11.62 和 9.39。當(dāng)原子數(shù)目為 4 096 000 時,單精度固定鄰居算法的加速效果較優(yōu),HOOMD-blue 軟件與單精度固定鄰居算法執(zhí)行時間的比值為12.18。優(yōu)化后的固定鄰居算法相對LAMMPS 和HOOMD-blue 的加速比如圖9 所示。隨著原子數(shù)目的增加,雙精度固定鄰居算法與兩個 MD 開源模擬軟件相比,其加速比呈現(xiàn)上升趨勢。HOOMD-blue 單精度/單精度算法加速比呈現(xiàn)先升高后略微下降的趨勢。
圖9 優(yōu)化后的固定鄰居算法相對LAMMPS 和HOOMD-blue 的加速比Fig.9 Acceleration ratio of the optimized fixed neighbor algorithm to LAMMPS and HOOMD-blue
針對固態(tài)共價晶體硅的分子動力學(xué)模擬,本文結(jié)合GPU 的硬件架構(gòu)特性,提出分子動力學(xué)模擬固定鄰居算法的優(yōu)化設(shè)計,采用數(shù)據(jù)結(jié)構(gòu)和分支結(jié)構(gòu)優(yōu)化方法,減少多線程訪存和分支耗時,提升晶體硅模擬的計算效率。與主要的GPU 加速MD 模擬軟件LAMMPS 和HOOMD-blue 相比,本文固定鄰居算法具有較優(yōu)的加速性能。下一步將對多GPU 并行的分子動力模擬算法進行研究,有效提升固相晶體硅體系的計算規(guī)模,實現(xiàn)數(shù)十億原子以上的模擬。