孟小華,覃大勝,鄭冬琴,周玉宇
(暨南大學 a.計算機科學系;b.物理系,廣州510632)
自從文獻[1]首次應用分子動力學研究硬球系統(tǒng)以來,分子動力學廣泛應用于晶體生長、壓痕試驗、摩擦學以及低壓金剛石合成等領域。在新興的納米工程領域,建立在連續(xù)介質(zhì)基礎上的宏觀機理很難解釋納米工程中出現(xiàn)的一些特殊現(xiàn)象,因而分子動力學方法成為重要的研究手段之一[2]。文獻[3]通過分子動力學模擬了Ni針尖在Au(001)表面的納米壓痕過程,其模擬結(jié)果發(fā)表在《科學》上,成為該領域的標志性成果。文獻[4]通過分子動力學仿真算法,研究了不同晶粒尺寸的多晶銅塑性變形過程,從而發(fā)現(xiàn)晶粒尺寸與強度不完全滿足Hall-Petch關系。
分子動力學仿真算法,有著理論分析方法和實驗方法無法比擬的優(yōu)點,但是其計算要求極高,往往受限于計算機的計算能力。為了提高仿真算法規(guī)模,國內(nèi)外許多專家學者做了大量研究[5-7]。仿真算法已由最初提高單機規(guī)模的串行算法[8],發(fā)展到如今通過將計算任務分配給多個CPUs,從而提高模擬規(guī)模的并行算法[9],規(guī)模也由最初的幾千個原子提高到上百萬、上千萬甚至到上億個原子。由于碳納米管[10]系統(tǒng)中含有大量粒子,結(jié)構(gòu)和受力復雜[11],其分子動力學模擬需要極強的計算能力,對CPU處理能力也有較高的要求。而利用CPU實現(xiàn)仿真串行算法,計算量極大,運算時間長,并且模擬粒子數(shù)也受到很大的限制。
本文基于CUDA平臺,通過調(diào)度GPU多個流處理單元[12-13],利用GPU強大的并行計算能力和高效的數(shù)據(jù)傳輸能力,對碳納米管分子動力學的并行仿真算法進行研究和實現(xiàn),通過把碳納米管系統(tǒng)分成多層可并行處理的單元進行并行計算。
分子動力學是一套分子模擬方法,是一種重要且有效的原子尺度計算機模擬手段。該方法主要依靠牛頓力學來模擬分子體系的運動,以在由分子體系的不同狀態(tài)構(gòu)成的系統(tǒng)中抽取樣本,從而計算體系的構(gòu)型積分,并以構(gòu)型積分的結(jié)果為基礎進一步計算體系的熱力學量和其他宏觀性質(zhì)。因而,用分子動力學是碳納米管仿真的重要途徑。
在分子動力學Verlet仿真算法中,粒子位置的泰勒展開式如下:
(1)將x(t+Δt)和x(t-Δt)進行泰勒展開如下式所示:
其中,x(t-Δt)表示為前一時刻的位置;x(t+Δt)表示為后一時刻的位置。
(2)將式(1)和式(2)相加,得到位置表達式如下:
(3)將式(1)和式(2)相減,再兩邊同時除以2Δt,獲得速度的表達式如下:
(4)由式(1)~式(4),在已知粒子 t-2Δt時刻的位置 x(t-2Δt)、t-Δt時刻的位置 x(t-Δt)和t-Δt時刻的加速度a(t-Δt)的情況下,啟動Verlet算法進行積分:根據(jù) x(t-2Δt),x(t-Δt)和 a(t-Δt),將 t=t-Δt代入式(3),獲得 x(t);根據(jù) x(t),基于一定的勢函數(shù)更新a(t);同時,根據(jù)x(t)和x(t-2Δt),將t=t-Δt代入式(4),更新v(t-Δt);即得到粒子x(t),v(t-Δt)和a(t),重復執(zhí)行該步驟。
由上,可利用Verlet列表法設計串行仿真算法的實現(xiàn)步驟如下:
(1)先構(gòu)造如圖1所示的碳納米管結(jié)構(gòu)模型,并設置各種參值。
(2)設計前端熱浴和后端熱浴。
(3)構(gòu)造并調(diào)用計算碳納米管結(jié)構(gòu)內(nèi)部的力的結(jié)構(gòu)函數(shù)及計算碳納米球C60與碳納米管間范德華力的函數(shù)等。
(4)計算所有的力,解牛頓運動方程,分層計算并積分(算法的核心)。
(5)模擬出碳納米球C60隨著熱浴時間的運行軌跡,根據(jù)運算結(jié)果,可畫出和碳納米管內(nèi)隨著能量傳送的溫度變化曲線圖。
在數(shù)據(jù)庫文件中讀取碳納米管系統(tǒng)的參數(shù)條件,如初始溫度、粒子數(shù)、密度時間等,構(gòu)造碳納米管系統(tǒng)模型MOD1及參數(shù)模型MOD2,模型MOD2可以設置各種參數(shù)在模型MOD1上使用,模型MOD1如圖1~圖3所示,由碳納米管和設置在碳納米管內(nèi)部的一個足球狀C60分子(由60個碳原子構(gòu)成的分子,形似足球,又名足球烯)組成,所述碳納米管為層狀中空結(jié)構(gòu),管身是準圓管結(jié)構(gòu),由多個六邊形碳環(huán)結(jié)構(gòu)單元組成,其直徑一般在一到幾十個納米之間,長度則遠遠大于其直徑,對該模型進行初始化,讀取模型中所有粒子的速度和位置坐標。
圖1 碳納米管
圖2 C60分子
圖3 碳納米管系統(tǒng)模型
在碳納米管中,除管口層粒子外,其余粒子相互組成正六邊形。如圖4所示,非管口粒子都有3個近鄰粒子(ip,jp,kp均為idx粒子的近鄰粒子)。
圖4 碳納米結(jié)構(gòu)
在碳納米管系統(tǒng)中,因為C60分子受到碳納米管管壁上碳粒子產(chǎn)生的力,在平衡位置上來回振動,當C60分子處于平衡狀態(tài)時,采用Gaussian熱浴法(約束溫度調(diào)節(jié)方法),在碳納米管兩端同時分別設置不同溫度的熱浴條件,并設定熱浴時間,使碳納米管管壁上碳粒子和C60分子因平衡狀態(tài)被打破而發(fā)生位置和速度的變化。而這些變化是相互影響的,具有傳遞性,伴隨著這些粒子位置和速度的變化,能量在長管狀的碳納米管系統(tǒng)內(nèi)發(fā)生傳遞。由于C60分子位于碳納米管內(nèi),其受力會隨著碳納米管組成粒子的位置改變而改變,因此伴隨著能量傳遞的過程,C60分子的運動狀態(tài)和能量也會發(fā)生改變;采用Gaussian熱浴法的基本原理是在運動方程中加入摩擦力fi,并將其與粒子速度vi聯(lián)系起來,其受力公式為:
當平衡態(tài)時,系統(tǒng)溫度不變,因此有dEk/dt=0。即有:,由此可得
在CUDA平臺上對碳納米管進行分割,分割成多層大小適宜、相互獨立的計算單元,分割原則是在避免過多重復計算的前提上提高并行度,即分割的計算單元必須合理粗細化。因為計算單元過大,則并行不明顯,計算單元過小,會導致太多不必要的重復計算。由于粒子間的作用力間距影響較大,必須設定一個距離臨界值,以判斷其位置關系,當粒子間的間距小于距離臨界值,為近鄰粒子;當粒子間的間距大于距離臨界值,為次鄰粒子,則認為其相互間的分子作用力可以忽略不計。
可以預料到的是每個中心粒子最多只會對其3個近鄰及3個近鄰的各2個近鄰,即6個次鄰施力??梢赃@樣來考慮,以中心粒子以及3個近鄰6個次鄰劃分到一組計算。分批次并行計算。在某一批次計算中,未當成中心粒子計算的一個粒子加入可并行計算隊列中,由于競爭關系,其3個近鄰6個次鄰在本批次中將不能加入隊列,同樣的,對這3個近鄰6個次鄰施力的其他所有未計算粒子也不能加入隊列中。這樣,同一批次的粒子將可以并行計算而不產(chǎn)生競爭關系。
設計分裂算法,采用CPU遍歷計算單元得到n個可并行運算隊列:
(1)若所有粒子均已作為中心計算粒子,跳到步驟(8)。
(2)按次序找到第一個非競爭粒子并沒有按中心計算的粒子,加入可并行運算隊列。
(3)標記該粒子的所有近鄰為本并行隊列一度不可并行粒子,若近鄰已經(jīng)是二度不可并行粒子,則提升為一度不可計算。
(4)標記該粒子所有次鄰為本并行隊列二度不可并行粒子,若次鄰為一度不可并行粒子,則不修改其度數(shù)。
(5)標記本粒子為已按中心計算粒子。
(6)判斷是否已遍歷到粒子隊尾。若是繼續(xù)執(zhí)行,若否跳到步驟(1)。
(7)下一個可并行隊列開始,返回步驟(2)。
(8)結(jié)束。
分裂算法在CPU預處理階段完成,只計算一次得出的可并行隊列可以在后面的循環(huán)計算中一直使用。當循環(huán)次數(shù)較多時,分裂算法所使用的時間占總時間比例將大大減少,對程序總體性能影響也降到比較低的程度。但是由于算法核心是批次計算,因此計算力時將多次啟動內(nèi)核計算函數(shù),一定程度上加大了運行開銷。
如圖5和圖6所示,調(diào)度GPU的流處理單元,采用Verlet算法進行并行運算和處理:
(1)按碳納米管系統(tǒng)模型所有粒子的位置,計算近鄰粒子以及次鄰粒子間的鍵關系和角度關系。
(2)調(diào)度GPU的流處理單元,并行計算被分割在不同計算單元里的粒子,并積累計算碳納米管管壁上每個粒子與其近鄰粒子的相互作用力。
(3)根據(jù)C60分子所在區(qū)域,計算C60分子內(nèi)每個粒子與其所在碳納米管區(qū)域中的每個粒子的相互作用力(范德華力)。
(4)根據(jù)粒子所受到的力及其速度,更新粒子位置,再執(zhí)行步驟(2)和步驟(3)。
(5)根據(jù)粒子所受到的力,計算粒子本次速度及碳納米管前、后兩端熱浴的熱流值。
(6)若達到循環(huán)頻數(shù),則計算結(jié)束;否則,在CPU端間隔保存粒子的數(shù)據(jù),返回步驟(4)。
圖5 碳納米管分子動力學并行仿真中的調(diào)度模式
圖6 碳納米管分子動力學并行仿真算法的流程
Verlet仿真程序是一個計算力→計算位置→計算速度→計算力……的不斷循環(huán)過程,每個粒子獨立運行一個線程,然后循環(huán)成千上萬次甚至更多,因此小部分的程序優(yōu)化也會給整個程序帶來極大的性能提升。對于CUDA的程序而言,優(yōu)化主要集中在指令優(yōu)化、存儲器訪問優(yōu)化、網(wǎng)格優(yōu)化以及資源均衡等。
目前,GPU單精度計算性能要遠遠超過雙精度計算性能,因此,程序中對精度要求不高的部分可使用單精度運算代替雙精度運算。在GPU中整數(shù)乘法、除法、求模運算的指令吞吐量較為有限,而控制流指令由于會影響指令發(fā)射器的并行發(fā)射,打斷指令流開銷巨大。此時,在可知整數(shù)范圍中以經(jīng)優(yōu)化的庫24位整數(shù)運算__mul24,__umul24代替32位整數(shù)運算,以__sinf(x),__fdivide(x,y)等代替相應的運算。盡量避免整數(shù)除和模運算,或以(i&(n-1))代替(i%n),除數(shù)是2的冪次方時以移位運算代替除法。
合理使用share memory資源,減少不必要的臨時變量,使用精簡函數(shù)以減少register使用。同時,將程序中的常量以及計算中可合并常量保存于constant memory常量存儲器中。Global memory是顯存中容量最大的存儲器,可被任意GPU線程讀寫存儲器任意位置,但由于全局存儲器沒有緩存,具有較高的訪存延遲。因此,程序中在減少不必要全局存儲讀寫的同時,增加處理器占用率,即當一個block訪存時,另一個block可切入運算隱藏延遲。此外,改變顯存的存儲結(jié)構(gòu),盡量滿足合并訪問存儲器。
考慮到GPU并行處理要發(fā)揮其強大的并行計算能力,要有足夠多的線程并行運算。因此,把部分較小計算量的計算交給CPU執(zhí)行。如計算C60與炭納米管間范德華力函數(shù),由于C60只有60個粒子,相對來說并行度小,將其計算放到CPU里。由于內(nèi)核的啟動是異步的,因此把CPU函數(shù)放到內(nèi)核啟動函數(shù)下面,CPU啟動內(nèi)核執(zhí)行GPU運算時可立即返回執(zhí)行CPU程序,實現(xiàn)CPU與GPU的并行處理。
本文測試在 Widows 7平臺下進行,GPU1和GPU2是2個不同型號的顯卡,每次實驗只選取其中一塊顯卡。為準備記錄運行時間并方便CPU與GPU并行運行,測試時間以CPU時間為記錄。部分硬件參數(shù)如表1所示。
表1 測試平臺部分硬件參數(shù)
分別使用CPU串行、GeForce GT 430顯卡和Quadro 2000顯卡按不同粒子數(shù)計算,測試三者耗時,如表2所示。
表2 串行并行程序運行于不同硬件平臺時的數(shù)據(jù)
為使數(shù)據(jù)形象化,將GPU1與GPU2數(shù)據(jù)繪制成條狀圖。由圖7可得,前期粒子數(shù)較少時,2種顯卡耗時相差不大,但隨著粒子數(shù)的增加,GeForce GT 430顯卡耗時增長速度逐漸快于Quadro 2 000顯卡。而當粒子數(shù)量增加到120 000個時,前者耗時幾乎是后者的2倍。
圖7 2種并行顯卡運行程序耗時對比
由表1可知,Quadro 2 000有著2倍于GeForce GT 430的 CUDA Cores(即流處理器 stream processor)。結(jié)合GPU的線程運行方式,當GPU的處理單元增多,同時處于運行狀態(tài)的線程也會相應增加。因此,當無其他硬件限制時,GPU并行計算的性能幾乎隨著GPU處理器的增加而線性增加,使用擁有更多處理器的GPU能給Verlet帶來更多的性能提升。
對比CPU串行程序與GPU并行程序的性能,計算不同粒子數(shù)下耗時,GPU使用Quadro 2000。
表格數(shù)據(jù)繪制成串并行程序執(zhí)行時間對比如圖8所示??梢悦黠@看出,當要計算的粒子數(shù)不斷增加時,并行程序相比于串行程序的加速比有一定程度的增加。粒子數(shù)量在1 000~4 000時,由于線程數(shù)據(jù)不足,不能很好地隱藏GPU程序加載時間、數(shù)據(jù)傳輸時間以及訪存時間。當粒子數(shù)量不斷增加時,加速效果越來越好。
圖8 串并行程序執(zhí)行時間對比
為分析并行算法的加速比,將表2中串并行程序運算時間按公式:
計算,GPU使用Quadro 2000,并將結(jié)果繪制成GPU并行計算加速比,如圖9所示。
圖9 GPU并行計算加速比
可以看出,隨著粒子數(shù)的增加,GPU并行算法相比于串行算法的加速比也隨之增大。當粒子數(shù)到達120 000個時,加速比達到1 400%,即14倍的加速效果。
為分析影響GPU并行程序性能因素,使用性能評估工具NVIDIA Visual Profiler對程序進行分析。目標程序計算600層12 000個粒子,每個block中有128個線程。首先分析各函數(shù)運行時間占總運算時間的比例。通過多次運算求平均值,如圖10所示。
圖10 各個函數(shù)的比重
在對碳納米管分子動力學進行仿真時,原來的分子動力學仿真算法有良好的效果,但因碳納米管系統(tǒng)粒子數(shù)的規(guī)模太大、計算量多、運行時間長,缺乏實用性。而由用Verlet算法實現(xiàn)的分子動力學仿真算法具有易并行性,可以在具有強大圖像處理和浮點計算能力的GPU上并行處理,從而大大提高了碳納米管分子動力學仿真算法的效率。
實驗結(jié)果表明,對小規(guī)模僅有120 000個粒子組成的碳納米管系統(tǒng),在擁有192個流處理器的NVIDIA Quadro 2000顯卡上實現(xiàn)的GPU并行仿真算法,與在擁有4核8線程的Intel Xeon W3550上實現(xiàn)的基于CPU的串行仿真算法相比,其加速比即可達到十多倍,從而有效地解決了算法速度上的缺陷。
然而由于本次實驗的硬件條件有限,所能計算的碳納米管系統(tǒng)規(guī)模有限,因此無法深度發(fā)掘CUDA的能力。而由實驗結(jié)果所示,碳納米管系統(tǒng)的粒子數(shù)規(guī)模越大,GPU并行算法對于串行算法的加速比就越明顯。在實際運用中,碳納米管系統(tǒng)的粒子數(shù)的規(guī)模都很大,超過千萬乃至十億,可以預估該基于GPU的分子動力學并行仿真算法,對碳納米管的研究會發(fā)揮至關重要的作用。
[1] Alder B J,Wainwright T E.Phase Transition for a Hard Sphere System[J].The Journal of Chemical Physics,1957,27(5):1208-1209.
[2] 唐玉蘭,胡 適,王東旭,等.納米工程中大規(guī)模分子動力學仿真算法的研究進展[J].機械工程學報,2008,44(2):8-15.
[3] Landman U,Luedtke W D,Burnham N,et al.Atomistic Mechanisms and Dynamics of Adhesion,Nanoindentation,and Fracture[J].Science,1990,248(4954):454-461.
[4] Schi?tz J,Jacobsen K W.A Maximum in the Strength of Nanocrystalline Copper[J].Science,2003,301(5638):1357-1359.
[5] 徐 偉,王 麗.MD算法在集群系統(tǒng)中的并行實現(xiàn)研究[J].計算機工程,2002,28(3):103-105.
[6] Proctor A J,Lipscomb T J,Zou Anqi,et al.Performance Analyses of a Parallel Verlet Neighbor List Algorithm for GPU-optimized MD Simulations[C]//Proceedings of International Conference on Biomedical Computing.Washington D.C.,USA:IEEE Press,2012:14-19.
[7] Anderson J A,Lorenz C D,TravessetA.General Purpose Molecular Dynamics Simulations Fully Implemented on Graphics Processing Units[J].Journal of Computational Physics,2008,227(10):5342-5359.
[8] Verlet L.Computer Experiments on Classical Fluids.I.Thermodynamical Properties of Lennard-Jones Molecules[J].Physical review,1967,159(1):98-103.
[9] Tamayo P,MesirovJP,BoghosianB M.Parallel Approaches to ShortRange Molecular Dynamics Simulations[C]//Proceedings of Conference on Supercomputing.New York,USA:ACM Press,1991:462-470.
[10] 曹 偉,宋雪梅,王 波,等.碳納米管的研究進展[J].材料導報,2007,21(5):77-82.
[11] 劉文志,李曉霞,余 翔,等.分子動力學模擬中基于GPU的范德華非鍵作用計算[J].計算機與應用化學,2010,27(12):1607-1612.
[12] 孟小華,劉堅強,區(qū)業(yè)祥,等.基于CUDA的拉普拉斯邊緣檢測算法[J].計算機工程,2012,38(18):190-193.
[13] 孟小華,黃叢珊,朱麗莎.基于CUDA的熱傳導GPU并行算法研究[J].計算機工程,2014,40(5):41-44,48.