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

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

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

1 概述

自從文獻［1］首次應用分子動力學研究硬球系統以來，分子動力學廣泛應用于晶體生長、壓痕試驗、摩擦學以及低壓金剛石合成等領域。在新興的納米工程領域，建立在連續介質基礎上的宏觀機理很難解釋納米工程中出現的一些特殊現象，因而分子動力學方法成為重要的研究手段之一［2］。文獻［3］通過分子動力學模擬了Ni針尖在Au(001)表面的納米壓痕過程，其模擬結果發表在《科學》上，成為該領域的標志性成果。文獻［4］通過分子動力學仿真算法，研究了不同晶粒尺寸的多晶銅塑性變形過程，從而發現晶粒尺寸與強度不完全滿足Hall-Petch關系。

分子動力學仿真算法，有著理論分析方法和實驗方法無法比擬的優點，但是其計算要求極高，往往受限于計算機的計算能力。為了提高仿真算法規模，國內外許多專家學者做了大量研究［5-7］。仿真算法已由最初提高單機規模的串行算法［8］，發展到如今通過將計算任務分配給多個CPUs，從而提高模擬規模的并行算法［9］，規模也由最初的幾千個原子提高到上百萬、上千萬甚至到上億個原子。由于碳納米管［10］系統中含有大量粒子，結構和受力復雜［11］，其分子動力學模擬需要極強的計算能力，對CPU處理能力也有較高的要求。而利用CPU實現仿真串行算法，計算量極大，運算時間長，并且模擬粒子數也受到很大的限制。

本文基于CUDA平臺，通過調度GPU多個流處理單元［12-13］，利用GPU強大的并行計算能力和高效的數據傳輸能力，對碳納米管分子動力學的并行仿真算法進行研究和實現，通過把碳納米管系統分成多層可并行處理的單元進行并行計算。

2 分子動力學Verlet仿真算法

分子動力學是一套分子模擬方法，是一種重要且有效的原子尺度計算機模擬手段。該方法主要依靠牛頓力學來模擬分子體系的運動，以在由分子體系的不同狀態構成的系統中抽取樣本，從而計算體系的構型積分，并以構型積分的結果為基礎進一步計算體系的熱力學量和其他宏觀性質。因而，用分子動力學是碳納米管仿真的重要途徑。

在分子動力學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算法進行積分:根據 x(t－2Δt)，x(t－Δt)和 a(t－Δt)，將 t=t－Δt代入式(3)，獲得 x(t);根據 x(t)，基于一定的勢函數更新a(t);同時，根據x(t)和x(t－2Δt)，將t=t－Δt代入式(4)，更新v(t－Δt);即得到粒子x(t)，v(t－Δt)和a(t)，重復執行該步驟。

由上，可利用Verlet列表法設計串行仿真算法的實現步驟如下:

(1)先構造如圖1所示的碳納米管結構模型，并設置各種參值。

(2)設計前端熱浴和后端熱浴。

(3)構造并調用計算碳納米管結構內部的力的結構函數及計算碳納米球C60與碳納米管間范德華力的函數等。

(4)計算所有的力，解牛頓運動方程，分層計算并積分(算法的核心)。

(5)模擬出碳納米球C60隨著熱浴時間的運行軌跡，根據運算結果，可畫出和碳納米管內隨著能量傳送的溫度變化曲線圖。

3 基于CUDA并行仿真算法的研究與實現

3．1 碳納米管系統模型的建立

在數據庫文件中讀取碳納米管系統的參數條件，如初始溫度、粒子數、密度時間等，構造碳納米管系統模型MOD1及參數模型MOD2，模型MOD2可以設置各種參數在模型MOD1上使用，模型MOD1如圖1～圖3所示，由碳納米管和設置在碳納米管內部的一個足球狀C60分子(由60個碳原子構成的分子，形似足球，又名足球烯)組成，所述碳納米管為層狀中空結構，管身是準圓管結構，由多個六邊形碳環結構單元組成，其直徑一般在一到幾十個納米之間，長度則遠遠大于其直徑，對該模型進行初始化，讀取模型中所有粒子的速度和位置坐標。

圖1 碳納米管

圖2 C60分子

圖3 碳納米管系統模型

在碳納米管中，除管口層粒子外，其余粒子相互組成正六邊形。如圖4所示，非管口粒子都有3個近鄰粒子(ip，jp，kp均為idx粒子的近鄰粒子)。

圖4 碳納米結構

在碳納米管系統中，因為C60分子受到碳納米管管壁上碳粒子產生的力，在平衡位置上來回振動，當C60分子處于平衡狀態時，采用Gaussian熱浴法(約束溫度調節方法)，在碳納米管兩端同時分別設置不同溫度的熱浴條件，并設定熱浴時間，使碳納米管管壁上碳粒子和C60分子因平衡狀態被打破而發生位置和速度的變化。而這些變化是相互影響的，具有傳遞性，伴隨著這些粒子位置和速度的變化，能量在長管狀的碳納米管系統內發生傳遞。由于C60分子位于碳納米管內，其受力會隨著碳納米管組成粒子的位置改變而改變，因此伴隨著能量傳遞的過程，C60分子的運動狀態和能量也會發生改變;采用Gaussian熱浴法的基本原理是在運動方程中加入摩擦力fi，并將其與粒子速度vi聯系起來，其受力公式為:

當平衡態時，系統溫度不變，因此有dEk/dt=0。即有:，由此可得

3．2 并行仿真算法研究

在CUDA平臺上對碳納米管進行分割，分割成多層大小適宜、相互獨立的計算單元，分割原則是在避免過多重復計算的前提上提高并行度，即分割的計算單元必須合理粗細化。因為計算單元過大，則并行不明顯，計算單元過小，會導致太多不必要的重復計算。由于粒子間的作用力間距影響較大，必須設定一個距離臨界值，以判斷其位置關系，當粒子間的間距小于距離臨界值，為近鄰粒子;當粒子間的間距大于距離臨界值，為次鄰粒子，則認為其相互間的分子作用力可以忽略不計。

可以預料到的是每個中心粒子最多只會對其3個近鄰及3個近鄰的各2個近鄰，即6個次鄰施力?？梢赃@樣來考慮，以中心粒子以及3個近鄰6個次鄰劃分到一組計算。分批次并行計算。在某一批次計算中，未當成中心粒子計算的一個粒子加入可并行計算隊列中，由于競爭關系，其3個近鄰6個次鄰在本批次中將不能加入隊列，同樣的，對這3個近鄰6個次鄰施力的其他所有未計算粒子也不能加入隊列中。這樣，同一批次的粒子將可以并行計算而不產生競爭關系。

設計分裂算法，采用CPU遍歷計算單元得到n個可并行運算隊列:

(1)若所有粒子均已作為中心計算粒子，跳到步驟(8)。

(2)按次序找到第一個非競爭粒子并沒有按中心計算的粒子，加入可并行運算隊列。

(3)標記該粒子的所有近鄰為本并行隊列一度不可并行粒子，若近鄰已經是二度不可并行粒子，則提升為一度不可計算。

(4)標記該粒子所有次鄰為本并行隊列二度不可并行粒子，若次鄰為一度不可并行粒子，則不修改其度數。

(5)標記本粒子為已按中心計算粒子。

(6)判斷是否已遍歷到粒子隊尾。若是繼續執行，若否跳到步驟(1)。

(7)下一個可并行隊列開始，返回步驟(2)。

(8)結束。

分裂算法在CPU預處理階段完成，只計算一次得出的可并行隊列可以在后面的循環計算中一直使用。當循環次數較多時，分裂算法所使用的時間占總時間比例將大大減少，對程序總體性能影響也降到比較低的程度。但是由于算法核心是批次計算，因此計算力時將多次啟動內核計算函數，一定程度上加大了運行開銷。

如圖5和圖6所示，調度GPU的流處理單元，采用Verlet算法進行并行運算和處理:

(1)按碳納米管系統模型所有粒子的位置，計算近鄰粒子以及次鄰粒子間的鍵關系和角度關系。

(2)調度GPU的流處理單元，并行計算被分割在不同計算單元里的粒子，并積累計算碳納米管管壁上每個粒子與其近鄰粒子的相互作用力。

(3)根據C60分子所在區域，計算C60分子內每個粒子與其所在碳納米管區域中的每個粒子的相互作用力(范德華力)。

(4)根據粒子所受到的力及其速度，更新粒子位置，再執行步驟(2)和步驟(3)。

(5)根據粒子所受到的力，計算粒子本次速度及碳納米管前、后兩端熱浴的熱流值。

(6)若達到循環頻數，則計算結束;否則，在CPU端間隔保存粒子的數據，返回步驟(4)。

圖5 碳納米管分子動力學并行仿真中的調度模式

圖6 碳納米管分子動力學并行仿真算法的流程

4 GPU并行算法的優化

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

4．1 指令優化

目前，GPU單精度計算性能要遠遠超過雙精度計算性能，因此，程序中對精度要求不高的部分可使用單精度運算代替雙精度運算。在GPU中整數乘法、除法、求模運算的指令吞吐量較為有限，而控制流指令由于會影響指令發射器的并行發射，打斷指令流開銷巨大。此時，在可知整數范圍中以經優化的庫24位整數運算__mul24，__umul24代替32位整數運算，以__sinf(x)，__fdivide(x，y)等代替相應的運算。盡量避免整數除和模運算，或以(i＆(n－1))代替(i%n)，除數是2的冪次方時以移位運算代替除法。

4．2 存儲器訪問優化

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

4．3 CPU+GPU 處理

考慮到GPU并行處理要發揮其強大的并行計算能力，要有足夠多的線程并行運算。因此，把部分較小計算量的計算交給CPU執行。如計算C60與炭納米管間范德華力函數，由于C60只有60個粒子，相對來說并行度小，將其計算放到CPU里。由于內核的啟動是異步的，因此把CPU函數放到內核啟動函數下面，CPU啟動內核執行GPU運算時可立即返回執行CPU程序，實現CPU與GPU的并行處理。

5 實驗結果與分析

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

表1 測試平臺部分硬件參數

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

表2 串行并行程序運行于不同硬件平臺時的數據

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

為使數據形象化，將GPU1與GPU2數據繪制成條狀圖。由圖7可得，前期粒子數較少時，2種顯卡耗時相差不大，但隨著粒子數的增加，GeForce GT 430顯卡耗時增長速度逐漸快于Quadro 2 000顯卡。而當粒子數量增加到120 000個時，前者耗時幾乎是后者的2倍。

圖7 2種并行顯卡運行程序耗時對比

由表1可知，Quadro 2 000有著2倍于GeForce GT 430的 CUDA Cores(即流處理器 stream processor)。結合GPU的線程運行方式，當GPU的處理單元增多，同時處于運行狀態的線程也會相應增加。因此，當無其他硬件限制時，GPU并行計算的性能幾乎隨著GPU處理器的增加而線性增加，使用擁有更多處理器的GPU能給Verlet帶來更多的性能提升。

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

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

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

圖8 串并行程序執行時間對比

5．3 并行算法性能分析

為分析并行算法的加速比，將表2中串并行程序運算時間按公式:

計算，GPU使用Quadro 2000，并將結果繪制成GPU并行計算加速比，如圖9所示。

圖9 GPU并行計算加速比

可以看出，隨著粒子數的增加，GPU并行算法相比于串行算法的加速比也隨之增大。當粒子數到達120 000個時，加速比達到1 400%，即14倍的加速效果。

為分析影響GPU并行程序性能因素，使用性能評估工具NVIDIA Visual Profiler對程序進行分析。目標程序計算600層12 000個粒子，每個block中有128個線程。首先分析各函數運行時間占總運算時間的比例。通過多次運算求平均值，如圖10所示。

圖10 各個函數的比重

6 結束語

在對碳納米管分子動力學進行仿真時，原來的分子動力學仿真算法有良好的效果，但因碳納米管系統粒子數的規模太大、計算量多、運行時間長，缺乏實用性。而由用Verlet算法實現的分子動力學仿真算法具有易并行性，可以在具有強大圖像處理和浮點計算能力的GPU上并行處理，從而大大提高了碳納米管分子動力學仿真算法的效率。

實驗結果表明，對小規模僅有120 000個粒子組成的碳納米管系統，在擁有192個流處理器的NVIDIA Quadro 2000顯卡上實現的GPU并行仿真算法，與在擁有4核8線程的Intel Xeon W3550上實現的基于CPU的串行仿真算法相比，其加速比即可達到十多倍，從而有效地解決了算法速度上的缺陷。

然而由于本次實驗的硬件條件有限，所能計算的碳納米管系統規模有限，因此無法深度發掘CUDA的能力。而由實驗結果所示，碳納米管系統的粒子數規模越大，GPU并行算法對于串行算法的加速比就越明顯。在實際運用中，碳納米管系統的粒子數的規模都很大，超過千萬乃至十億，可以預估該基于GPU的分子動力學并行仿真算法，對碳納米管的研究會發揮至關重要的作用。

［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］ 唐玉蘭，胡 適，王東旭，等．納米工程中大規模分子動力學仿真算法的研究進展［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算法在集群系統中的并行實現研究［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］ 孟小華，劉堅強，區業祥，等．基于CUDA的拉普拉斯邊緣檢測算法［J］．計算機工程，2012，38(18):190-193．

［13］ 孟小華，黃叢珊，朱麗莎．基于CUDA的熱傳導GPU并行算法研究［J］．計算機工程，2014，40(5):41-44，48．