材料微观结构演化大规模分子动力学软件比较

聂宁明 胡长军 张云泉 贺新福 张博尧 李士刚

1

2

3

4

1

3

123

中国科学院计算机网络信息中心，北京 100190 北京科技大学计算机与通信工程学院，北京 100083 中国科学院计算技术研究所，北京 100190

4中国原子能科学院研究院，北京 102413 ([email protected])

Comparison of Large-scale Molecular Dynamics Software for Materials Computing

Ningming Nie1 Changjun Hu2 Yuquan Zhang3 Xinfu He4 Boyao Zhang1 Shigang Li3

1

(Computer Network Information Center, Chinese Academy of Science, Beijing 100190)

(School of Computer &Communication Engineering, University of Science &Technology Beijing, Beijing, 100083) (Institute of Computing Technology, Chinese Academy of Science, Beijing 100190) (China Institute of Atomic Energy, Beijing 102413)

2

3 4

Abstract Molecular dynamics simulation is an important tool for studying materials microstructure evolution under radiation effects. Four kinds of mainstream large-scale molecular dynamics open source software for materials computing, LAMMPS, Ls1-MarDyn, IMD and CoMD, are introduced in detail in this paper, with data structures, computational methods, parallel decomposition methods and atomic storage, etc. are compared and analyzed. By the design examples, computational efficiency, parallel performance and memory usage have been tested. It provides the basis for subsequent research to realize large temporal and spatial scale molecular dynamics simulation.

Key words Molecular Dynamics; Large-scale; Data Structure, Material Irradiation Damage

摘要 分子动力学模拟是研究材料辐照效应下微观结构演化的重要工具。本文对四款主流材料计算大规模分子

动力学开源软件LAMMPS，Ls1-MarDyn，IMD和CoMD进行了详细的介绍，从数据结构、计算方法、并行分解方式、原子存储等多个方面进行横向分析比较。通过设计算例测试了各软件的计算效率，并行性能和内存使用情况。为后续实现大时空尺度分子动力学模拟研究提供了研究基础。 关键词 分子动力学模拟；大规模；数据结构；材料辐照损伤 中图法分类号 TP391

材料辐照损伤研究涉及从原子尺度到宏观尺度的9个数量级的跨越，以及皮秒量级原子键断裂过程到几十年工程结构失效和破坏的跨时间尺度的非线性过程。这一演化和发展过程在实验很难观测和发现。因此，目前为止，材料工作者对材料辐照损伤过程中内部微观过程和演化机理尚缺乏深入的了解。高性能计算机技术的发展使人们从微观层次深入理解材料成为可能，成为当今国际材料界的重要研究领域。分子动力学模拟作为原子尺度模拟的经典计算方法，在材料微观结构演化规律研究中处于重要地位。

分子动力学模拟(moleculardynamics simulation , MD)是通过利用计算机求解体系内所有粒子的运动方程来模拟粒子的运动轨迹，从而获得系统的温度、体积、压力、应力等宏观量和微观过程量。从上世纪五十年代发展以来[1,2]，已经在包括物理，计算化学，计算生物，材料科学，药物设计等多个领域有了十分

丰富的应用。分子动力学模拟作为一种非常有效的材料计算技术，已成为与实验同等重要的科学研究的方法。

由于计算能力的限制，目前分子动力学通常模拟的是几万到几千万个原子条件下材料的内部演化过程，而通常材料的微观结构及缺陷尺寸远远大于现有计算能力所实现的计算规模，计算模拟结果与实际材料的微观演化过程及宏观性能等还有很大的差异。因此分子动力学模拟的发展的关键在于空间规模的扩大和时间尺度的延长。只有在原子个数无限多，计算的时间足够长时，才能真实地反映材料性质的宏观行为。

近些年，随着高性能计算机技术的飞速发展，分子动力学计算技术被注入了新的活力。2015年7月，

我国研制的“天河二号”超级计算机再度荣登全球超级计算机500强排行榜榜首，连续五次称冠世界[3]。高性能计算机的发展为大规模分子动力学模拟提供了发展的肥沃土壤。制约实际应用的瓶颈在于大规模分子动力学软件的研发。开展大规模分子动力学模拟软件的研究成为国际上的研究热点。科学家们基于混合架构高性能计算机开发了新的分子动力学程序，例如LLNL与IBM的团队开发了面向高性能计算机的经典分子动力学代码--ddcMD代码。2005年，他们首次从原子层次原子尺度研究了金属钽和铀在超高温高压下的凝固过程，在IBM BlueGene/L的131072个CPU上，达到101.7万亿次/秒的持续计算性能，获得了Gordon Bell奖，该研究对评估核武器储备的可靠性具有重要意义[4]。2007年，他们又实现了第一个基于分子动力学的微米尺度原子级K-H不稳定性数值模拟，在IBM BlueGene/L 全系统212,992个CPU上，获得了54.3万亿次/秒的持续计算性能，再次获得了Gordon Bell奖[5]。近年来美国启动了名为“Spatio-Temporal Frontiers of Atomistic Simulations in the Petaflop Computational World”[6]的科研项目，旨在基于高性能计算机开发可扩展的分子动力学程序，研究晶界与位错的相互作用以及材料中空洞的形核与长大过程，在此基础上采用可扩展并行短程MD程序SPaSM在Roadrunner上对106-1012原子进行了大规模的仿真模拟，达到了369 万亿次/秒持续计算性能[7]。除此之外，针对不同领域应用的通用分子动力学软件也有许多，如LAMMPS[8],NAMD[9],GROMACS[10]等等。它们各有各的特点，在不同领域获得了广泛的应用。本文关注的是在材料计算领域广泛使用，可进行大规模并行计算的开源分子动力学软件。我们选取了LAMMPS，Ls1-MarDyn，IMD和CoMD四款开源分子动力学软件进行详细的比较说明，为我们的后续研发工作提供研究基础与支持。

本文的安排如下：首先是分子动力学软件的介绍，将对选取的LAMMPS，Ls1-MarDyn，IMD和CoMD四款大规模并行计算的分子动力学软件进行详细介绍。接下来，对这四款软件在并行性能，数据结构，原子结构存储等几个方面进行横向比较，并通过数值算例进行验证和分析说明。最后是全文总结，通过总结前文对四款MD软件在大规模并行计算上特点的分析，提出针对辐照效应下金属材料计算特点的数据结构设计思想，为下一步研究工作指出方向。

1. 大规模材料计算分子动力学软件介绍

1.1 LAMMPS（Large-scale Atomic/Molecular Massively Parallel Simulator）软件

LAMMPS软件[8]是由美国Sandia国家实验室开

发的一款通用的大规模分子动力学并行计算软件，是世界上使用者最广泛的几款分子动力学软件之一。从2004年到2011年之间，就达到了8万9千多次的下载量。

LAMMPS主要用于分子动力学相关的一些计算和模拟工作，可在包括气态，液态或者固态相形态、各种系综下，采用不同的力场和边界条件来模拟全原子，聚合物，生物，金属，粒状和粗粒化体系。可实现从几个粒子到上亿粒子体系的分子动力学计算。LAMMPS可提供并支持多种势函数计算，包括对势，如Lennard-Jones(L-J)势（多用于气体，液体分子间作用力计算）、Morse势等，多体势，如Embedded Atom Model(EAM)势（多用于单一金属或合金体系计算）、 modified EAM (MEAM)势等，应用十分广泛。

LAMMPS软件为开源软件，以C++编写，支持用户的自行修改和扩展。LAMMPS支持串行与并行计算。并行LAMMPS采用MPI并行，针对并行计算特点设计数据结构与并行策略，具有良好的并行扩展性。其软件结构可由如图1所示的类结构关系图来展示。

在计算粒子间相互作用时，LAMMPS通过cutoff来设定每个粒子需计算的相互作用粒子范围大小，采用neighbor list，也即近邻列表的数据结构来记录每个粒子的邻居粒子信息，如图2所示。粒子信息存储在数组中，每个粒子的邻近表数组中存储粒子的邻居粒子索引。 在计算时，通过neighbor list来追踪所有与该粒子有相互作用的粒子间作用力等信息。 LAMMPS采用空间分解方法来并行分配模拟区域。将整个模拟区域分解成若干子区域分配到各个处理器上。在每个处理器的计算子区域上设置ghost区域来存储子区域边界原子信息，以便并行计算时各个处理器间的相互通信与计算。LAMMPS提供了两种

图 1 LAMMPS软件结构示意图[8]

Fig. 1 Software Structure Chart of LAMMPS[8]

图 2 LAMMPS neighbor list 数据结构展示 Fig. 2 Data Structure of Neighbor List for LAMMPS

方法来实现负载均衡，Shift方法和RCB方法，分别如图3a和图3b所示。Shift方法通过改变处理器间的分割线，来调整分给每个处理器的空间。最后划分的子区域还是网格的形式。RCB方法是一种类似“砌瓷砖”的方法，它不是把模拟区域划成网格状，而是大小不等的四边形区域，来让每个区域有相等数目的原子。

a) shift方法实现负载平衡 b) RCB方法实现负载平衡

图 3 LAMMPS的负载均衡[8] Fig. 3 Load Balancing Method of LAMMP[8]

此外，LAMMPS已有支持其他架构计算的版本，

包括GPU（支持CUDA和OpenCL），和Intel Xeon Phi等架构，以及OPENMP版本的代码。

1.2. IMD(ITAP Molecular Dynamics)软件 IMD是一款由德国University of Stuttgart于1997年开始研发的用于经典分子动力学模拟的开源软件包[11]，支持包括对势，适用于金属材料的EAM势，适用于共价系统的Stillinger-Weber 和Tersoff 势，以及适用于液体的 Gay-Berne势等多种势函数的分子动力学计算。1999年，IMD在德国Julich超算中心的超级计算机T3E-1200-262上用512核实现了5x109个粒子的分子模拟，取得了当时的世界纪录[12]。

IMD采用的模块化设计，C语言编写，使用标准的MPI库来实现并行消息传递。其软件结构如图4所示。

IMD采用link cell的数据结构来存储数据。先将整个box分成若干cell，每个cell是个结构体，每一

项存储指向原子相关信息数组的指针，具体如图5所示。此外，IMD还提供可选项，对于3D的EAM势计算算例可提供与LAMMPS类似的neighbor list的数据结构进行数据存储[13]。

并行分解时，IMD将处理器排列成三维的笛卡尔网格，然后按区域将box里cells数组细分成相同大小的块，然后将它们分配给不同的处理器。各处理器设置ghost区域用作接收与发送粒子数据的缓冲区。IMD 不提供负载平衡的策略。

IMD除支持MPI的版本之外，还有OPENMP、MPI+OPENMP和GPU版本的实现。

Typesallocate memorydomain decompositionGeometry routines for MPIdirect schemeSet_bufferBuffered schemepotentialIMDforcesdirect cell communicationcommunication schemeFix cellsBuffered schemeintegratetimeioSocket_ioOnline Visualizationpictures 图 4 IMD 软件结构示意图 Fig. 4 Software Structure Chart of IMD

1.3. Ls1-MarDyn(Large systems 1: molecular dynamics)软件

以IMD的开发团队为基础，由德国University of

datafarcort012...n-2n-1eam_rho012...n-2n-1eam_dF012...n-2n-1kraft012...n-2n-1impuls012...n-2n-1sorte012...n-2n-1vsorte012...n-2n-1mass012...n-2n-1Pot_eng012...n-2n-1nummer012...n-2n-1nn_max 图 5 IMD 中cell的数据结构展示 Fig. 5 Data sStructure of Cell in IMD

Stuttgart，University of Kaiserslautern， Technische Universit?t München和University of Paderborn几家大学联合开发的旨在挑战万亿级别的分子模拟的MD代码—“Ls1-MarDyn”[14]，保持了目前最大规模的分子动力学模拟记录：在德国莱布尼兹超级计算中心的超级计算机SuperMUC上用14万核实现了4x1012个粒子的分子模拟[15]。Ls1-MarDyn只支持刚性粒子

和恒体积系宗下的MD计算，适用于短程力计算，不适用于离子计算。

Ls1-MarDyn是以C++编写的开源代码，基于MPI实现并行。Ls1-MarDyn采用模块化设计，将代码分为物理模型计算模块和并行算法模块两类，其软件结构如图6所示。目前，Ls1-MarDyn只支持Lennard-Jones(L-J)势，Tersoff势等对势计算，无法实现EAM等多体势计算。

simulationLinked cellcellAdaptive Linked celldomaindecompositiondecompositionmardyndomainIntegratormolecularinputoutputTimerleapfrogKDdecompositionFig.6 Software Structure Chart of Ls1-MarDyn

针对其适用于气体流体之类分子运动变化频繁的计算的特点，Ls1-MarDyn代码设计linked cell数据结构来进行数据存储，将模拟区域以网格划分成小的单元格cell，单元格长为截断半径，依据分子位置信息，将每个分子对应到相应的cell中，如图7a所示。针对不均匀的分子分布，Ls1-MarDyn对linked cell 数据结构做出改进，对包含粒子数目较多的单元格进行2次分割，设计了自适应linked cell的数据结构，如图7b所示。

Ls1-MarDyn采用了两种并行区域分解的方式。一种是按进程数等分模拟区域，即将n个进程以尽量分配给每个维度相同进程数的原则建立笛卡尔坐标，将模拟区域等分为n个等体积的小区域分到每个进程。这种区域分解方法简单易操作，只需初始化时建立一次即可。但是当计算到后期粒子变化移动较多时，会

图 6 Ls1-MarDyn软件结构示意图

a) 2-d linked cell示意图 b) 2-d adaptive linked cell示意图

图 7 Linked cell数据结构示意图

Fig.7 Data Structure of Linked Cell for Ls1-MarDyn

造成计算负载的不均衡。另一种是基于KD树的区域分解方法，也即用垂直于坐标轴的线将模拟区域递归地划分为两个负载相近的子区域，直到分割的子区域被分配到处理器上（负载近似于1/p），其中负载的计算以cell为基本单位进行估算。Ls1-MarDyn 暂时还未有支持GPU, openmp，intel MIC等异构计算版本代码。

1.4. COMD: Classical Molecular Dynamics Proxy Application

COMD[17]是由美国能源部ASCR资助的Exascale Co-Design Center for Materials in Extreme Environment (ExMatEx)中心设计开发的MD代码，由LANL的 Jamal Mohd-Yusof, Chris Sewell 和Sriram Swaminarayan 基于SPaSM的分子动力学代码来进行开发维护。ExMatEx中心旨在极端环境下材料特性的E级超大规模计算模拟研究，该项目强调多尺度模型，计算机科学，以及通过协同设计实现百亿亿次的超大

规模材料模拟。因此COMD的研发是作为协同设计的工具，方便整个ExMatEx根据需要对新的体系结构、编程模型等进行多尺度建模计算、功能扩展和重新实现，以及性能测试等而开发实现的。COMD是开源代码，目前发布的是1.1版本[18]，只能计算LJ势和EAM势，EAM势的计算只支持fcc晶格结构，只能建立规则的长方体box。COMD是以C语言编写的基础MD程序, 功能不多，主要结构见图8所示。 COMD基于MPI实现并行，采用link cell的数据结构存储数据，空间分解方法进行任务划分。COMD根据cutoff大小划分cell，按坐标XYZ方向进行处理器核数的划分，再将每个方向上的cell均分到每个处理器上，将整个计算区域进行空间的分解。COMD暂时还未考虑负载平衡策略的实现。COMD除MPI版本之外，还实现了OPENMP并行。对于不同计算

架构，COMD还实现了针对各计算架构的版本，如AMD下的OpenCL版，基于NVDIA的GPU下的CUDA 版，正在调试Intel MIC (aka Xeon Phi)和 Intel's FF candidate architecture emulator架构下的众核计算版本[19]。

constantsLj forcememUtilparallelCoMDeamlinkcellstimestepdecompositionperformanceinitAtomsoutputHaloExchange2. 大规模材料计算分子动力学软件分析与比较

2.1 大规模并行分子动力学软件比较

通过上一节的介绍，对LAMMPS，IMD，Ls1-MarDyn以及COMD这四种MD软件的大致情况有了初步的了解，本节将对这四种软件从数据结构、计算方法、与并行分解方式等角度进行横向的详细比较。详见表1。

四种软件都采用空间分解的方式进行并行分解任务。只有LAMMPS和Ls1-MarDyn两种提供了负载平衡技术。相比较而言，基于自适应link cell结构的Ls1-MarDyn在负载均衡的调节上会比LAMMPS更加灵活。

由于我们考虑的都是设置有cutoff的短程力计算，对于此类型计算，通常有两类数据结构用于粒子信息

图 8 CoMD软件结构示意图 Fig. 8 Software structure chart of CoMD

表 1四款软件比较总结

Table 1 Summary of Comparison for Four kinds Software

软件名称 软件类型 编程语言 支持势函数类型

LAMMPS 通用型 C++

对势（LJ等）,多体势（EAM等）等多种

势函数计算方式 数据结构存储类型 数据划分方法

是否实现负载平衡 不同计算架构版本

提供Shift方法和RCB方法两种方法 OPENMP版，GPU（

CUDA

和

OPENMP版，无 OPENMP版，GPU（CUDA）版，OpenCL版，Intel MIC版等

不提供 基于KD树的区域分解方法 暂时未提供 空间分解 Neighbor list

Link cell

(Neighbor list 可选) 空间分解

空间分解

空间分解

Link cell

Link cell

查表插值

查表插值

参数化计算

查表插值

IMD 通用型 C

对势（LJ等）,多体势（EAM等）等多种

Ls1-MarDyn

专用型（气体流体计算） C++

仅支持对势（LJ势）

COMD

通用型（软件模块） C

对势（LJ）,多体势（EAM）

MPI+OPENMP版 GPU版

OpenCL）版，Intel MIC版

的存储，分别为neighbor list和link cell。neighbor list

结构是每个粒子都根据cutoff来确定并存储自己的邻居粒子，每次迭代时只遍历自己的邻居粒子。这种存储方式的好处是计算时很容易找到需计算的对象，并且不会多出额外粒子的计算，计算效率很高。但是存储代价高，每个粒子都要维护自己的list，在运行一定步数之后还需更新自己的list。因此这类数据结构更适合于邻居粒子变化缓慢的应用，如固体材料等的计算。LAMMPS采用的便是neighbor list的数据结构存储原子，IMD也提供此类数据结构的存储选项。Link cell是将计算区域划分成若干cell，按cell来存储粒子信息。计算时，按照cutoff来确定每个cell的邻居cell信息，按cell遍历粒子进行计算。这种方式计算效率低，每次迭代计算时都需遍历本身以及周围邻近cell，会多出非邻居粒子的计算，计算量大。但是好处是存储量较少，且易于做负载平衡。因此此类结构更适于类似气体，流体材料等邻居粒子频繁变化应用的计算。IMD, Ls1-MarDyn和COMD都采用了link cell的数据结构存储粒子信息。

MD计算中对于粒子间作用力计算的关键在于势

相关推荐：