一种在GPU上利用多体作用模型进行粒子计算的方法

摘要

本发明提供一种在GPU上利用多体作用模型进行粒子计算的方法，包括：在加载有GPU的计算机上，将粒子系统中粒子的属性信息以及每个粒子的邻近粒子的标号信息保存到系统内存中；然后再传输到GPU的全局内存中，并在全局内存中为每个粒子对其邻近粒子施加的作用力开辟存储数组；为每一个粒子单独分配一个GPU中的计算线程；根据建立的粒子间多体相互作用模型记下本线程计算粒子的所有受力，同时保存其对邻近粒子的受力；把保存的每一个粒子邻近粒子的受力从全局内存中读入，并加和到相应标号的粒子身上，得到粒子所受到的完整作用力；在GPU上计算每个粒子的势能，并将结果输出到系统内存上，由CPU统计并计算单粒子平均势能。

说明书

技术领域

本发明涉及粒子计算，特别涉及在GPU上利用多体作用模型进行粒子计算的方法。

背景技术

多体相互作用的计算广泛存在于材料、化学、生物、天体物理等多个领域。在材料科学中，多体相互作用的计算是材料原子模拟的一种典型应用。所述的材料原子模拟是材料计算的一种研究方法，该模拟方法通过研究系统中每一个原子的运动变化来研究物质的一些性质。在现有技术中存在多种对物质进行原子模拟的方法，例如，量子力学中的从头算法(abinitio)，嵌入原子方法EAM(Embeded Atom Method)，分子动力学(Molucular Dynamics，MD)模拟方法(如LJ(Lennard-Jones)势函数、Morse势函数等)等。但上述方法在材料的原子模拟中都各有局限之处。

量子力学中的从头算法能够获得势能的精确值，但由于它依赖原子的局域密度函数，因此计算量大、计算时间长，仅适合模拟规模较小的系统，在工程模拟中的实用性不高。嵌入原子方法EAM是基于有效介质方法而来，适合于没有成键取向结构的密堆金属，如过渡金属元素，并且也因为依赖原子局域密度函数，使得计算规模受到限制。在一些工程模拟中，由于采用经验势能代替原子间的实际相互作用势能可能已经足够，因此，可以采用经验势能的计算模型如LJ势函数、Morse势函数等。LJ势是分子动力学模拟中最简单和最常用的一种经验势函数，原子间作用不考虑方向属性，是一种各向同性的二体相互作用。因此常适合于模拟稀有气体、简单金属等一些简单系统。Morse势也是二体相互作用，亦不能准确描述具有很强方向性的共价晶体，只适合于描述一些简单金属。对于像单晶硅那样具有很强方向性的共价晶体上述现有技术实际上并不适合。

在具有很强方向性的共价晶体中，要正确描述原子间的相互作用，不仅要考虑两个原子间的作用，而且要体现成键取向的变化所产生的影响，因此一般要考虑三体及三体以上间的相互作用，多体作用模型就是基于此提出的。在多体作用模型中，原子间作用势函数的通用形式如下：

在上述通用势函数的基础上，对于晶体硅等半导体材料，本领域的技术人员提出了具体的经验势函数，如Stillinger-Weber势函数或Tersoff势函数等。但无论多体作用模型采用了何种势函数，为了加快计算过程、保证计算结果的准确性，都需要在计算机上实现对势函数的计算，从而实现对原子间相互作用的准确求解。

在现有技术中，在计算机上对离散粒子系统中多体作用模型的计算一般是依赖于CPU完成的，由于CPU采用单线程进行计算，因此线程在计算粒子与其邻近粒子的受力时，可以同时在系统内存中对该粒子的邻近粒子受力信息进行读写，从而在一个步骤内得到粒子的完整作用力。但正因为CPU采用单线程进行计算，因而对于离散粒子系统中的多个粒子，只能一个一个依次进行计算，并且由于CPU计算主频的有限能力，使粒子系统的计算规模和计算效率受到了极大的限制。鉴于图形处理单元GPU(Graphics Processing Unit)采用了多通道计算技术，在浮点密集型运算和大量数据的并行处理上有极大的优势，因此，可以考虑在GPU上实现粒子系统中对多体作用模型的计算。现有技术中，在GPU上已经实现了遵从简单运动方程的粒子系统的渲染以及物理模拟中二体相互作用的计算等多种应用。关于上述现有技术的具体内容可见以下参考文献：专利公告号为CN100354894C的专利申请；以及J.A.Anderson，et al.，Journal ofComputational Physics，227，5342-5359，2008。

在所述已知的对二体相互作用在GPU上实现的计算方法中，由于粒子之间的作用是对势，两个粒子之间的作用与第三个粒子无关，所以可以在GPU单个计算线程上对每一个粒子与其周边所有邻近粒子的作用全部计算以求得该粒子所受到的全部作用力。但在多体相互作用中，由于两个粒子之间的作用与第三个粒子有关，所以在计算两粒子之间作用时必须考虑第三个粒子的影响，也即由于牛顿第三定律力的作用是相互的，在计算二粒子受力时，第三个粒子的受力信息也必须被改变。但由于图形处理单元GPU上计算线程同时对系统中的所有粒子进行并行处理，所以在某一线程计算粒子与其邻近粒子的受力时，存在不能同时对该粒子其邻近粒子受力信息进行读写的缺点。

发明内容

本发明的目的是克服现有技术的上述不足，针对采用多体相互作用模型的离散粒子系统，提供一种基于图形处理单元GPU的加速计算实现方法。

为了实现上述目的，本发明提供了一种在GPU上利用多体作用模型进行粒子计算的方法，该方法在加载有图形处理单元GPU的计算机上实现，所述GPU包括计算内核、全局内存、共享内存、常量内存以及纹理内存；该方法包括：

步骤1)、在加载有GPU的计算机上，将粒子系统中粒子的属性信息以及每个粒子的邻近粒子的标号信息保存到所述计算机的系统内存中，所述粒子的属性信息包括初始的位置、速度、加速度、势能，以及粒子标号、粒子种类；

步骤2)、将所述系统内存中的粒子属性信息以及每个粒子的邻近粒子标号信息传输到所述GPU的全局内存中，并在所述全局内存中为每个粒子对其邻近粒子施加的作用力开辟存储数组；

步骤3)、为所述粒子系统中的每一个粒子单独分配一个所述GPU中的计算线程；

步骤4)、在为粒子所分配的线程上，根据所述的多体作用模型计算线程所对应粒子的势能，然后通过对势能求导的力计算模型获取该粒子所受到的作用力，并将粒子对其各个邻近粒子的作用力存储在所述GPU的共享内存中，待与所有邻近粒子相互作用计算完毕后，再将所述粒子对邻近粒子的作用力写回到步骤2)所得到的在全局内存中开辟的力存储数组中；

步骤5)、在为粒子所分配的线程上，根据本线程计算的对应粒子标号，把步骤4)中该粒子邻近粒子所在线程计算的对本粒子的受力从所述全局内存的存储数组中读出，然后再累加到步骤4)计算得到的粒子所受到的作用力上，得到粒子所受到的完整作用力；

步骤6)、根据步骤5)计算得到的粒子所受到的完整作用力确定粒子的加速度，然后由运动方程结合所述属性信息中初始的位置、速度信息，计算下一时刻的速度与位置；

步骤7)、重复执行上述的粒子作用力计算、加速度计算、速度与位置计算操作，直至整个粒子系统达到稳定状态或满足用户所要求的重复计算次数。

上述技术方案中，在步骤7)的重复计算过程中，在执行一定次数的重复计算后，对所述系统各个粒子在各次计算过程中所得到的势能进行累加，然后根据累加的结果计算所述系统粒子的平均势能。

上述技术方案中，在所述的步骤3)中，所述的为每一个粒子单独分配一个计算线程包括：

为所述粒子系统做网格单元的划分，然后为每一个网格单元分配一个线程块，使得所述线程块中的线程数多于所述网格单元中的粒子数。

上述技术方案中，在所述的步骤3)中，所述的为每一个粒子单独分配一个计算线程包括：

根据GPU中线程块所开辟的线程数目，对所述的粒子系统做网格单元的划分，使得所述网格单元中的粒子数少于所述线程块中的线程数。

上述技术方案中，所述的多体作用模型包括Tersoff多体作用模型、Tersoff-Brenner多体作用模型、Stillinger-Weber多体作用模型。

上述技术方案中，在所述的步骤4)中，所述的力计算模型包括：

$F=-▿E$

其中的E表示粒子的势能，表示对粒子位置的求导计算。

上述技术方案中，在所述的步骤5)中，所述的粒子邻近粒子所在线程在步骤4)计算的对本粒子的受力从所述全局内存的存储数组中读出后，存入到共享内存，在所述共享内存中实现与步骤4)计算得到的粒子所受到的作用力的累加，然后将粒子所受到的完整作用力再写回到全局内存中。

上述技术方案中，在所述的步骤5)中，所述的粒子邻近粒子所在线程在步骤4)计算的对本粒子的受力从所述全局内存的存储数组中读出后，直接在所述的全局内存中实现与步骤4)计算得到的粒子所受到的作用力的累加，然后将所得到的结果在全局内存中存储。

上述技术方案中，在所述的步骤6)中，所述的运动方程包括：

$m\frac{\mathrm{dv}}{\mathrm{dt}}=\mathrm{\Sigma F}+F_{\mathrm{external}}$

其中，m表示粒子质量；v表示粒子速度；F表示粒子所受邻近粒子的作用力，F_external表示粒子所受的外部作用力。

本发明的优点在于：

1、本发明利用GPU可对大量数据进行并行处理的特点，对粒子系统中的各个粒子同时进行计算，加快了整体的运算效率，扩大了计算规模。

2、本发明在计算时把粒子的属性信息保存于全局内存，并在共享内存中对其进行计算或读写，提高了算法的通用性和可扩展性，程序开发方便，计算效率明显提高。

3、本发明可根据线程块中拟开辟的线程数目合理设置每个线程块对应网格的粒子数目，或根据系统每个网格中的粒子数目确定每个线程块中需开辟的线程数目，这可以最大程度发挥GPU计算设备的性能。

4、本发明克服了GPU在并行计算采用多体相互作用模型的离散粒子系统时不能同时对每个粒子的邻近粒子受力信息进行读写的缺点，维护了计算的完整性和准确性。

5、可推广到其它采用多体相互作用模型的离散粒子系统基于图形处理单元GPU的加速求解。

附图说明

图1为本发明方法的流程图；

图2为采用本发明的方法与采用现有方法进行原子模拟计算的计算性能比较图；

图3为GPU上线程块中线程开辟数目一定时，原子模拟计算时间与网格尺寸间的变化关系图；

图4为当网格尺寸恒定时，线程块中开辟的线程数目与原子模拟计算效率间的变化关系图；

图5为本发明在一个实施例中对硅晶材料的计算模拟结果图；

图6为本发明在一个实施例中对硅晶体系单原子平均势能随时间变化的计算模拟结果图。

具体实施方式

下面结合附图和具体实施方式对本发明加以说明。

在一个实施例中，在加载有图形处理单元GPU的计算机上通过Tersoff多体作用模型实现了对单晶硅体系的原子模拟。在该实施例中，计算机系统所采用的GPU具有如下配置：集成128颗计算内核，包括1.6G字节的全局内存(Global Memory)，64K字节的常量内存(Constant Memory)，以及共享内存(Shared Memory)和纹理内存(Texture Memory)在内的多种内存，通过CUDA(Compute Unified Device Architecture)编译环境控制图形处理单元GPU上的计算行为。

本实施例中所涉及的Tersoff多体作用模型是现有技术，为了方便理解以及下文中说明的需要，首先对Tersoff多体作用模型加以说明。Tersoff多体作用模型属于多体作用模型中的键序(Bond Order)模型，该模型的特点是原子之间的相互作用依赖于原子所处的局域环境，周边配位原子数目愈多，原子之间键合作用愈弱。下面所述的Tersoff多体作用模型公式用于计算一个原子系统中的总势能，其表达如下：

$E_{\mathrm{total}}=\frac{1}{2}\underset{i\ne j}{\Sigma }{V}_{\mathrm{ij}}\left(r_{\mathrm{ij}}\right)=\frac{1}{2}\underset{i\ne j}{\Sigma }{f}_C\left(r_{\mathrm{ij}}\right)[V_R\left(r_{\mathrm{ij}}\right)+B_{\mathrm{ij}}{V}_A\left(r_{\mathrm{ij}}\right)]$

其中，V_R(r_ij)＝Aexp(-λ₁r_ij)

V_A(r_ij)＝-Bexp(-λ₂r_ij)

$f_C\left(r\right)=\left(\begin{array}{cc}1,& r<R-D\\ \frac{1}{2}-\frac{1}{2}\sin [\frac{\pi }{2}(r-R)/D],& R-D<r<R+D\\ 0,& r>R+D\end{array}\right)$

B_ij＝(1+βⁿζ_ijⁿ)^-1/2n

${\zeta }_{\mathrm{ij}}=\underset{k\ne i,j}{\Sigma }{f}_C\left(r_{\mathrm{ik}}\right)g\left({\theta }_{\mathrm{ijk}}\right)\exp \left[{{\lambda }_3}^3{(r_{\mathrm{ij}}-r_{\mathrm{ik}})}^3\right]$

$g\left(\theta \right)=1+\frac{c^2}{d^2}-\frac{c^2}{d^2+{(h-\cos \theta )}^2}$

上述公式中，V_R是排斥项，V_A是吸引项，B_ij是一个与键角θ_ijk有关的系数，包含了原子间多体相互作用，反映了局域原子环境，下标i、j、k表示原子的标号，r_ij为原子之间的距离，θ_ijk为键ij和键ik之间的夹角。

上述公式中各个因子的具体取值会根据材料的不同而发生变化，以本实施例中所采用的单晶硅为例，各个因子的具体取值如下：

A＝3264.7eV，B＝95.373eV，

β＝0.33675，n＝22.956，R＝3.0，D＝0.2，

c＝4.8381，d＝2.0417，λ₃＝0.0，h＝0.0

关于Tersoff多体作用模型的具体内容在参考文献1“J.Tersoff，Physical Review B，37(12)，6991-7000，1988”中有详细说明。

在本实施例中，假设所要模拟的单晶硅原子体系为立方空间，其三个方向空间特征长度3.2nm，包含1728个原子，初始体系温度为292K，每个模拟时间步大约为0.17fs。周期性边界条件控制，时间步长积分格式采用蛙跳算法。

下面结合图1，对如何在GPU上通过Tersoff多体作用模型实现对单晶硅体系的原子模拟加以说明。

对具有1728个原子的单晶硅体系进行原子模拟，首先需要将系统原子的属性信息，如位置、速度、加速度、势能、粒子标号、粒子种类等，通过计算机系统上的CPU存储到计算机的系统内存中，同时系统内存中存储每个原子的邻近原子标号信息。这些存储数据在后续的操作中会被用到。需要说明的是，在本实施例中，由于模拟系统设定的材料温度没有达到熔点，所以系统中硅原子能够保持完整的晶胞结构，每个原子具有四个最邻近原子，它们的标号信息能够预先被指出。

在计算机的系统内存中完成对原子的属性信息、邻近原子的标号信息的存储后，在GPU的全局内存上同样为计算所需的所有粒子属性信息、及每个粒子的邻近粒子标号信息开辟空间，然后把系统内存中所存储的原子的属性信息、邻近原子的标号信息传输到GPU全局内存对应的数组中，并在全局内存中为每个原子的所有邻近原子开辟与本原子发生作用的作用力储存数组，该数组按照原子顺序标号依次存放每个原子的四个邻近原子受到的作用力，并在初始时对此数组赋值为0。原子的属性信息存入GPU的全局内存后，在后续的势函数计算以及原子间作用力的计算中会被用到。

在利用GPU进行原子模拟的过程中，需要保证一个原子对应GPU上一个单独的线程以实现模拟计算，因此在前一步骤已经给出所要模拟的原子系统的包括原子个数在内的基本信息的基础上，在此实现计算系统中原子与GPU中线程间的对应。从前面的说明中已经知道，本实施例所要模拟的原子系统包括1728个原子，对此原子系统做网格单元的划分，即将用于表示原子系统的立方空间按照空间体积划分为多个网格单元，每个网格单元内有一定数量范围的原子，并且让每个网格单元对应GPU上开辟的一个线程块(Thread block)。对网格单元的划分应当满足一定的标准，如果网格单元中的粒子数目过低，则不能充分发挥GPU计算设备的性能，如果网格单元中的粒子数目过高，则每个线程块中的线程不足以计算所有的粒子，会影响计算的准确性。本实施例中，网格单元在三个空间方向上的特征尺寸取1.04nm，即将包含1728个原子所在的原子系统在空间上分为27个网格单元，则每个网格单元内平均含有64个原子。实现网格单元的划分后，就可以将网格单元与GPU中的线程块进行对应，并在每个线程块内为相应网格单元中的原子开辟线程数。例如，对于具有27个网格单元的原子系统，为其分配27个线程块，然后在每个线程块内开辟线程。本实施例中由于每个网格单元内平均有64个原子，为了保证系统中每一个原子有一个单独的计算线程，因此线程块内所开辟的线程数目应大于64，本实施例中开辟128个。在上述的描述中，先对原子系统进行网格单元的划分，然后以划分后的网格单元为基础在线程块内设定线程数目。在实际应用中，也可以根据线程块中拟开辟的线程数对原子系统进行合理的网格单元的划分，但最终要能够保证一个原子对应一个单独的计算线程。

实现原子系统中的原子与线程间的对应后，就可以对各个原子所受到的作用力进行计算。在计算时，首先在GPU上执行一个内核(kernel)控制函数，该控制函数可根据多体作用模型计算原子的受力情况。由于在本实施例中，所采用的多体作用模型为Tersoff势函数，根据前面所介绍的关于Tersoff势函数的相关内容，可以得到原子系统中各个原子的势能E(可由势函数中的V_ij(r_ij)计算得到)以及整个原子系统的总体势能E_total，然后根据力计算模型对系统中原子之间的相互作用进行计算。其中，所述的力计算模型公式如下：

$F=-▿E$

表示对原子位置的求导。

但要注意的是，此时所得到的原子受力情况并不完整。Tersoff模型中的B_ij≠B_ji表示势函数在两原子之间是不对称分布的。因此系统中每个原子的受力必须单独计算，但对邻近原子的受力亦需同时进行读写，这样才能保证原子作用力计算的完整性与准确性。下面结合一个具体的实例加以说明。假设一个原子a周围有四个邻近原子，分别用b、c、d、e表示。原子a受到邻近原子b的作用力F_ab由Tersoff势函数对原子位置r_a求导得出，但与原子a对应的的计算线程只能改变原子a受到的作用力。由于力的作用是相互的，实际上原子b也应受到一个相反的力-F_ab，并需要加和到自己身上。但由于GPU采用并行多线程计算，原子b同时也正由另一个计算线程计算其与周边原子的受力，所以这个相反的力-F_ab在相同的时刻不能加和在其身上，否则易导致系统内存读写冲突。因此，此一步骤所得到的原子受力情况是不完整的。要克服这一问题首先需要记录本线程所计算的对应原子受到其所有邻近原子的全部作用力，同时将它与各个邻近原子间的作用力分别写到GPU的共享内存中，待与所有邻近原子的受力都计算完毕后，再将原子与每一个邻近原子间的作用力从共享内存写到前述的邻近原子的作用力存储数组，该存储数组在GPU全局内存中开辟。

最后，在本步骤由GPU开辟的线程实现对相应原子的完整受力进行计算。在计算时，把本线程计算原子的受力信息从全局内存读入共享内存，同时根据本线程计算的对应原子标号，把上一步骤中该原子邻近原子所在线程计算的对本原子的受力亦从全局内存中读入共享内存，并将这些作用力与上一步骤记录的原子所受到的作用力进行累加。待所有邻近原子的受力累加完毕后再把该原子的受力写回到全局内存，更新每个原子的受力信息，完成每个原子作用力的计算。将上述的累加过程在共享内存中实现有助于提高粒子作用力计算的效率，但如果对于计算效率的要求不是特别高，在具体实现时也可以将上述的累加过程直接在全局内存中实现。

在原子的完整作用力计算完毕以后，就可以根据每个原子受力的大小求取原子的加速度，再由原子的加速度以及原子的当前位置、速度求出原子在下一时刻的位置和速度，如此循环迭代计算，直至整个原子系统达到定态，此时系统中单原子平均势能在某一恒定值附近波动。

上述过程中求取原子的完整作用力以后，可由牛顿经典力学获取原子的加速度等信息，从而确定原子的运动方程：

$m\frac{\mathrm{dv}}{\mathrm{dt}}=\mathrm{\Sigma F}+F_{\mathrm{external}}$

其中，m表示原子质量；v表示原子速度；F表示原子所受邻近原子的作用力，F_external表示原子所受的外部作用力，此处不考虑其影响。对上述运动方程进行时间积分可以求出原子在下一时刻的位置和速度。

在运动方程中所涉及的与原子相关的初始信息可以从前述的原子属性信息中得到。

在上述实施例中，所要模拟的单晶硅原子体系具有1728个原子，而在另一个实施例中，在其他控制条件不变的前提下，可以扩大所要模拟的单晶硅原子体系的规模，使系统中原子个数为13824个，三个方向的空间特征长度为6.4nm。在此实施例中，在计算相同的1000个时间步下，频率为2.67GHz的CPU需要1450.2s，而利用图形处理单元GPU则只需要13.76s，系统运算效率明显提升。

在图2中给出了单晶硅原子系统在CPU环境与GPU环境中进行模拟计算所需花费时间的比较，以及不同规模的单晶硅原子系统在CPU环境与GPU环境中做模拟计算时在性能提升上的比较。其中，VC6.0环境中单核CPU计算频率为2.0GHz，gcc、icc编译环境中单核CPU计算频率为2.67GHz，其中的CPU+gcc/icc就是本发明中所提到的加载有图形处理单元GPU的计算系统。从图2的模拟计算结果可以看出，在相同的单晶硅原子系统的模拟规模下，GPU环境下模拟计算所需花费的时间要比CPU环境下计算所需花费时间更短。此外，当原子模拟系统规模较大时，GPU环境下进行模拟计算相比CPU环境性能提升更为明显。

在图3和图4中对线程块中所开辟的线程数目与原子系统所划分的网格尺寸之间的关系进行了说明，从图3中可以看出，在图形处理单元GPU上当线程块中拟开辟线程数目一定时，系统计算时间将会随着网格尺寸的增大，亦即单位网格粒子数目的增加而减少。但当线程块对应的计算网格达到一定尺寸时，计算效率不再有较大的提升，说明存在一个较优的网格尺寸或单位网格粒子数目值。

图4表明当网格尺寸恒定，亦即单位网格粒子数目恒定时，线程块中线程开辟数目亦会对计算效率产生影响。可以发现，随着线程开辟数目的增多，计算时间反而稍微升高，计算效率略微下降。这说明在一个可计算的有效范围内，线程数目尽量较少开辟，可能计算系统线程的开辟也需要微量时间。

图5、图6是基于图形处理单元GPU对上述实施例中含有1728个原子的单晶硅原子系统的模拟计算结果。从图中可以看到，计算10万个时间步亦即0.017ns时系统晶格结构依然保持完整，并且从单原子平均势能的波动和其平均值看到，计算过程是准确和有效的，这同时证实了本发明提供的粒子计算方法的可靠性。

上述说明中都以Tersoff多体作用模型为例，对如何在GPU上实现粒子系统的模拟计算进行了说明，但本领域的普通技术人员应当了解，本发明并不局限于前述的Tersoff多体作用模型，现有技术中的Tersoff-Brenner多体作用模型、Stillinger-Weber多体作用模型等都可适用于本发明。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。