高 萍，段曉輝，劉衛(wèi)國(guó)，郭佳旭，柳 嘉，萬(wàn)吳兵，甘 霖，楊廣文

（1.清華大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)系，北京 100084；2.國(guó)家超級(jí)計(jì)算無(wú)錫中心，江蘇 無(wú)錫 214072；3.山東大學(xué) 軟件學(xué)院，山東 濟(jì)南 250101；4.吉林大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，吉林 長(zhǎng)春 130012）

0 引言

分子動(dòng)力學(xué)（Molecular Dynamics,MD）模擬是一種計(jì)算機(jī)模擬技術(shù)，它基于經(jīng)典的牛頓力學(xué)原理，計(jì)算系統(tǒng)中粒子的受力、速度和坐標(biāo)，從而得到粒子在運(yùn)動(dòng)過(guò)程中的詳細(xì)信息，進(jìn)而獲得系統(tǒng)演化的微觀細(xì)節(jié)和系統(tǒng)宏觀量。作為原子尺度的模擬技術(shù)，MD 是一種探索微觀世界的有力工具，廣泛地應(yīng)用于物理、化學(xué)、生物和材料科學(xué)等領(lǐng)域。

近年來(lái)，二維材料的研究呈現(xiàn)指數(shù)增長(zhǎng)，成為材料科學(xué)領(lǐng)域中的研究熱點(diǎn)。安德烈·海姆因分離出備受關(guān)注的石墨烯而獲得2010 年的諾貝爾物理獎(jiǎng)。他在2019年的中國(guó)科幻大會(huì)論壇上曾說(shuō)，二維材料將成為未來(lái)人類文明的代表性材料。石墨是三維的，而二維的石墨烯只有一層原子。當(dāng)維度降低后，材料的特性和功能會(huì)發(fā)生變化，因此，揭示其物理和化學(xué)性質(zhì)是二維材料的重要研究?jī)?nèi)容之一。

然而，對(duì)于二維材料的研究，傳統(tǒng)的實(shí)驗(yàn)方法無(wú)法獲得單個(gè)原子和分子的動(dòng)態(tài)信息，而MD 模擬通過(guò)全原子力場(chǎng)計(jì)算能夠突破這一局限。目前，國(guó)內(nèi)外的研究者們已經(jīng)開(kāi)發(fā)了很多MD 軟件包，如AMBER[1]、GROMACS[2]、LAMMPS[3]、NAMD[4]等。其中，LAMMPS是材料科學(xué)領(lǐng)域最流行的開(kāi)源MD 軟件之一。它包含多種勢(shì)函數(shù)，支持不同種類材料體系的模擬，如用于模擬二維材料層間作用力的Kolmogorov-Crespi（KC）[5]力場(chǎng)、Inter-Layer Potential（ILP）[6-8]力場(chǎng)等。

當(dāng)前，我們正在邁向百億億次超級(jí)計(jì)算時(shí)代，越來(lái)越多超級(jí)計(jì)算機(jī)的處理器采用異構(gòu)眾核設(shè)計(jì)，如國(guó)產(chǎn)的神威·太湖之光[9]和新一代神威超級(jí)計(jì)算機(jī)[10]。在眾多的應(yīng)用中，MD 模擬是超級(jí)計(jì)算機(jī)上的主要數(shù)值模擬方法之一，消耗了大量的計(jì)算資源。雖然超級(jí)計(jì)算機(jī)具有強(qiáng)大的計(jì)算能力，但由于分子動(dòng)力學(xué)軟件的算法設(shè)計(jì)存在缺陷，無(wú)法充分發(fā)揮處理器的性能，導(dǎo)致計(jì)算效率偏低，難以滿足科研工作中的模擬需求。

近年來(lái)，為了提高M(jìn)D 模擬在國(guó)產(chǎn)超算上的計(jì)算效率，科研人員在異構(gòu)眾核超級(jí)計(jì)算機(jī)上進(jìn)行了一系列嘗試。在神威·太湖之光上，Dong 等人[11]對(duì)Lennard-Jones（L-J）力場(chǎng)進(jìn)行優(yōu)化，利用雙緩沖技術(shù)實(shí)現(xiàn)DMA 訪存和受力計(jì)算的重疊；Duan 等人提出軟件Cache 方案用于L-J 力場(chǎng)中的離散訪存問(wèn)題[12]，提出混合內(nèi)存更新策略用于三體作用力Tersoff 力場(chǎng)的優(yōu)化[12]；Gao 等人優(yōu)化了反應(yīng)力場(chǎng)ReaxFF，先提出隔離計(jì)算和更新策略[13]，又從數(shù)據(jù)局部性角度提出了行加鎖Cache 策略解決寫(xiě)沖突問(wèn)題[14]。除了LAMMPS，科研人員也對(duì)AMBER[15]、GROMACS[16]、NAMD[17]等MD 軟件開(kāi)展了移植和性能優(yōu)化工作。然而，已有工作中，并未有針對(duì)層間力場(chǎng)KC 的性能優(yōu)化工作。

1 新一代神威超級(jí)計(jì)算機(jī)

新一代神威超級(jí)計(jì)算機(jī)[10]繼承了神威·太湖之光[9]的體系結(jié)構(gòu)設(shè)計(jì)思想，專門(mén)為大規(guī)模線程和數(shù)據(jù)并行而設(shè)計(jì)。它的每個(gè)計(jì)算節(jié)點(diǎn)搭載的是最新的高性能申威處理器（SW39000），結(jié)構(gòu)如圖1 所示。

與申威處理器SW26010 類似，每個(gè)SW39000 處理器包含6 個(gè)核組（Core Group，CG），每個(gè)CG 包含1 個(gè)管理處理單元（通常稱為主核）和64 個(gè)計(jì)算處理單元（通常稱為從核），以及主核和從核可以共享的16 GB DDR4主存空間。

主核和從核均采用我國(guó)自主研發(fā)的申威64 指令集。主核擁有32 KB 的L1 指令緩存、32 KB 的L1 數(shù)據(jù)高速緩存和512 KB 的L2 高速緩存。從核有32 KB 的L1 指令緩存和256 KB 的本地?cái)?shù)據(jù)存儲(chǔ)（Local Data Memory，LDM）空間。其中，LDM 的部分空間可以配置成由硬件支持的自動(dòng)Cache，這一特性提升了內(nèi)存訪問(wèn)性能，是SW39000 處理器相對(duì)于SW26010 的改進(jìn)設(shè)計(jì)。主存和LDM 之間可以通過(guò)直接內(nèi)存訪問(wèn)（Direct Memory Access，DMA）完成連續(xù)數(shù)據(jù)塊傳輸，也可以通過(guò)全局讀?。℅lobal Load，GLD）和全局存儲(chǔ)（Global Store，GST）指令傳輸小而離散的數(shù)據(jù)。同一核組內(nèi)的每?jī)蓚€(gè)從核之間可以通過(guò)遠(yuǎn)程內(nèi)存訪問(wèn)（Remote Memory Access，RMA）實(shí)現(xiàn)數(shù)據(jù)傳輸。

2 層間作用力KC 力場(chǎng)

最初，KC 力場(chǎng)是由Kolmogorov 和Crespi 為石墨系統(tǒng)開(kāi)發(fā)的。由于該力場(chǎng)的勢(shì)函數(shù)能夠捕捉到層狀結(jié)構(gòu)的各向異性性質(zhì)，具有描述層間力學(xué)和摩擦學(xué)特性的能力，已被用于碳納米管和層狀石墨烯結(jié)構(gòu)的MD 研究中[18-20]。

KC 的勢(shì)函數(shù)公式由三部分構(gòu)成：長(zhǎng)范圍各向同性的L-J 吸引項(xiàng)、短范圍各向同性的類似Morse 勢(shì)的排斥項(xiàng)和用于修正各項(xiàng)異性效應(yīng)的橫向的高斯型斥力項(xiàng)。表達(dá)式如下：

其中，z0，A，λ，δ和C2n是擬合的參數(shù)。rij和ρij分別是中心原子i和j之間的全距離和橫向距離。橫向距離是指一個(gè)石墨烯層上的原子j與相鄰石墨烯層上的原子i的表面法線之間的最短距離。

由于KC 力場(chǎng)同時(shí)涉及長(zhǎng)范圍和短范圍的相互作用，且截?cái)喟霃皆介L(zhǎng)，中心原子的鄰居原子越多，計(jì)算密度越大，導(dǎo)致利用KC 力場(chǎng)在模擬二維材料體系時(shí)占用了絕大多數(shù)的計(jì)算時(shí)間，因此，優(yōu)化KC 力場(chǎng)的性能將會(huì)有效提升二維材料的模擬效率。

3 算法設(shè)計(jì)與優(yōu)化

在二維材料的MD 模擬過(guò)程中，層間作用力的計(jì)算時(shí)間占比通常在95%以上，因此屬于計(jì)算密集型問(wèn)題。而采用KC 力場(chǎng)描述層間作用力（模擬流程如圖2 所示）時(shí)，時(shí)間占比達(dá)到97%以上，因此，KC 力場(chǎng)的性能直接決定了整體的模擬效率。

圖2 KC 力場(chǎng)計(jì)算流程

本文利用熱點(diǎn)采樣工具GPTL 收集到KC 力場(chǎng)的熱點(diǎn)函數(shù)占比如表1 所示，其中，斥力計(jì)算（calc_FRep）占比98%，這是因?yàn)槌饬τ?jì)算公式更加復(fù)雜，需要計(jì)算中心原子與短鄰居原子之間的相互作用。本節(jié)從數(shù)據(jù)布局、算法優(yōu)化、訪存策略和并行計(jì)算四個(gè)角度對(duì)KC 力場(chǎng)進(jìn)行性能優(yōu)化。

表1 熱點(diǎn)函數(shù)分布

3.1 數(shù)據(jù)布局

法向量用于計(jì)算原子間的橫向距離。KC 力場(chǎng)的原版實(shí)現(xiàn)采用預(yù)處理方案，提前遍歷短鄰接表計(jì)算出每個(gè)中心原子i的法線并保存在多維數(shù)組中。由于受力計(jì)算僅需要中心原子i的法線，不依賴于其他計(jì)算過(guò)程，因此本工作采用“即算即用”策略，存儲(chǔ)結(jié)果僅需常量數(shù)組，避免多次遍歷短鄰接表，從而節(jié)省了帶寬和內(nèi)存資源。同時(shí)，也對(duì)數(shù)組的訪問(wèn)順序進(jìn)行了優(yōu)化。這一過(guò)程如圖3 所示。

圖3 法向量循環(huán)計(jì)算

另外，本文工作還對(duì)計(jì)算受力所需的原子信息進(jìn)行結(jié)構(gòu)體打包操作，包括原子坐標(biāo)、類型和tag 標(biāo)記，共32字節(jié)，這是計(jì)算每對(duì)原子受力時(shí)需要頻繁讀取的數(shù)據(jù)，打包操作極大提升了軟件Cache 讀取數(shù)據(jù)時(shí)的訪存效率。

3.2 算法優(yōu)化

3.2.1 消除冗余計(jì)算

因?yàn)閷娱g作用力只涉及計(jì)算本地中心原子i的短鄰居，而鏡像原子不參與層間相互作用的計(jì)算，所以本工作在構(gòu)建短鄰居列表時(shí)只需搜索本地原子，從而消除搜索鏡像原子的短鄰居帶來(lái)的冗余計(jì)算。

3.2.2 多核心融合

由于計(jì)算法向量需要遍歷短鄰居列表，因此可以將短鄰居列表的構(gòu)建與法向量的計(jì)算融合。計(jì)算長(zhǎng)范圍作用力時(shí)，引力項(xiàng)和斥力項(xiàng)均需遍歷長(zhǎng)鄰接表，兩者之間沒(méi)有依賴關(guān)系，因此可以融合（如圖4 所示）。經(jīng)過(guò)兩次融合后的算法，均需遍歷本地中心原子i，為了復(fù)用中心原子信息，兩者可以放在同一循環(huán)中計(jì)算。

圖4 受力計(jì)算融合

3.2.3 設(shè)置緩存區(qū)

短鄰居受力是嵌套在長(zhǎng)鄰接表的循環(huán)中計(jì)算的，設(shè)長(zhǎng)鄰接表長(zhǎng)度為N，短鄰接表長(zhǎng)度為n，則累加短鄰居受力時(shí)的寫(xiě)內(nèi)存操作為N×n次。由于短鄰居列表只由中心原子決定，因此本文為當(dāng)前中心原子設(shè)置一個(gè)Buffer緩沖區(qū)，暫存短鄰居受力結(jié)果，從而將從核寫(xiě)主存的次數(shù)減少到N×1 次，節(jié)約帶寬資源。

3.3 并行計(jì)算

3.3.1 線程級(jí)并行

LAMMPS 中已實(shí)現(xiàn)基于MPI 的進(jìn)程級(jí)并行，每個(gè)進(jìn)程的任務(wù)是計(jì)算nlocal 個(gè)本地中心原子的受力。在此基礎(chǔ)上，本工作充分利用從核資源實(shí)現(xiàn)線程級(jí)并行，將nlocal 個(gè)中心原子劃分為更細(xì)粒度的從核任務(wù)，分配給64 個(gè)從核并行計(jì)算，以提升模擬效率。

3.3.2 行加鎖軟件Cache 累加受力

由于SW39000 不支持主從核內(nèi)存的寫(xiě)回一致性檢查，多個(gè)從核同時(shí)并行寫(xiě)回受力結(jié)果時(shí)會(huì)產(chǎn)生寫(xiě)沖突問(wèn)題。本文采用基于原子操作FAAL 指令的行加鎖Cache策略，實(shí)現(xiàn)多個(gè)從核并行累加中心原子及其鄰居原子的受力至主存中。

3.3.3 從核通信累加能量

能量的計(jì)算需要對(duì)全局變量累加，為避免多個(gè)從核產(chǎn)生寫(xiě)沖突，本文的實(shí)現(xiàn)是每個(gè)從核先計(jì)算各自的能量，在整個(gè)循環(huán)結(jié)束后，利用從核通信機(jī)制，將64 個(gè)從核的結(jié)果歸約到從核0，最后由從核0 寫(xiě)回至主存。

3.3.4 利用C++特性

新一代神威的C++特性為從核編程提供方便。主核調(diào)用從核函數(shù)時(shí)，參數(shù)的傳遞僅需一個(gè)this 指針。對(duì)于代碼中的常量if 分支判斷，采用基于模板參數(shù)的編程方法，在預(yù)編譯時(shí)生成多個(gè)不同函數(shù)，從而避免每次循環(huán)迭代的判斷開(kāi)銷(xiāo)。

3.4 訪存策略

3.4.1 硬件Cache

新一代神威超級(jí)計(jì)算機(jī)支持硬件Cache 功能，可以方便地讀取離散的數(shù)據(jù)。硬件Cache 的大小有三種不同配置：0 KB，32 KB 和128 KB。本工作充分利用自動(dòng)硬件Cache 這一特性，解決常量數(shù)組和訪問(wèn)頻率較低數(shù)組的離散讀取問(wèn)題。

3.4.2 軟硬件協(xié)同Cache 策略

LAMMPS 采用鄰接表數(shù)據(jù)結(jié)構(gòu)存儲(chǔ)原子間的關(guān)系，這將導(dǎo)致離散訪存問(wèn)題。本文提出一種軟硬件Cache 協(xié)同策略解決KC 力場(chǎng)計(jì)算過(guò)程中的離散訪存問(wèn)題。

對(duì)于頻繁訪問(wèn)的原子信息，本文采用軟件Cache 策略，一次讀取一個(gè)Cache 行，原子信息以結(jié)構(gòu)體的形式存儲(chǔ)，進(jìn)一步提升軟件Cache 通過(guò)DMA 訪存的效率。

4 整體性能評(píng)估

本章利用5016 原子的六層石墨烯算例Case0（如表2 所示）對(duì)性能進(jìn)行整體評(píng)估。通過(guò)LAMMPS 中的replicate 命令將此算例進(jìn)行不同規(guī)模擴(kuò)展（如表3 所示），分別用于加速比、可擴(kuò)展性和不同平臺(tái)的性能對(duì)比測(cè)試。

表2 六層石墨烯算例Case0

表3 不同規(guī)模的算例

4.1 加速比

圖5 展示了各個(gè)優(yōu)化版本的加速比。通過(guò)數(shù)據(jù)組織和算法優(yōu)化，重構(gòu)代碼相對(duì)原版實(shí)現(xiàn)8 倍加速比。利用從核進(jìn)行線程級(jí)并行后，整體性能加速20 倍，此時(shí)硬件Cache 的大小配置為0，即從核通過(guò)GLD 操作訪問(wèn)主存空間。當(dāng)硬件Cache 的大小配置為32 KB 和128 KB 時(shí)，整體性能分別加速109 倍和120 倍。采用軟硬件協(xié)同Cache 策略時(shí)，整體性能實(shí)現(xiàn)155 倍加速，此時(shí)硬件Cache 配置為32 KB，軟件Cache 配置為128 KB（256 行，16 列）。

圖5 不同優(yōu)化版本的加速比

4.2 其他平臺(tái)對(duì)比

本節(jié)在不同的平臺(tái)（如表4 所示）對(duì)比了KC 力場(chǎng)模擬二維材料體系時(shí)的整體性能，如圖6 所示。本工作（SW-OPT）在SW39000 單核組上具有良好的整體性能表現(xiàn)，比單個(gè)E5 2680-v3 核心快3.2～3.7 倍，比單個(gè)Gold 6138 核心快1.6～1.85 倍。

表4 測(cè)試平臺(tái)參數(shù)

圖6 不同平臺(tái)性能對(duì)比

4.3 可擴(kuò)展性

本節(jié)測(cè)試了弱強(qiáng)擴(kuò)展性以及并行效率，計(jì)算方式如下：

其中，TN表示N個(gè)進(jìn)程并行時(shí)的模擬性能，即每秒模擬的時(shí)間步（Timesteps/Second），base=6。

弱擴(kuò)展性的表現(xiàn)如圖7 所示。以6 進(jìn)程為基準(zhǔn)，使用Case0 和S0 兩種算例，每次測(cè)試將進(jìn)程數(shù)和原子數(shù)目同時(shí)等比例擴(kuò)展2 倍。Case0 擴(kuò)展到1 536 進(jìn)程時(shí)，原子總數(shù)為1 284 096，并行效率為92.7%，模擬性能約為2 ns/day；S0 擴(kuò)展到1 536 進(jìn)程時(shí)，原子總數(shù)為10 272 768，并行效率為99%。由此可得，單進(jìn)程負(fù)載的原子數(shù)目越多，計(jì)算密度越大，并行效率越高。

圖7 弱擴(kuò)展性

強(qiáng)擴(kuò)展性的表現(xiàn)如圖8 所示。使用S1 和S2 算例測(cè)試從6 進(jìn)程擴(kuò)展至1 536 進(jìn)程時(shí)的并行效率分別為48.3%和78.1%。從圖8 中可以看出，算例S2 的并行效率高于S1。這是因?yàn)榭傇訑?shù)不變時(shí)，隨著進(jìn)程數(shù)目增加，單進(jìn)程分配的原子數(shù)目減少，計(jì)算密度降低，性能更容易受到其他因素影響，如進(jìn)程間通信。

圖8 強(qiáng)擴(kuò)展性

5 結(jié)論

本文基于新一代神威超級(jí)計(jì)算機(jī)，對(duì)二維材料模擬中的層間力場(chǎng)KC 開(kāi)展移植、算法優(yōu)化和線程級(jí)從核并行等工作，針對(duì)層間力場(chǎng)的熱點(diǎn)函數(shù)進(jìn)行了深入分析，采用一系列的優(yōu)化方案，有效提高了二維材料的模擬效率。研究表明：

（1）通過(guò)重新組織數(shù)據(jù)和算法優(yōu)化，使整體性能相對(duì)原版加速8 倍；

（2）利用從核實(shí)現(xiàn)線程級(jí)并行，在不使用任何訪存策略時(shí)，實(shí)現(xiàn)了20 倍的整體加速；

（3）采用軟硬協(xié)同Cache 策略優(yōu)化訪存性能，實(shí)現(xiàn)了155 倍的整體加速。

與其他平臺(tái)性能相比，本工作具有更為優(yōu)異的性能表現(xiàn)，比E5 2680-v3 單核心快3.2～3.7 倍，比Gold 6138單核心快1.6～1.85 倍。通過(guò)對(duì)整體性能進(jìn)行可擴(kuò)展性測(cè)試發(fā)現(xiàn)，弱擴(kuò)展性的并行效率為92.7%～99%，本工作有助于提高更大規(guī)模的二維材料體系的模擬效率，同時(shí)也為解決神威上的類似問(wèn)題提供參考方案。