基于GPU的碳納米管分子動力學并行仿真

孟小華，覃大勝，鄭冬琴，周玉宇

(暨南大學 a．計算機科學系;b．物理系，廣州510632)

1 概述

自從文獻［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)分成多層可并行處理的單元進行并行計算。

2 分子動力學Verlet仿真算法

分子動力學是一套分子模擬方法，是一種重要且有效的原子尺度計算機模擬手段。該方法主要依靠牛頓力學來模擬分子體系的運動，以在由分子體系的不同狀態(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)隨著能量傳送的溫度變化曲線圖。

3 基于CUDA并行仿真算法的研究與實現(xiàn)

3．1 碳納米管系統(tǒng)模型的建立

在數(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。即有:，由此可得

3．2 并行仿真算法研究

在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 碳納米管分子動力學并行仿真算法的流程

4 GPU并行算法的優(yōu)化

Verlet仿真程序是一個計算力→計算位置→計算速度→計算力……的不斷循環(huán)過程，每個粒子獨立運行一個線程，然后循環(huán)成千上萬次甚至更多，因此小部分的程序優(yōu)化也會給整個程序帶來極大的性能提升。對于CUDA的程序而言，優(yōu)化主要集中在指令優(yōu)化、存儲器訪問優(yōu)化、網(wǎng)格優(yōu)化以及資源均衡等。

4．1 指令優(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的冪次方時以移位運算代替除法。

4．2 存儲器訪問優(yōu)化

合理使用share memory資源，減少不必要的臨時變量，使用精簡函數(shù)以減少register使用。同時，將程序中的常量以及計算中可合并常量保存于constant memory常量存儲器中。Global memory是顯存中容量最大的存儲器，可被任意GPU線程讀寫存儲器任意位置，但由于全局存儲器沒有緩存，具有較高的訪存延遲。因此，程序中在減少不必要全局存儲讀寫的同時，增加處理器占用率，即當一個block訪存時，另一個block可切入運算隱藏延遲。此外，改變顯存的存儲結(jié)構(gòu)，盡量滿足合并訪問存儲器。

4．3 CPU+GPU 處理

考慮到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的并行處理。

5 實驗結(jié)果與分析

本文測試在 Widows 7平臺下進行，GPU1和GPU2是2個不同型號的顯卡，每次實驗只選取其中一塊顯卡。為準備記錄運行時間并方便CPU與GPU并行運行，測試時間以CPU時間為記錄。部分硬件參數(shù)如表1所示。

表1 測試平臺部分硬件參數(shù)

分別使用CPU串行、GeForce GT 430顯卡和Quadro 2000顯卡按不同粒子數(shù)計算，測試三者耗時，如表2所示。

表2 串行并行程序運行于不同硬件平臺時的數(shù)據(jù)

5．1 不同GPU間并行算法性能對比

為使數(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帶來更多的性能提升。

5．2 串并行算法性能對比

對比CPU串行程序與GPU并行程序的性能，計算不同粒子數(shù)下耗時，GPU使用Quadro 2000。

表格數(shù)據(jù)繪制成串并行程序執(zhí)行時間對比如圖8所示?？梢悦黠@看出，當要計算的粒子數(shù)不斷增加時，并行程序相比于串行程序的加速比有一定程度的增加。粒子數(shù)量在1 000～4 000時，由于線程數(shù)據(jù)不足，不能很好地隱藏GPU程序加載時間、數(shù)據(jù)傳輸時間以及訪存時間。當粒子數(shù)量不斷增加時，加速效果越來越好。

圖8 串并行程序執(zhí)行時間對比

5．3 并行算法性能分析

為分析并行算法的加速比，將表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ù)的比重

6 結(jié)束語

在對碳納米管分子動力學進行仿真時，原來的分子動力學仿真算法有良好的效果，但因碳納米管系統(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．