晶硅分子动力学模拟的GPU加速算法优化_林琳.pdf

第 49卷 第 4期2023年 4月Computer Engineering 计算机工程晶硅分子动力学模拟的 GPU加速算法优化林琳1，祝爱琦2，赵明璨2，张帅2，叶炎昊2，徐骥2，韩林3，赵荣彩3，侯超峰2（1.郑州大学 信息工程学院，郑州 450001；2.中国科学院过程工程研究所，北京 100190；3.郑州大学 国家超级计算郑州中心，郑州 450001）摘要：分子动力学(MD)模拟是研究硅纳米薄膜热力学性质的主要方法，但存在数据处理量大、计算密集、原子间作用模型复杂等问题，限制了 MD 模拟的深入应用。针对晶硅分子动力学模拟算法中数据访问不连续和大量分支判断造成并行资源浪费、线程等待等问题，结合 Nvidia Tesla V100 GPU 硬件体系结构特点，对晶硅 MD 模拟算法进行设计。通过全局内存的合并访存、循环展开、原子操作等优化方法，利用 GPU 强大并行计算和浮点运算能力，减少显存访问及算法执行过程中的分支冲突和判断指令，提升算法整体计算性能。测试结果表明，优化后的晶硅 MD模拟算法的计算速度相比于优化前提升了 1.691.97 倍，相比于国际上主流的 GPU 加速 MD 模拟软件 HOOMD-blue和 LAMMPS分别提升了 3.203.47倍和 17.4038.04倍，具有较好的模拟加速效果。关键词：分子动力学；图形处理器；合并访存；循环展开；原子操作；性能优化开放科学（资源服务）标志码（OSID）：中文引用格式：林琳，祝爱琦，赵明璨，等.晶硅分子动力学模拟的GPU加速算法优化 J.计算机工程，2023，49（4）：166-173.英文引用格式：LIN L，ZHU A Q，ZHAO M C，et al.GPU-accelerated algorithm optimization for molecular dynamics simulation of crystalline silicon J.Computer Engineering，2023，49（4）：166-173.GPU-Accelerated Algorithm Optimization for Molecular Dynamics Simulation of Crystalline SiliconLIN Lin1，ZHU Aiqi2，ZHAO Mingcan2，ZHANG Shuai2，YE Yanhao2，XU Ji2，HAN Lin3，ZHAO Rongcai3，HOU Chaofeng2（1.School of Information Engineering，Zhengzhou University，Zhengzhou 450001，China；2.Institute of Process Engineering，Chinese Academy of Sciences，Beijing 100190，China；3.National Supercomputing Center in Zhengzhou，Zhengzhou University，Zhengzhou 450001，China）【Abstract】Molecular Dynamics（MD）is one of the main methods used to study the thermodynamic properties of silicon nano-films；however，these studies have problems such as processing massive amounts of data，computational intensity，and complex interatomic interaction，which limit the comprehensive application of MD simulations.To address discontinuities in data access and branch judgments causing the parallel waste of resources and thread waiting in the crystalline silicon MD simulation algorithm，this study combines the characteristics of Nvidia Tesla V100 Graphics Processor Unit（GPU）hardware architecture to design the crystalline silicon MD simulation algorithm.Global memory optimization methods such as coalesced access，loop unrolling，and atomic operation are designed for the MD simulation algorithm.The combination of optimization design and a GPU with powerful parallel and floating-point computing capabilities reduces branch conflicts and judgment instructions during memory access and algorithm execution and improves the overall computing performance of the algorithm.The test results show that the optimized crystal silicon MD simulation algorithm is 1.69-1.97 times faster than the unoptimized algorithm.The optimized algorithm performs 3.20-3.47 and 17.40-38.04 times better than the GPU-accelerated MD simulation software HOOMD-blue and LAMMPS，respectively.The simulations achieve good computation performance.【Key words】Molecular Dynamics（MD）；Graphics Processor Unit（GPU）；coalesced access；loop unrolling；atomic operation；performance optimizationDOI：10.19678/j.issn.1000-3428.0064457基金项目：国家自然科学基金（21776280，22073103）；北京市自然科学基金（JQ21034）；河南省重大科技专项（201400211300）。作者简介：林琳（1997），女，硕士研究生，主研方向为分子动力学模拟；祝爱琦、赵明璨、张帅、叶炎昊，博士研究生；徐骥，副研究员、博士；韩林，副教授、博士；赵荣彩，教授、博士；侯超峰（通信作者），副研究员、博士。收稿日期：2022-04-13 修回日期：2022-06-02 Email：先进计算与数据处理文章编号：1000-3428（2023）04-0166-08 文献标志码：A 中图分类号：TP391.9第 49卷 第 4期林琳，祝爱琦，赵明璨，等：晶硅分子动力学模拟的 GPU加速算法优化0概述分子动力学（Molecular Dynamics，MD）模拟是研究纳微系统热力学平衡态和非平衡态过程的主要方法，应用牛顿定律积分来模拟原子运动，包括计算各个原子的受力、更新原子位置与速度等信息，最终结合统计物理研究体系的热力学物理性质1。MD模拟能够细致刻画原子振动周期内的微观过程，已成为研究诸如硅等半导体介电材料及低维纳米结构导热性质的主要手段。目前，在晶硅半导体材料的MD 模拟中，原子间作用势函数经常采用多体势2，由于原子数目多、计算密集、原子间作用力复杂，因此需要处理巨量的计算负载3-5。近年来，随着异构计算系统的兴起，利用图形处理器（Graphics Processing Unit，GPU）进行大规模高性能计算展现出了巨大的行业发展潜力。通过大量线程的集中并行计算和高速的存储带宽，GPU 相比中央处理器（Central Processing Unit，CPU）呈现出更高的峰值计算能力6-7。MD 模拟中多体原子间作用力的计算具有计算密集的特征，同时作为该模拟算法的核心部分，对于数亿原子规模以上的处理计算量巨大。因此，基于 GPU 的高性能计算将提供一种有效的解决方案。当前，LAMMPS8、HOOMD-blue9、NAMD10、GROMACS11等多个国际上主要的 MD 模拟软件均支 持 GPU 加 速 计 算12。晶 硅 MD 模 拟 采 用 多 体Tersoff13势函数描述 Si-Si原子之间的作用力，多体Tersoff 势函数的力求解是晶硅 MD 模拟算法的核心。目前，关于面向 GPU的 MD模拟的优化工作，国内外已针对其算法、内存读取等方面进行了设计。例如：文献 6 阐述了访存效率、硬件资源特征对性能的影响及优化方式；文献 2-3 介绍内存资源利用率、算法设计的优化方式；文献 14 针对软件缓存和最小化主存储器访问进行优化设计。本文针对晶硅分子动力学模拟中计算负载量大的问题，利用具有强大并行计算能力和高浮点性能的图形处理器，对晶硅 MD 模拟算法进行性能优化。利用合并访存的优化方法，改进 GPU 加速算法的访存效率和计算吞吐量。采用循环展开的优化方法，提高 GPU 并行计算资源利用率。通过原子操作以减少显存访问，提高算法整体计算效率。1算法原理与实现分子动力学模拟采用Tersoff势函数描述晶硅的原子间作用。算法实现包括硅原子邻居列表的建立、原子势能与力的计算、原子位置与速度的更新计算等过程。1.1Tersoff势函数Tersoff 势函数从量子力学理论出发，基于共价键体系中键级与周围原子环境的关系，建立的一种原子间多体相互作用势模型。该类势函数能较好地描述晶硅原子间的作用力，主要包括 T115、T216、T3133个模型。本文主要面向晶体硅的热传递和热物性的变化，采用 T3势函数来模拟晶体硅。T3势函数的数学模型13如式（1）所示：E=iEi=12i=1Nj iVij（1）其中：E为体系总势能；Ei为原子 i的位置势能；Vij为原子 i与原子 j之间的相互作用能。Vij的表达式如下：Vij=fc(rij)fR(rij)+bijfA(rij)（2）其中：rij是原子 i与 j之间的键长；fc是截断函数，随着原 子 间 距 离 的 增 大 逐 渐 减 弱 降 低 为 0，表 达 式 如式（3）所示；R、S是截断半径；fR是排斥项对势，fA是吸引项对势，表达式分别如式（4）和式（5）所示；bij反映了原子间的多体作用，表达式如式（6）所示。fc(r)=|1,rij Rij12+12cos|(rij-Rij)Sij-Rij,Rij rij Sij 0,rij Sij（3）fR(r)=A exp(-1r)（4）fA(r)=-B exp(-2r)（5）bij=(1+nnij)-0.5/n（6）其中：为键级系数；ij描述了原子 i与周围原子间的角势能，表达式如式（7）所示；c、d、h 为常量系数；ijk为 rij与 rik之间的键角，表达式如式（8）所示。ij=k jfC(rik)g(ijk)（7）g()=1+c2d2-c2d2+(h-cos ijk)2（8）部分常量参数的值如表 1所示。1.2晶硅 MD模拟算法异构实现NVIDIA Tesla V100 GPU 处理器采用 NVIDIA Volta 架构，相较于 CPU 更擅长处理数据规模大、计算量大、计算密集的并行计算任务。结合 CPU+GPU的异构计算系统特性，在晶硅 MD 模拟算法的执行中，原子信息的初始化、邻居列表的建立在 CPU上执表 1Tersoff势函数参数设置 Table 1Parameter setting of Tersoff potential function参数名A/eVB/eV1/12/1R/S/106参数值1 830.8471.182.479 91.732 22.73.01.099 91672023年 4月 15日Computer Engineering 计算机工程行，对于算法中邻居列表的更新、Tersoff势函数和原子间作用力的计算、原子位置与速度的更新部分映射在 GPU 上并行执行，最后将 GPU 上并行计算的结果传输到 CPU内存。图 1为晶硅分子动力学模拟算法流程。首先，读取原子的参数信息到CPU内存，在CPU上对其进行内存分配及数据处理；其次，通过 cudaMemcpy（）函数将数据信息传输到 GPU全局内存；之后，通过 GPU并行计算出各原子邻居信息，更新不断运动原子的新邻居列表数据；然后，在 GPU上实现 Tersoff势函数的数学计算，得到原子势能与原子间的作用力，并通过运动方程积分获得原子新位置与新速度；最后，将更新后各原子信息传输到 CPU，并输出查看和统计。随着算法时间步长的逐渐迭代，重复上述 GPU上的并行计算，更新原子邻居列表、位置、速度、作用力等数据，达到热力学平衡态。由于采用邻居列表算法，因此整体模拟的计算量由全模拟域原子相互作用的O（n2）降低至当前的 O（n），其中 n为模拟区域的原子总数。Tersoff势函数作为晶硅分子动力学模拟中的核心计算函数，计算原子 i、j 与原子 i、k 之间的键角为ijk时，涉及 3个原子，在实现 CPU 串行代码时需要通过三重 for循环来进行计算，且每一层循环都需要读取全部的原子数据，计算量大。因此，实现 Tersoff势函数的并行计算来提高算法的执行效率是十分必要的。算法 1 为 Tersoff 势函数在 GPU 上运行的伪代码，虽然基于 GPU 实现了多线程并行计算，去掉了三重 for 循环的最外层循环，但是仍然存在双重 for循环和大量的密集计算，例如、bij、fc、fR、fA等变量的计算过程，使得并行计算资源利用率低、线程等待问题显著及并行模拟效率差，这些都将成为接下来优化的重点。算法 1 Tersoff势函数的计算过程输入 硅原子的位置信息、模型参数信息输出 硅原子总势能、原子间作用力1.int i=blockIdx.x*blockDim.x+threadIdx.x；/i代表/原子的 id号，即一个线程处理一个原子2.i_x=x i；/原子 i的位置信息3.for j to num_size i do /遍历原子 i的所有邻居原子，/计算其原子势能与原子间作用力4.j_x=neighbor（i，j）；/获取原子 i的邻居原子 j5.rij=distance（i_x，j_x）；/计算原子 i与 j间的距离6.for k to num_size i do7.if（k=j）continue；8.k_x=neighbor（i，k）；/获取原子 i的邻居原子 k9.rik=distance（i_x，k_x）；/计算原子 i与 k间的距离10.（i，j）+=fc（rik）*g（ijk）；/计算原子 i 与周围原子的/角势能11.end pute（bij）；/计算 bij13.PotEnergy+=fc（rij）fR（rij）+bij*fA（rij）；/计算原子i的/势能14.nbf_list j =compute（Fij）；/将原子 i 与邻居原子 j/间的相互作用力存储在数组 nbf_list 中15.nbf_list k =compute（Fik）；/将原子 i与邻居原子 k/间的相互作用力存储在数组 nbf_list 中16.end for17.ppotent i =PotEnergy；/原子i的势能值存储在ppotent/数组中2晶硅 MD模拟算法优化针 对 晶 硅 MD 模 拟 算 法 在 GPU 上 的 并 行 计算17-19面临的原子数量多、计算密集与指令执行复杂问题，将 GPU 硬件体系结构与 MD 模拟算法相结合，基于全局内存的合并访存、循环展开、原子操作的优化方法，充分利用 GPU 并行计算能力，提高整体 MD模拟算法的效率。2.1合并访存优化合并访存指的是相邻的线程访问段对齐的地址，避免线程访问不连续地址造成内存带宽的浪费。晶硅 MD 模拟算法采用 CUDA C 编程语言，对于处理一个包含顶点数据 x、y、z、w 的矢量指令，需要用图 1晶硅分子动力学模拟算法流程Fig.1Procedure of crystalline silicon MD simulation algorithm168第 49卷 第 4期林琳，祝爱琦，赵明璨，等：晶硅分子动力学模拟的 GPU加速算法优化4个流处理器（Stream Processor，SP）并行计算完成矢量操作，虽然这样操作比较灵活，但是对于 4维矢量指令，CUDA会将其先转换为相互独立的 4条标量指令，再将其分配到不同的 SP上并行执行20-21。如图2所示，在MD模拟算法中，位置x、速度v、力f采用 float4矢量类型的结构阵列（Array of Structures，AOS）内存布局方式在GPU上进行存储。根据 CUDA的执行模式，对 float4语句进行计算时，先转换成 4条标量指令，再分配到不同的 SP上执行，对于读或写 x、y、z、w任一字段时，因为剩余3个字段在缓存行上都将被隐式加载，所以这种方式在非必需的数值上浪费缓存空间。结合 GPU体系结构的特点，将 float4进行分解重新分配，使用阵列结构（Structure of Arrays，SOA）内存布局，实现合并访存，如图 3所示。针对 GPU 上内存带宽的使用，将 MD 模拟中原子位置、速度等定义方式采用新的数据结构，定义为float*x，*y，*z，*w，该 SOA 模式存储数据有利于合并内存访问和发挥 GPU 的内存带宽，对全局内存实现更高效的利用。SOA 通过对 MD 模拟中数据结构的更改，达到了内存访问对齐的效果，不仅可减少线程访问资源的浪费，而且可提高 GPU 内存资源利用率，进而提高算法执行效率。2.2循环展开优化循环展开是通过展开代码减少循环次数的指令级并行方法22。在使用循环展开时，将循环主体在代码中进行多次编写，降低代码的循环判断次数和分支出现频率，提高线程的并行效率。通常地，循环体的复制数量被称为循环展开因子，展开因子越大，计算速度会有一定提升，但展开因子过大，意味着程序代码越长，消耗的寄存器越多，反而可能造成执行效率下降，因此选取合理的展开因子，可以减少分支判断消耗，提升算法性能22。当系统初始温度为 300 K 时，在模拟中统计多数原子的邻居个数为 4，通过在作用力计算核函数中加入判断语句，将邻居原子个数为 4 的原子间作用力计算过程采用循环展开的方式，消除嵌套 for循环结构。该优化方法减少了分支判断消耗，增加了更多的独立执行指令，提高了流处理器运算单元的活跃度，进一步改善了力计算核函数的并行执行效率。如图 4所示为部分循环展开过程（即算法 1中第 3行第 11 行的 for 循环代码段如图 4 左侧所示），在未使用循环展开前，当线程并行计算时需执行 4 次分支图 2AOS内存布局Fig.2AOS memory layout图 3SOA内存布局Fig.3SOA memory layout1692023年 4月 15日Computer Engineering 计算机工程判断语句，得到原子 i与邻居原子的角势能。在使用循环展开优化后，如图 4右侧所示，删除了 for循环语句，增加了变量 j 分别在 03 下变量 k 由 03 的独立执行指令，直接计算出原子 i 与各个邻居原子的 值，避免了分支判断指令的消耗，提高了线程的执行效率。2.3原子操作优化在CUDA执行过程中，原子操作可确保多个线程对同一块内存地址操作的写保护，其中每次只允许一个线程对内存中变量进行写操作，并且当一个线程正在对变量进行操作时，其他线程若想对该变量进行操作，则只能等待前一线程执行完成，这样的执行方式可确保变量更新正确23。在晶硅 MD模拟算法中，针对多个原子间作用力计算，使用原子操作atomicAdd（）24-25来计算原子受到的合力并存储到 GPU显存中。在 MD 模拟算法的核心函数 Tersoff原子间力的计算过程中，引入原子操作 atomicAdd（），并使用存取速度快的寄存器暂存原子间的作用力，不仅避免了多次向全局内存（nbf_list 数组）存储的操作，而且避免了在更新原子速度与位移时对原子合力的计算。使用原子操作优化的具体实现步骤为：首先定义临时寄存器变量；然后搜索、计算完成所有邻居原子之间的受力后，将计算结果先存储在临时寄存器中；最后使用原子操作 atomicAdd（），以串行的思想计算原子受到的合力，写回 GPU 全局内存中，便于之后原子位置与速度更新时对数据的读取。通过atomicAdd（）原子操作，Tersoff原子间力的计算过程虽然增加了原子所受合力的计算，但是减少了核心函数中反复向 GPU全局内存中存储各个邻居原子对该原子的单独作用力造成的内存访问和运动积分函数实现过程中对原子间部分作用力的读取及合力的计算所造成的时间损耗，使耦合优化算法性能达到最佳26-27。3计算验证3.1计算环境晶硅 MD 模拟算法通过 CPU+GPU 耦合的异构模式进行计算，设计算法提高计算性能。算法运行的实验环境如表 2所示。3.2计算测试3.2.1正确性分析基于第 2节所述的优化方法，在 Tesla V100 GPU上分别对原子数目为 32 768、262 144、1 000 000、2 097 152、4 096 000的算例进行模拟算法的测试，如图 5所示。由图 5可以看出，当原子数目为 32 768时，在单精度和初始温度为 300 K 的条件下，优化前每个原子平均势能值与优化后每原子平均势能值在4.610 19 eV 平衡值左右波动，证明了晶硅 MD 模拟算法优化前后结果的正确性。图 4的循环展开Fig.4Loop unrolling of 表 2实验环境 Table 2Experimental environment处理器CPUGPU硬件配置Intel Xeon Gold 5218 CPU 2.30 GHz(16 cores per socket,2 CPU sockets per node)NVIDIA Tesla V100,PCIe编译器GCC 4.8.5CUDA 9.0170第 49卷 第 4期林琳，祝爱琦，赵明璨，等：晶硅分子动力学模拟的 GPU加速算法优化3.2.2性能测试与分析在装备 Tesla V100 GPU 的计算节点上对优化前后的算法进行性能测试与分析。分别选取原子数目为 32 768、262 144、1 000 000、2 097 152、4 096 000，测试在单精度和初始温度为 300 K 的条件下，优化前后每 1 000个时间步长（每个时间步对应 1.27 帧/s）的耗时并进行对比分析，如表 3所示。由表 3可以看出，晶硅 MD 模拟算法在使用合并访存的优化方法后，不同原子数目下优化前后算法的加速比分别达到 1.21、1.30、1.29、1.29、1.30 倍，验证了采用的数据结构重构方式，实现了算法的合并访存，进而提升了数据的访存效率和计算的吞吐量。在合并访存优化的基础上，实现循环展开的优化算法，解决算法中 Tersoff核函数循环规模大的问题，即将函数中原子邻居个数为 4 的循环执行完全展开。表 4为采用合并访存和循环展开相结合优化方法前后的算法性能对比。由表 4 可以看出，当原子数目为 1 000 000 时，加速比最大，达到了 1.79 倍，相比于仅使用合并访存优化，性能得到了进一步的提升。可见，采用循环展开的优化方式，降低了算法执行过程中的分支冲突和判断指令的消耗，提高了算法的并行计算效率。在合并访存和循环展开优化方法的基础上，加入原子操作方法，优化前后的算法性能对比如表5所示。由表 5可以看出，当原子数目为 32 7684 096 000时，算法优化前后的加速比分别为 1.69、1.79、1.97、1.87、1.90倍，性能得到很大的提升。因此，在计算原子间作用力时使用原子操作指令，不仅防止了多线程并行计算时对同一地址空间的写冲突，而且减少了原子位置和速度更新时访问全局内存的次数，提高了算法的整体性能。综上所述：在合并访存、循环展开、原子操作 3种优化方法相结合后，当原子数目为 1 000 000时，加速比最大，可达到1.97倍；3种优化方法的结合不仅提高了算法的访存效率、计算吞吐量，而且降低了分支指令的消耗和显存的访问频率，使晶硅 MD模拟算法的整体加速效果达到最佳。3.2.3与LAMMPS和HOOMD-blue软件的性能对比LAMMPS8和 HOOMD-blue9作为应用广泛的开源 MD 模拟软件，采用空间分解并行算法，可支持千万原子以上规模的原子分子体系以及多种常见的势能函数，并行性能较好。LAMMPS8利用 MPI 执行消息传递的 MD 程序，可使用 GPU 加速包来进行加速模拟。HOOMD-blue9是专门支持 GPU 的分子动力学模拟的软件包，能够运用多种力场执行模拟，计算的每一步都在 GPU 上执行。本文在 Tesla V100 GPU 计算节点上分别安装 HOOMD-blue V2.6.0 和lammps-stable_5Jun2019，运行相同规模的晶硅，统计测试结果。图 6为使用 GPU加速包的 LAMMPS和 HOOMD-blue模拟软件与优化后的晶硅MD模拟算法的性能对图 5晶硅 MD模拟算法优化前后每个原子的平均势能变化Fig.5Variation of mean potential energy per atom before and after optimization of the crystalline silicon MD simulation algorithm表 3合并访存优化前后性能对比 Table 3Performance comparison before and after coalesced access optimization原子数目/个32 768262 1441 000 0002 097 1524 096 000优化前耗时/s0.0760.4301.5403.1006.290优化后耗时/s0.0630.3301.1902.4004.830加速比/倍1.211.301.291.291.30表5合并访存、循环展开和原子操作优化前后性能对比 Table 5Performance comparison before and after coalesced access，loop unrolling and atomic operation optimizations原子数目/个32 768262 1441 000 0002 097 1524 096 000优化前耗时/s0.0760.4301.5403.1006.290优化后耗时/s0.0450.2400.7801.6603.310加速比/倍1.691.791.971.871.90表 4合并访存和循环展开优化前后性能对比 Table 4Performance comparison before and after coalesced access and loop unrolling optimizations原子数目/个32 768262 1441 000 0002 097 1524 096 000优化前耗时/s0.0760.4301.5403.1006.290优化后耗时/s0.0520.2800.8601.8403.720加速比/倍1.461.541.791.681.691712023年 4月 15日Computer Engineering 计算机工程比，其中，灰色线条表示使用 GPU加速包的 LAMMPS软件与优化后的晶硅 MD模拟算法模拟的加速比，黑色线条表示使用HOOMD-blue软件与优化后的晶硅MD模拟算法模拟的加速比。由图 6可以看出，针对不同原子数目，优化后算法的计算速度与国际主流 MD模拟软件 HOOMD-blue和使用 GPU加速包的 LAMMPS相比，分别提升了 3.203.47倍和 17.4038.04倍，表明了本文算法的高效性，进一步验证了耦合优化能显著提升算法性能。4结束语本文针对晶硅半导体材料分子动力学模拟算法在GPU上的并行计算问题，设计并实现合并访存、循环展开和原子操作的优化方法，改进了模拟算法的数据访存效率、计算吞吐量和并行计算资源利用率，进而使得算法整体性能达到最佳。实验结果表明，优化后算法的计算速度相比于优化前提升了1.691.97倍，相比于使用 GPU加速包的 LAMMPS和 HOOMD-blue分别提升了17.4038.04倍和3.203.47倍，耦合优化后的算法在GPU上的实现取得了显著的加速效果，为其他分子动力学模拟在 GPU加速器上的优化提供了借鉴。但由于本文设计的晶硅分子动力学模拟算法仅基于单个GPU进行优化，因此后续将考虑多GPU的并行分子动力学模拟实现，进一步改进晶硅分子动力学模拟的时间和空间分辨率。参考文献 1 赵成龙，施慧彬，俞忻峰.基于 OpenCL的 Lammps短程力算法优化研究 J.计算机工程与科学，2015，37（9）：1614-1620.ZHAO C L，SHI H B，YU X F.Short-range force algorithm optimization in Lammps based on OpenCL J.Computer Engineering&Science，2015，37（9）：1614-1620.（in Chinese）2 HOU C F，XU J，WANG P，et al.Efficient GPU-accelerated molecular dynamics simulation of solid covalent crystals J.Computer Physics Communications，2013，184（5）：1364-1371.3 HOU C F，XU J，WANG P，et al.Petascale molecular dynamics simulation of crystalline silicon on Tianhe-1A J.International Journal of High Performance Computing Applications，2013，27（3）：307-317.4 侯超峰，高国贤，徐骥.纳米材料制备及物性测量的虚拟过程工程初探 J.计算机与应用化学，2016，33（9）：1003-1007.HOU C F，GAO G X，XU J.Explorations to the virtual process engineering offabrication and property measurement of nanomaterials J.Computers and Applied Chemistry，2016，33（9）：1003-1007.（in Chinese）5 HOU C F，ZHANG C L，GE W，et al.Record atomistic simulation of crystalline silicon：bridging microscale structures and macroscale propertiesJ.Journal of Computational Chemistry，2020，41（7）：731-738.6 MINKIN A S，KNIZHNIK A A，POTAPKIN B V.GPU implementations of some many-body potentials for molecular dynamics simulationsJ.Advances in Engineering Software，2017，111：43-51.7 潘林.基于GPU的分子动力学势函数及结构特征量模拟系统 D.哈尔滨：哈尔滨工业大学，2014.PAN L.GPU-based molecular dynamics potential function and structural characteristic simulation system D.Harbin：Harbin Institute of Technology，2014.（in Chinese）8 THOMPSON A P，AKTULGA H M，BERGER R，et al.LAMMPSa flexible simulation tool for particle-based materials modeling at the atomic，meso，and continuum scales J.Computer Physics Communications，2022，271：108171.9 ANDERSON J A，GLASER J，GLOTZER S C.HOOMD-blue：a Python package for high-performance molecular dynamics and hard particle Monte Carlo simulations J.Computational Materials Science，2020，173：109363.10 PHILLIPS J C，BRAUN R，WANG W，et al.Scalable molecular dynamics with NAMDJ.Journal of Computational Chemistry，2005，26（16）：1781-1802.11 ABRAHAM M J，MURTOLA T，SCHULZ R，et al.GROMACS：high performance molecular simulations through multi-level parallelism from laptops to supercomputers J.SoftwareX，2015，1/2：19-25.12 胡蓉，阳王东，王昊天，等.基于 GPU加速的并行 WMD算法 J.计算机科学，2021，48（12）：24-28.HU R，YANG W D，WANG H T，et al.Parallel WMD algorithm based on GPU acceleration J.Computer Science，2021，48（12）：24-28.（in Chinese）13 TERSOFF J.Empirical interatomic potential for silicon with improved elastic properties J.Physical Review B，1988，38（14）：9902-9905.14 LIANG J，HUA R，ZHANG H，et al Accelerated molecular dynamics simulation of silicon crystals on TaihuLight using OpenACC J.Parallel Computing，2020，99：102667.15 TERSOFF J.New empirical model for the