基于GPU 的固態(tài)晶體硅分子動(dòng)力學(xué)算法優(yōu)化

李 靖，祝愛(ài)琦，韓 林，侯超峰

（1.鄭州大學(xué) 信息工程學(xué)院，鄭州 450001；2.中國(guó)科學(xué)院過(guò)程工程研究所，北京 100190；3.鄭州大學(xué) 國(guó)家超級(jí)計(jì)算鄭州中心，鄭州 450001）

0 概述

近年來(lái)，利用圖形處理器（Graphics Processing Unit，GPU）加速分子動(dòng)力學(xué)（Molecular Dynamics，MD）模擬已成為計(jì)算化學(xué)、計(jì)算材料學(xué)等領(lǐng)域的研究熱點(diǎn)［1-2］。分子動(dòng)力學(xué)模擬通常從原子、分子尺度上研究凝聚態(tài)材料的熱力學(xué)性質(zhì)和行為，被廣泛應(yīng)用于微觀計(jì)算模擬領(lǐng)域中［3］。晶體硅是信息技術(shù)產(chǎn)業(yè)中應(yīng)用最廣泛和最重要的材料之一，其納微結(jié)構(gòu)的熱傳遞性質(zhì)和導(dǎo)熱性能會(huì)顯著影響電子設(shè)備的穩(wěn)定性，導(dǎo)致器件性能降低。因此，分子動(dòng)力學(xué)模擬方法對(duì)接近真實(shí)尺度硅納微結(jié)構(gòu)的直接模擬的實(shí)現(xiàn)與其導(dǎo)熱機(jī)制的研究具有重要意義。

模擬體系的計(jì)算量和時(shí)空分辨率的提高是分子動(dòng)力學(xué)模擬領(lǐng)域的研究熱點(diǎn)。高性能加速芯片（如GPU）為處理模擬體系的巨量計(jì)算提供契機(jī)。GPU作為一種具有良好并行性能且適應(yīng)密集計(jì)算的新型加速計(jì)算工具，在通用科學(xué)計(jì)算領(lǐng)域具有較大的應(yīng)用潛力。分子動(dòng)力學(xué)模擬優(yōu)化的研究大多集中在通用分子動(dòng)力學(xué)模擬算法優(yōu)化和程序并行算法設(shè)計(jì)方面，很少專(zhuān)門(mén)針對(duì)固體材料并結(jié)合GPU 硬件結(jié)構(gòu)特性對(duì)MD 算法進(jìn)行設(shè)計(jì)和優(yōu)化。文獻(xiàn)［4］提出一種高效的多粒子碰撞動(dòng)力學(xué)算法，并在多個(gè)GPU 上實(shí)現(xiàn)。針對(duì)多體勢(shì)的高效分子動(dòng)力學(xué)模擬，文獻(xiàn)［5］提出一種基于GPU 的力求解算法。文獻(xiàn)［6］實(shí)現(xiàn)分子動(dòng)力學(xué)三體勢(shì)的GPU 加速計(jì)算。文獻(xiàn)［7］提出不同于對(duì)勢(shì)分子動(dòng)力學(xué)模擬的方法，根據(jù)復(fù)雜的多體勢(shì)描述一種GPU 加速的固態(tài)共價(jià)晶體分子動(dòng)力學(xué)模擬的高效和可擴(kuò)展算法。

針對(duì)固態(tài)共價(jià)晶體硅的分子動(dòng)力學(xué)模擬，本文結(jié)合原子間多體勢(shì)函數(shù)與CUDA 編程模型，通過(guò)對(duì)GPU 并行計(jì)算特征和模擬熱點(diǎn)進(jìn)行分析，提出面向GPU 計(jì)算平臺(tái)的固定鄰居算法設(shè)計(jì)與優(yōu)化。根據(jù)GPU 硬件架構(gòu)特性，將并行線(xiàn)程處理的float4 數(shù)據(jù)結(jié)構(gòu)重新映射為連續(xù)float 布置的形式，采用三目運(yùn)算符控制指令排隊(duì)，消除條件分支語(yǔ)句造成的流水線(xiàn)停頓，減少全局訪存消耗和分支結(jié)構(gòu)耗時(shí)，提升算法的計(jì)算效率。

1 Tesla V100 GPU 硬件架構(gòu)

為充分發(fā)揮GPU 的計(jì)算性能，從軟件與硬件相結(jié)合的角度分析固態(tài)材料分子動(dòng)力學(xué)模擬算法。以Tesla V100 GPU 為實(shí)驗(yàn)平臺(tái)描述GPU 硬件架構(gòu)［8-10］。Tesla V100 GPU 架構(gòu)如圖1 所示，由多個(gè)圖形處理簇（Graphics Processing Cluster，GPC）、紋理處理簇（Texture Processing Cluster，TPC）、流式多處理器（Streaming Multiprocessor，SM）以及內(nèi)存控制器組成。每個(gè)SM 內(nèi)有64 個(gè)流處理器（Streaming Processor，SP），其中，每個(gè)SP 具有自己的緩存空間和寄存器組，并通過(guò)線(xiàn)程束（warp）調(diào)度器的方式控制SP 的并行執(zhí)行，將L1 高速緩存和共享內(nèi)存組合到單個(gè)內(nèi)存塊中，內(nèi)存總?cè)萘繛?28 KB，在SM 上運(yùn)行的線(xiàn)程對(duì)L1 高速緩存和共享內(nèi)存的訪問(wèn)相對(duì)于全局內(nèi)存具有更快的訪問(wèn)速度，在不同SM 之間的L1 高速數(shù)據(jù)緩存和共享內(nèi)存相互獨(dú)立。Tesla V100 GPU 總共擁有80 個(gè)SM、5 120 個(gè)32 位浮點(diǎn)處理單元、2 560 個(gè)64 位浮點(diǎn)處理單元、640 個(gè)Tensor 核心以及320 個(gè)紋理單元（Tex）［8］。

圖1 Tesla V100 GPU 架構(gòu)Fig.1 Architecture of Tesla V100 GPU

在硬件與執(zhí)行程序的對(duì)應(yīng)關(guān)系上，一個(gè)SP 控制一個(gè)線(xiàn)程的執(zhí)行。為提高并行效率，NVIDIA 把32 個(gè)線(xiàn)程組織成線(xiàn)程束的形式，并在warp 內(nèi)以單指令多線(xiàn)程（Single Instruction Multiple Threads，SIMT）的方式執(zhí)行。在線(xiàn)程塊中的多個(gè)warp 由線(xiàn)程束調(diào)度器調(diào)度進(jìn)入SM 執(zhí)行［11］。在Tesla V100 GPU 中大量的SM 雖然提高GPU 硬件的并行性，但是其具有較高并行性的復(fù)雜架構(gòu)給MD 模擬算法的優(yōu)化帶來(lái)挑戰(zhàn)。在線(xiàn)程并行優(yōu)化方面，不僅減少線(xiàn)程之間的相互依賴(lài)關(guān)系，還需要充分利用其硬件特性以減少訪存開(kāi)銷(xiāo)。

2 固定鄰居算法在GPU 上的多體勢(shì)實(shí)現(xiàn)

2.1 Tersoff 多體勢(shì)的分子動(dòng)力學(xué)模擬

MD 模擬是通過(guò)刻畫(huà)微觀原子的運(yùn)動(dòng)過(guò)程，由相關(guān)的數(shù)據(jù)統(tǒng)計(jì)來(lái)研究模擬體系的結(jié)構(gòu)和熱力學(xué)性質(zhì)的方法。其中，在周?chē)従釉觿?shì)場(chǎng)作用下每個(gè)原子根據(jù)牛頓定律運(yùn)動(dòng)［12］。原子間的相互作用是MD 模擬算法的主要組成部分。與原子間對(duì)勢(shì)相比，原子間多體勢(shì)具有更復(fù)雜的模型形式，但通常能更準(zhǔn)確、更有效地體現(xiàn)共價(jià)材料模擬體系的物理特性［13］。

本文針對(duì)固態(tài)共價(jià)晶體硅，基于固定鄰居算法，利用Tersoff 多體勢(shì)實(shí)現(xiàn)晶體硅的MD 模擬計(jì)算。Tersoff 勢(shì)函數(shù)的數(shù)學(xué)表達(dá)式如式（1）和式（2）所示：

其中：UR為排斥勢(shì)項(xiàng)；UA為吸引勢(shì)項(xiàng)；fC為截?cái)嗪瘮?shù)；bij為鍵級(jí)項(xiàng)。UR、UA、fC、bij、ζij和g（θ）的數(shù)學(xué)表達(dá)式如式（3）～式（8）所示：

其中：ζij為原子i與原子j間的角勢(shì)能。上述公式的模型參數(shù)見(jiàn)文獻(xiàn)［14-15］。固態(tài)共價(jià)晶體硅的MD 模擬計(jì)算首先需要生成模擬體系晶體硅原子的初始位置，并在特定溫度下給出原子的初始速度，通過(guò)Tersoff 多體勢(shì)計(jì)算每個(gè)硅原子受其他硅原子的作用合力，然后選擇蛙跳算法積分牛頓運(yùn)動(dòng)方程，求解硅原子的新速度和新位置。

2.2 固定鄰居算法的GPU 實(shí)現(xiàn)

固態(tài)共價(jià)晶體硅具有規(guī)則的晶胞結(jié)構(gòu)，在低溫狀態(tài)下每個(gè)原子周?chē)? 個(gè)鄰居原子。在Tesla V100 GPU 平臺(tái)實(shí)現(xiàn)算法時(shí)，模擬體系中所有的硅原子按超胞［7］的順序在分配的設(shè)備存儲(chǔ)器上排列，并且在后續(xù)的迭代模擬中按超胞順序分配的內(nèi)存位置保持不變。超胞中原子的分組情況如圖2 所示。根據(jù)每個(gè)原子的鄰居原子在超胞內(nèi)的數(shù)目，把超胞內(nèi)的原子分為三組：第一組原子在超胞內(nèi)有一個(gè)鄰居原子，如圖2（a）中的灰色原子；第二組原子在超胞內(nèi)有兩個(gè)鄰居原子；第三組原子在超胞內(nèi)有四個(gè)鄰居原子，如圖2（a）中的黑色原子。算法執(zhí)行時(shí)將超胞內(nèi)上述三組原子按順序布置，如圖2（b）所示。此外，為提高緩存效率，采用靜態(tài)重排方法代替動(dòng)態(tài)更新鄰居列表的過(guò)程，根據(jù)超胞內(nèi)原子分組情況，把原子在超胞內(nèi)和超胞外的鄰居原子信息分組聚合并進(jìn)行連續(xù)排列。由于固定鄰居算法沒(méi)有鄰居列表動(dòng)態(tài)更新的過(guò)程，因此原子間作用力的計(jì)算、原子位置和速度的更新占大量的模擬總時(shí)間。

圖2 每個(gè)超胞中計(jì)算原子的分組示意圖Fig.2 Schematic diagram of grouping of calculated atoms in each supercell

固定鄰居算法流程如圖3 所示，在CPU 端通過(guò)控制模擬的原子數(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ù)“一個(gè)線(xiàn)程模擬一個(gè)原子”的分配原則，使超胞中的原子數(shù)等于線(xiàn)程塊中的線(xiàn)程數(shù)，并利用線(xiàn)程分支控制超胞內(nèi)的分組情況。在GPU設(shè)備端并行計(jì)算每個(gè)超胞中每個(gè)原子和其鄰居原子之間的距離，進(jìn)一步計(jì)算每個(gè)原子受鄰居原子的全部作用力。最后，通過(guò)積分時(shí)間計(jì)算原子的新速度和新位置。當(dāng)算法沒(méi)有達(dá)到模擬設(shè)定的迭代步數(shù)時(shí)，GPU 端繼續(xù)執(zhí)行MD 模擬，直到達(dá)到設(shè)定的模擬迭代步數(shù)，程序結(jié)束。

圖3 固定鄰居算法流程Fig.3 Procedure of fixed neighbor algorithm

雖然固定鄰居算法采用超胞分組的方式提高緩存命中率，但是在設(shè)備端并行執(zhí)行原子受力計(jì)算時(shí)也增加了if 分支，為GPU 程序并行執(zhí)行帶來(lái)了額外的分支消耗。Tersoff 多體勢(shì)比對(duì)勢(shì)具有更復(fù)雜的模型形式，當(dāng)計(jì)算原子間作用力時(shí)使用數(shù)據(jù)類(lèi)型（float4）并未實(shí)現(xiàn)計(jì)算吞吐量和數(shù)據(jù)全局訪存的最佳匹配。因此，為提高程序的性能，本文結(jié)合GPU硬件架構(gòu)特性，對(duì)固定鄰居算法進(jìn)行設(shè)計(jì)優(yōu)化。

3 固定鄰居算法優(yōu)化實(shí)現(xiàn)

在對(duì)Tesla V100 GPU 平臺(tái)的硬件特性進(jìn)行分析之后，本文分別對(duì)固定鄰居算法的數(shù)據(jù)結(jié)構(gòu)和分支結(jié)構(gòu)進(jìn)行優(yōu)化，充分利用GPU 計(jì)算資源來(lái)提高線(xiàn)程并行效率，進(jìn)一步提升固定鄰居算法的性能。

3.1 數(shù)據(jù)結(jié)構(gòu)優(yōu)化

早期的GPU 適用于圖形渲染設(shè)計(jì)，其運(yùn)算對(duì)象是表示顏色的R、G、B、A 和表示坐標(biāo)的X、Y、Z、W，都是對(duì)一個(gè)四維向量進(jìn)行操作［16］。隨著GPU 應(yīng)用越來(lái)越廣泛，數(shù)據(jù)結(jié)構(gòu)可能不僅具有表示坐標(biāo)的X、Y、Z、W和表示顏色的R、G、B、A，還有可能存在一些非向量的數(shù)據(jù)。因此，本文對(duì)不同的數(shù)據(jù)做不同的處理，而不同的數(shù)據(jù)結(jié)構(gòu)方式會(huì)影響數(shù)據(jù)的內(nèi)存局部性。為提高通用性，NVIDIA 引入標(biāo)量計(jì)算的概念，當(dāng)發(fā)展到Tesla V100 時(shí)，其內(nèi)部包含大量專(zhuān)門(mén)針對(duì)標(biāo)量的計(jì)算單元（CUDA CORE），大幅提高運(yùn)算性能。

在固定鄰居算法中，數(shù)據(jù)結(jié)構(gòu)主要包括各個(gè)原子的位置、速度、力等屬性，在三維空間中，它們都是以四維向量方式定義的，四個(gè)分量x、y、z、w分別表示在三個(gè)方向上的方位坐標(biāo)和空值。線(xiàn)程對(duì)數(shù)據(jù)結(jié)構(gòu)的重映射示意圖如圖4 所示。float4 數(shù)據(jù)結(jié)構(gòu)如圖4（a）所示，原子屬性以四維向量x、y、z、w的形式連續(xù)存儲(chǔ)在全局內(nèi)存中，從三個(gè)維度分別計(jì)算原子間的距離和受力。當(dāng)計(jì)算x方向上兩個(gè)原子間的距離和受力時(shí)，在warp 中，如Thread 0 訪問(wèn)x1，Thread 1訪問(wèn)x2，x1與x2之間相隔y1、z1、w1，兩個(gè)線(xiàn)程對(duì)原子數(shù)據(jù)的訪問(wèn)和計(jì)算會(huì)造成計(jì)算吞吐量的降低。在圖4（b）中，結(jié)合Tesla V100 GPU 的硬件特性，修改算法數(shù)據(jù)結(jié)構(gòu)布局，使原子屬性以單位標(biāo)量x、y、z、w的形式各自分開(kāi)存儲(chǔ)在全局內(nèi)存中，在計(jì)算原子間x方向上的距離和受力時(shí)，在warp 中，Thread 0 訪問(wèn)x1，Thread 1 訪問(wèn)x2，類(lèi)似操作實(shí)現(xiàn)了線(xiàn)程對(duì)原子數(shù)據(jù)的連續(xù)訪問(wèn)。這種新的數(shù)據(jù)結(jié)構(gòu)布局方式有助于充分利用硬件特性，提升計(jì)算吞吐量和訪存帶寬性能，提高算法的運(yùn)行效率。

圖4 線(xiàn)程對(duì)數(shù)據(jù)結(jié)構(gòu)的重映射示意圖Fig.4 Remapping schematic diagram of data structures by threads

3.2 分支結(jié)構(gòu)優(yōu)化

GPU 分支結(jié)構(gòu)執(zhí)行方式不同于CPU 分支結(jié)構(gòu)執(zhí)行方式。在CPU 中設(shè)置分支預(yù)測(cè)機(jī)制，有效縮短執(zhí)行時(shí)間。在GPU 執(zhí)行中，如果warp 中任何一個(gè)線(xiàn)程進(jìn)入到不同的分支路徑，則其余的線(xiàn)程都處于等待狀態(tài)，直到該線(xiàn)程執(zhí)行完分支程序。warp 中線(xiàn)程按順序串行通過(guò)多條不同分支路徑，當(dāng)所有分支路徑都執(zhí)行后，warp 中的所有線(xiàn)程才會(huì)回到同一條執(zhí)行路徑上。

warp 分支結(jié)構(gòu)如圖5 所示。一個(gè)線(xiàn)程束warp 包含32 個(gè)線(xiàn)程，在warp 內(nèi)以SIMT 的方式執(zhí)行［17］。在執(zhí)行if 分支時(shí)，當(dāng)前18 個(gè)線(xiàn)程指向的數(shù)據(jù)執(zhí)行時(shí)，剩余后14 個(gè)線(xiàn)程指向的數(shù)據(jù)都需要等待，造成線(xiàn)程資源浪費(fèi)。若程序中的分支路徑較多，且warp 中的線(xiàn)程都需要執(zhí)行時(shí)，會(huì)延長(zhǎng)程序的執(zhí)行時(shí)間，造成性能的嚴(yán)重?fù)p失。

圖5 warp 分支結(jié)構(gòu)Fig.5 Structure of warp branch

大多數(shù)對(duì)于分支結(jié)構(gòu)的優(yōu)化保留在線(xiàn)程層，線(xiàn)程交換是通過(guò)減少同一個(gè)warp 中不同分支的方式來(lái)提高程序并行性，根據(jù)分支控制數(shù)據(jù)的依賴(lài)關(guān)系，重新排列數(shù)據(jù)數(shù)組，建立線(xiàn)程與數(shù)據(jù)的重映射關(guān)系，從而提升性能［18-20］。然而，一些特定算法的分支結(jié)構(gòu)無(wú)法通過(guò)上述方式進(jìn)行優(yōu)化，本文引入三目運(yùn)算符的方法對(duì)分支結(jié)構(gòu)進(jìn)行優(yōu)化。分支結(jié)構(gòu)轉(zhuǎn)換為三目運(yùn)算符的示意圖如圖6 所示。

圖6 分支結(jié)構(gòu)轉(zhuǎn)換為三目運(yùn)算符的示意圖Fig.6 Schematic diagram of branch structure converted to the ternary operator

固定鄰居算法在超胞中的原子受力需要分組賦值的情況下，采用三目運(yùn)算符代替分支結(jié)構(gòu)。在執(zhí)行分支結(jié)構(gòu)代碼時(shí)，if-else 語(yǔ)句會(huì)對(duì)分支指令進(jìn)行排隊(duì)處理，而三目運(yùn)算符對(duì)分支的兩種情況同時(shí)進(jìn)行預(yù)處理，使三目運(yùn)算符代碼的耗時(shí)周期相較于分支結(jié)構(gòu)代碼更短，有效減少分支耗時(shí)并提升算法的運(yùn)行效率。

4 計(jì)算驗(yàn)證與分析

4.1 計(jì)算環(huán)境

本文固態(tài)共價(jià)晶體硅的MD 模擬是在Tesla V100 GPU 單卡上實(shí)現(xiàn)并行執(zhí)行。固定鄰居算法的應(yīng)用環(huán)境如表1 所示。采用算法性能最優(yōu)的超胞尺寸，其包含256 個(gè)硅原子，在x、y、z方向上長(zhǎng)度分別為2、4、4 個(gè)單位晶胞，模擬的初始溫度為300 K。

表1 固定鄰居算法的計(jì)算環(huán)境Table 1 Computational environment of fixed neighbor algorithm

4.2 計(jì)算測(cè)試分析

4.2.1 正確性分析

本文通過(guò)算例驗(yàn)證固定鄰居算法優(yōu)化前后的有效性和正確性，總的迭代步數(shù)設(shè)置為100 萬(wàn)步，原子模擬數(shù)目為32 768 個(gè)，每間隔1 000 迭代步，輸出一次統(tǒng)計(jì)結(jié)果。GPU 固定鄰居算法優(yōu)化前后每個(gè)原子的平均勢(shì)能隨時(shí)間變化曲線(xiàn)如圖7 所示，SP 表示單精度，DP 表示雙精度。

圖7 GPU 固定鄰居算法優(yōu)化前后每個(gè)原子的平均勢(shì)能隨時(shí)間變化曲線(xiàn)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)化前后原子的平均勢(shì)能變化穩(wěn)定，雙精度的原子勢(shì)能趨向于-4.610 24 eV，單精度的原子勢(shì)能趨向于-4.610 20 eV，模擬結(jié)果與文獻(xiàn)［7］一致，證明了固定鄰居算法優(yōu)化前后的計(jì)算正確可靠。

4.2.2 性能測(cè)試與分析

在超胞尺寸固定不變的情況下，本文考慮到Tesla V100 GPU 單卡硬件資源的限制，分別對(duì)不同原子數(shù)目進(jìn)行模擬，實(shí)驗(yàn)選擇原子數(shù)目為32 768、262 144、1 024 000、4 096 000、10 240 000，以模擬1 000 迭代步的運(yùn)行時(shí)間作為評(píng)價(jià)算法性能的指標(biāo)，通過(guò)對(duì)比算法優(yōu)化前后的耗時(shí)來(lái)驗(yàn)證優(yōu)化效果。

將數(shù)據(jù)結(jié)構(gòu)與分支結(jié)構(gòu)優(yōu)化方法相結(jié)合，固定鄰居算法優(yōu)化前后的性能對(duì)比如表2 所示。t1表示雙精度固定鄰居算法優(yōu)化前的運(yùn)行時(shí)間，t2表示雙精度固定鄰居算法優(yōu)化后的運(yùn)行時(shí)間，t3表示單精度固定鄰居算法優(yōu)化前的運(yùn)行時(shí)間，t4表示單精度固定鄰居算法優(yōu)化后的運(yùn)行時(shí)間。當(dāng)原子數(shù)目為4 096 000時(shí)，優(yōu)化后雙精度固定鄰居算法的性能最優(yōu)，與固定鄰居算法優(yōu)化前的運(yùn)行時(shí)間相比，雙精度固定鄰居算法與單精度固定鄰居算法的加速效果分別提高了27.7%與20.7%。

表2 固定鄰居算法優(yōu)化前后性能對(duì)比Table 2 Performance comparison of fixed neighbor algorithms before and after optimization

固定鄰居算法優(yōu)化后的加速比如圖8 所示。當(dāng)原子數(shù)目為32 768 時(shí)，優(yōu)化后雙精度固定鄰居算法加速比是優(yōu)化后單精度固定鄰居算法加速比的1 倍多。其原因?yàn)楫?dāng)原子規(guī)模較小時(shí)，優(yōu)化后雙精度固定鄰居算法計(jì)算1 000 迭代步的運(yùn)行時(shí)間較短。此時(shí)，單精度固定鄰居算法計(jì)算力的核函數(shù)耗時(shí)占GPU 端總耗時(shí)的比例較小，雙精度固定鄰居算法計(jì)算力的核函數(shù)耗時(shí)占GPU 端總耗時(shí)比例較大，因此，優(yōu)化后雙精度固定鄰居算法的加速比較高。固定鄰居算法從原子數(shù)目為262 144 開(kāi)始，優(yōu)化后單精度和雙精度固定鄰居算法的加速比隨著原子數(shù)目的增大而增大，并在硬件資源限制的條件下達(dá)到一個(gè)最高值，然后，加速比開(kāi)始下降。

圖8 固定鄰居算法優(yōu)化后的加速比Fig.8 Acceleration ratio of fixed neighbor algorithm after optimization

在Tesla V100 GPU 處理器資源可承受范圍內(nèi)，固定鄰居算法的優(yōu)化效果依賴(lài)于計(jì)算的浮點(diǎn)精度和計(jì)算規(guī)模。在GPU 計(jì)算資源和存儲(chǔ)資源相同的條件下，原子數(shù)目多的體系優(yōu)化效果更好。其原因?yàn)楫?dāng)原子數(shù)目較小時(shí)，浮點(diǎn)計(jì)算不能充分利用Tesla V100 GPU 處理器并行資源，計(jì)算性能不能達(dá)到最優(yōu)，當(dāng)原子數(shù)目較大時(shí)，Tesla V100 GPU 處理器才能充分發(fā)揮計(jì)算優(yōu)勢(shì)，加速效果明顯。

4.3 與LAMMPS和HOOMD-blue軟件的性能對(duì)比

LAMMPS 是一種支持多種原子間勢(shì)函數(shù)并被廣泛應(yīng)用的開(kāi)源MD 軟件，利用空間分解技術(shù)將三維模擬區(qū)域劃分為更小的三維子域，利用MPI 并行化提高計(jì)算規(guī)模和性能。KOKKOS 加速包是基于LAMMPS 的C++加速庫(kù)，可以有效地在GPU 硬件上運(yùn)行［21］。目前，KOKKOS 加速包沒(méi)有精度選項(xiàng)，所有編譯和計(jì)算都必須以雙精度執(zhí)行［22-23］。HOOMDblue 是專(zhuān)門(mén)支持GPU 加速的分子動(dòng)力學(xué)模擬軟件，能夠在多種力場(chǎng)下執(zhí)行MD 模擬計(jì)算［24-25］。本文使用LAMMPS 軟件KOKKOS 加速包和支持GPU 加速的HOOMD-blue 軟件在Tesla V100 GPU 上對(duì)Tersoff多體勢(shì)進(jìn)行單雙精度計(jì)算，測(cè)試環(huán)境與4.1 節(jié)相同，選擇模擬1 000 迭代步的運(yùn)行時(shí)間作為性能對(duì)比的基準(zhǔn)。

LAMMPS 和HOOMD-blue 與優(yōu)化后的固定鄰居算法性能對(duì)比如表3 所示，表中T1表示雙精度固定鄰居算法的執(zhí)行時(shí)間，T2表示單精度固定鄰居算法的執(zhí)行時(shí)間，T3表示LAMMPS 雙精度固態(tài)晶體硅分子動(dòng)力學(xué)模擬的執(zhí)行時(shí)間，T4表示HOOMD-blue 雙精度固態(tài)晶體硅分子動(dòng)力學(xué)模擬的執(zhí)行時(shí)間，T5表示HOOMDblue 單精度固態(tài)晶體硅分子動(dòng)力學(xué)模擬的執(zhí)行時(shí)間，加速比1 表示LAMMPS 雙精度固態(tài)晶體硅分子動(dòng)力學(xué)模擬與雙精度固定鄰居算法（LAMMPS 雙精度/雙精度算法）的執(zhí)行時(shí)間比值，加速比2 表示HOOMD-blue雙精度固態(tài)晶體硅分子動(dòng)力學(xué)模擬與雙精度固定鄰居算法（HOOMD-blue 雙精度/雙精度算法）的執(zhí)行時(shí)間比值，加速比3 表示HOOMD-blue 單精度固態(tài)晶體硅分子動(dòng)力學(xué)模擬與單精度固定鄰居算法（HOOMD-blue單精度/單精度算法）的執(zhí)行時(shí)間比值。

表3 LAMMPS 和HOOMD-blue 與優(yōu)化后的固定鄰居算法性能對(duì)比Table 3 Performance comparison of LAMMPS and HOOMD-blue with the optimized fixed neighbor algorithm

相比兩個(gè)MD 開(kāi)源模擬軟件，固定鄰居算法具有較優(yōu)的加速效果。當(dāng)原子數(shù)目為10 240 000 時(shí)，雙精度固定鄰居算法加速效果最優(yōu)，LAMMPS 和 HOOMD-blue 軟件與雙精度固定鄰居算法的執(zhí)行時(shí)間的比值分別為 11.62 和 9.39。當(dāng)原子數(shù)目為 4 096 000 時(shí)，單精度固定鄰居算法的加速效果較優(yōu)，HOOMD-blue 軟件與單精度固定鄰居算法執(zhí)行時(shí)間的比值為12.18。優(yōu)化后的固定鄰居算法相對(duì)LAMMPS 和HOOMD-blue 的加速比如圖9 所示。隨著原子數(shù)目的增加，雙精度固定鄰居算法與兩個(gè) MD 開(kāi)源模擬軟件相比，其加速比呈現(xiàn)上升趨勢(shì)。HOOMD-blue 單精度/單精度算法加速比呈現(xiàn)先升高后略微下降的趨勢(shì)。

圖9 優(yōu)化后的固定鄰居算法相對(duì)LAMMPS 和HOOMD-blue 的加速比Fig.9 Acceleration ratio of the optimized fixed neighbor algorithm to LAMMPS and HOOMD-blue

5 結(jié)束語(yǔ)

針對(duì)固態(tài)共價(jià)晶體硅的分子動(dòng)力學(xué)模擬，本文結(jié)合GPU 的硬件架構(gòu)特性，提出分子動(dòng)力學(xué)模擬固定鄰居算法的優(yōu)化設(shè)計(jì)，采用數(shù)據(jù)結(jié)構(gòu)和分支結(jié)構(gòu)優(yōu)化方法，減少多線(xiàn)程訪存和分支耗時(shí)，提升晶體硅模擬的計(jì)算效率。與主要的GPU 加速M(fèi)D 模擬軟件LAMMPS 和HOOMD-blue 相比，本文固定鄰居算法具有較優(yōu)的加速性能。下一步將對(duì)多GPU 并行的分子動(dòng)力模擬算法進(jìn)行研究，有效提升固相晶體硅體系的計(jì)算規(guī)模，實(shí)現(xiàn)數(shù)十億原子以上的模擬。