CN109885917A - 一种并行分子动力学模拟方法及系统 - Google Patents

Abstract

一种并行分子动力学模拟方法，包括：依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力；依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息；基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测。本发明提供的技术方案结合邻居列表矩阵特征，确定粒子对关系矩阵和计算索引结构，实现了不同并行度下结果的按位一致，实现了计算的可复现性，同时最小限度的影响原有模拟过程，便于研究者排除模拟算法上的错误。

Description

一种并行分子动力学模拟方法及系统

技术领域

本发明涉及高性能数值计算可复现领域，具体涉及一种并行分子动力学模拟方法及系统。

背景技术

数值计算可复现性是指当科学计算在具有不同类型和数量的处理单元的同一台或不同的机器上运行多次时，获得相同的结果，是程序调试，验证和确认的关键。在排除随机数因素影响外，数值计算不可复现的原因是计算机表示精度有限，造成浮点的加法运算不具有结合律，即(a+b)+c不等于a+(b+c)。并行计算因其计算量大、并行执行模式复杂加剧了数值计算不可复现的问题。这导致判断是否将程序正确的移植到新的系统越来越困难，因为计算结果很快偏离基准案例。

由于计算机只能表示有限精度，表达的计算数据本身存在截断误差，同时由于动态调度，并行执行的影响使得计算顺序存在不确定性。而大规模并行计算的计算量大，截断误差累积次数更多，并行也加剧了计算次序的不确定性。因此，现在针对并行数值计算可复现性的解决方法主要就是从降低截断误差影响和降低计算乱序的两个方向展开的，分别对应于核心计算部分高精度数据类型和排序方法。

1.高精度数据类型

对于影响浮点运算结果的核心计算部分使用高精度数据类型是研究者大多使用的计算方法。数值可复现性和稳定性在存在噪声的模拟应用中已有一定研究。任意精度数学库在上世纪70和80年代得到了广泛的应用。在科学工程计算中，使用更高的精度类型来完成数值模拟。更高精度的计算结果可以通过任意精度库来实现，这些库可以解决很多工程问题，例如当双精度浮点数仍然无法满足精度情况下，纠错数值不稳定性。

由于并行划分的引入，导致不同并行度下针对同一问题的浮点运算顺序不一致。一般情况，用更高精度的数据类型取代原有的double数据类型，例如128bit的浮点数，但是这种方法引入了非常大的内存和计算时间消耗。折中的方法是使用大多数平台都支持的long double数据类型，一种80bit的浮点数。这种数据类型可以弱化数据类型的提高对模拟效率造成的影响。

Kahan在1965年提出了一种利用进位位(carry bit)来保证截断部分的加法运算。Kahan加将一个浮点数表征为两个部分，值和误差，通过携带进位位，来保证浮点计算的准确性。另外Knuth同样提出类似的有限精度加法，利用进位位来保证浮点计算精度。之后，对于浮点计算可复现性的研究，多集中于采用不同种类的高精度的数据类型和改进的浮点运算中。

2.排序方法

最直接的方法就是将参与计算的数据按照某一种顺序排列后相加，例如，按照升序。通过排序的方式，计算的一致性可以在任何条件下保持一致。当然，这种排序的累积和的方式也是具有截断误差的，但是其在不同的工况下的截断误差是相同的，因此可以达到计算可复现性的目的。排序方法的目的就是要确保不同并行度下，全部的浮点运算是一致的。当前研究者认为，排序方法是非常耗时的计算方法，即时在并行计算的情况下。虽然排序方法可以带来一些好处，但是在实践中对于排序的使用较少。另外一种类似排序的方式，是通过按数对相加在累加的过程。这个方法可以将相近精度的数据先相加再进行累加，这种方法的速度要优于基于排序的方法。

分子动力学模拟计算中能量偏移是典型的研究数值计算可复现的场景。单个粒子所受合力计算过程是影响统计能量偏移的核心计算过程。该过程计算耗时占据分子动力学模拟的核心部分，如果将该过程采用高精度数据类型方法，会引入更大比例的计算耗时。而单个粒子合力计算过程较容易确定计算顺序，更适合采用排序的方法设计并实现并行计算可复现性。

数值计算的不可复现问题，给实际应用中程序的移植、调试、验证和确认带来了很大的困难。实际应用中，为充分利用硬件资源和提高计算效率，数值计算通常采用并行技术。大规模，高并行计算极大地放大了数字灵敏度，从而加剧了该问题。通常采用核心计算部分高精度数据类型和排序方法两个方向来解决大规模并行计算可复现性问题。在分子动力学模拟计算中，单个粒子所受合力计算过程是影响统计能量偏移这一典型数值计算可复现场景的核心计算过程。

发明内容

为了解决现有技术中所存在的在分子动力学模拟计算中，数值计算的不可复现带来的程序的移植、调试、验证和确认困难的问题，本发明提供一种并行分子动力学模拟方法，对于单原子模拟中核心浮点计算部分合力计算，确定邻居列表粒子关系矩阵和计算索引结构，实现了不同并行度下结果的按位一致，实现了计算的可复现性，同时最小限度的影响原有模拟过程，便于研究者排除模拟算法上的错误。

本发明提供的技术方案是：

一种并行分子动力学模拟方法，其改进之处在于，包括：

依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力；

依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息；

基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测。

优选地，所述全邻居列表粒子关系矩阵的建立包括：

基于所有粒子坐标，速度、受力关系和作用势建立系统中所有粒子的全邻居列表粒子关系矩阵。

优选地，建立所述全邻居列表粒子关系矩阵还包括确定矩阵中各元素的值：

若两个粒子间存在作用关系，则两个粒子在所述全邻居列表粒子关系矩阵中对应的元素的值为1；

若两个粒子间不存在作用关系，则两个粒子在所述全邻居列表粒子关系矩阵中对应的元素的值为0。

优选地，所述依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力包括：

依据所述全邻居列表粒子关系矩阵的顺序，利用所述矩阵中第i行元素计算系统中第i个粒子所受合力，所述i表示正整数，直至计算系统中所有粒子所受的合力。

优选地，所述利用所述矩阵中第i行元素计算系统中第i个粒子所受合力包括：

以粒子间距离或列表粒子序作为标识，建立计算索引结构；

依据所述计算索引结构的顺序，利用所述矩阵中第i行元素计算系统中第i个粒子所受合力。

优选地，所述第i个粒子所受合力按下式进行计算：

式中，N：系统中与粒子i存在作用势的粒子的个数；h(i,j)：粒子i和粒子j间的关系函数；f(i,j)：粒子i和粒子j间的作用势函数。

优选地，所述粒子i和粒子j间的关系函数h(i,j)如下式所示：

所述粒子i和粒子j间的作用势函数f(i,j)如下式所示：

式中：r：粒子i和粒子j间的距离。

优选地，所述依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息包括：

依据所述每个粒子所受合力，通过同一个随机数发生器生成相同的速度生成命令随机数，利用Verlet速度积分算法计算和更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息。

优选地，所述基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测包括：

利用调整函数对每个粒子的速度和位置信息进行排序，确定累加统计顺序；

基于所述累加统计顺序，进行系统宏观物理量的统计与观测。

优选地，所述宏观物理量包括：动能、温度、势能和守恒量。

一种并行分子动力学模拟系统，包括：合力模块、更新模块、统计与观测模块；

合力模块：用于依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力；

更新模块：用于依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息；

统计与观测模块：用于基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测。

与现有技术相比，本发明的有益效果为：

本发明提供的技术方案采用全邻居列表粒子关系矩阵实现了不同并行度下结果的按位一致，实现单原子模拟情况下并行计算可复现性，同时最小限度的影响原有模拟过程，便于排除模拟算法上的错误。

本发明提供的技术方案在通过对宏观物理量的统计、观测和粒子所受合力的计算过程中采用定序操作的方法，在规避采用高精度数据类型方法引入更大比例的计算耗时的情况下，确保数值计算可复现性，最大程度还原模拟过程，便于程序的移植、调试、验证和确认。

本发明提供的技术方案适用于多种分子动力学软件，实现过程简单、易理解，适用性强。

附图说明

图1为本发明的并行分子动力学模拟方法示意图；

图2为本发明的分子动力学模拟基本流程；

图3为本发明S粒子分布情况；

图4为本发明半邻居列表粒子关系矩阵；

图5为本发明全邻居列表粒子关系矩阵；

图6为本发明并行分子动力学模拟系统示意图。

具体实施方式

为了更好地理解本发明，下面结合说明书附图和实例对本发明的内容做进一步的说明。

实施例1：

一种并行分子动力学模拟方法，如图1所示，包括：

步骤1：依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力；

步骤2：依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息；

步骤3：基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测。

步骤1：依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力。

具体地，所述全邻居列表粒子关系矩阵的建立包括：

基于所有粒子坐标，速度、受力关系和作用势建立系统中所有粒子的全邻居列表粒子关系矩阵。

具体地，所述全邻居列表粒子关系矩阵的建立还包括确定矩阵中各元素的值；

若两个粒子间存在作用关系，则两个粒子在所述全邻居列表粒子关系矩阵中对应的元素的值为1；

若两个粒子间不存在作用关系，则两个粒子在所述全邻居列表粒子关系矩阵中对应的元素的值为0。

具体地，所述依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力包括：

依据所述全邻居列表粒子关系矩阵的顺序，利用所述矩阵中第i行元素计算系统中第i个粒子所受合力，所述i表示正整数，直至计算系统中所有粒子所受的合力。

具体地，所述利用所述矩阵中第i行元素计算系统中第i个粒子所受合力包括：

以粒子间距离或列表粒子序作为标识，建立计算索引结构；

依据所述计算索引结构的顺序，利用所述矩阵中第i行数据按粒子间距离序伪代码计算系统中第i个粒子所受合力。

具体地，所述第i个粒子所受合力按下式进行计算：

式中，N：系统中与粒子i存在作用势的粒子的个数；h(i,j)：粒子i和粒子j间的关系函数；f(i,j)：粒子i和粒子j间的作用势函数。

具体地，所述粒子i和粒子j间的关系函数h(i,j)如下式所示：

所述粒子i和粒子j间的作用势函数f(i,j)如下式所示：

式中：r：粒子i和粒子j间的距离。

步骤2：依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息；

具体地，所述依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息包括：

依据所述每个粒子所受合力，通过同一个随机数发生器生成相同的速度生成命令随机数，利用Verlet速度积分算法计算和更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息。

步骤3：基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测。

具体地，所述基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测包括：

利用调整函数对每个粒子的速度和位置信息进行排序，确定累加统计顺序；

基于所述累加统计顺序，进行系统宏观物理量的统计与观测。

具体地，所述宏观物理量包括：动能、温度、势能和守恒量。

实施例2：

基于同一种发明构思，本发明还提供了一种并行分子动力学模拟系统，如图6所示，包括：合力模块、更新模块、统计与观测模块；

合力模块：用于依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力；

更新模块：用于依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息；

统计与观测模块：用于基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测。

其中，合力模块中，具体地，所述全邻居列表粒子关系矩阵的建立包括：

基于所有粒子坐标，速度、受力关系和作用势建立系统中所有粒子的全邻居列表粒子关系矩阵。

具体地，所述全邻居列表粒子关系矩阵的建立还包括确定矩阵中各元素的值；

若两个粒子间存在作用关系，则两个粒子在所述全邻居列表粒子关系矩阵中对应的元素的值为1；

若两个粒子间不存在作用关系，则两个粒子在所述全邻居列表粒子关系矩阵中对应的元素的值为0。

具体地，所述依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力包括：

依据所述全邻居列表粒子关系矩阵的顺序，利用所述矩阵中第i行元素计算系统中第i个粒子所受合力，所述i表示正整数，直至计算系统中所有粒子所受的合力。

具体地，所述利用所述矩阵中第i行元素计算系统中第i个粒子所受合力包括：

以粒子间距离或列表粒子序作为标识，建立计算索引结构；

依据所述计算索引结构的顺序，利用所述矩阵中第i行数据按粒子间距离序伪代码计算系统中第i个粒子所受合力。

具体地，所述第i个粒子所受合力按下式进行计算：

式中，N：系统中与粒子i存在作用势的粒子的个数；h(i,j)：粒子i和粒子j间的关系函数；f(i,j)：粒子i和粒子j间的作用势函数。

具体地，所述粒子i和粒子j间的关系函数h(i,j)如下式所示：

所述粒子i和粒子j间的作用势函数f(i,j)如下式所示：

式中：r：粒子i和粒子j间的距离。

更新模块中，所述依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息包括：

依据所述每个粒子所受合力，通过同一个随机数发生器生成相同的速度生成命令随机数，利用Verlet速度积分算法计算和更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息。

统计与观测模块中，所述基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测包括：利用调整函数对每个粒子的速度和位置信息进行排序，确定累加统计顺序；

基于所述累加统计顺序，进行系统宏观物理量的统计与观测。

具体地，所述宏观物理量包括：动能、温度、势能和守恒量。

实施例3：

经典的分子动力学并行划分方法中，分为力分解方法、原子分解方法和空间分解方法。本发明提出的技术方案是基于空间分解的分子动力学划分方法。

利用分子动力学求解，原子之间的作用势场为：

式中，r：原子I和J之间的距离，r_c：作用势的截断半径。ε：特征能量；和σ：特征长度。

选取液态Ar作为模型流体，其中σ＝0.34nm，ε＝1.67×10^-21J，m＝6.63×10^-26kg，式中，m：原子质量。为了节省计算时间，截断半径r_c＝2.2σ。

在一个封闭环境中模拟，具有恒定的分子数目，恒定的模拟体积和恒定的能量的微正则系统是有统计力学参与的分子动力学模拟系统。

分子动力学模拟的基本流程如图2所示，选取要研究的系统及其边界，选取系统内粒子间的作用势能模型，设定系统中粒子的初始位置和初始动量。进而通过相应积分算法和势能函数，计算粒子间作用力及各粒子的速度和位置。最后当体系达到平衡后，依据相关的统计公式，获得各宏观参数和输运性质，用于模拟数据结果分析。

分子动力学模拟已有很多可用软件，具有代表性的是LAMMPS(large-scaleatomistic/molecular massively parallel simulator)，它是C++语言编写的一款基于空间分解划分方法的开源软件包。在并行执行中，待模拟区域通过空间划分的方法，将子区域分配到各个处理器上，处理器之间的通信通过调用MPI来完成，可以获得达到90％甚至更高的并行效率。针对分子动力学基础计算流程，在LAMMPS中按以下方式实现：

步骤1：依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力。

通过create_atom、velocity和pair命令进行初始化。create_atom可以通过给定模拟系统密度，模拟区域大小，晶格类型直接对系统粒子进行初始化，也可以通过读取已有数据文件对系统进行初始化，对于给定模拟系统中每个粒子可给予一个Tag标识作为粒子标识。该方法在不同并行度下，按照MPI_Cart分配的处理器编号不同，造成不同并行度下对粒子的编号方式不同。Velocity命令根据指定速度分布和随机数对系统中粒子速度进行初始化，可以通过读取文件初始化。pair命令对于系统中粒子对间的作用势进行设置。

半列表，利用牛顿第三定律，位于同一处理器每个粒子对(i，j)仅在列表中出现一次，仅被某一个处理器存储；

全列表，位于同一处理器每个粒子对在列表中出现两次，被i和j分别存储，默认的邻居列表的建立方式为半邻居列表。

在并行模拟情况，每个处理器维护本地的邻居列表结构。每个本地邻居列表的合集是该模拟系统的总邻居列表，但值得注意的是，半邻居列表的合集在不同并行度下是不等价的，而全邻居列表是等价的。对于解决并行可复现的方法主流分为高精度和排序方法，结合以上分析，解决合力计算部分的计算可复现方法可以综合为如下四种方法：

表1合力计算解决方案表

编号	方法名称
		M1	半列表+高精度
M2	全列表+高精度
		M3	半列表+排序
M4	全列表+排序

在分子动力学的模拟计算中，核心需要维护的数据结构就是粒子间关系的邻居列表。一个粒子i的邻居列表维护了与其发生作用势的其他粒子j的标号。由于粒子间作用势在3D情况下是维护球形区域的作用关系。根据粒子间作用势函数和计算模拟需求，一般会设置邻居列表维护的截断半径，从而确定列表中的j粒子。

那么，对于如图3所示的给定的N个粒子的模拟系统S来说，如果某一个粒子j在粒子i的截断半径作用势范围内，则认为j和i是有作用势关系的，j和i称为一个粒子作用势对。对于给定系统S中的全部粒子作用势对的合集，称为系统S的粒子关系矩阵M。基于牛顿第三定律，半邻居粒子关系矩阵如图4所示，全邻居列表料子关系矩阵如图5所示。

对于给定模拟系统的粒子位置初始分布，可以通过create_atom命令生成。该命令在不同并行度下，因为按照MPI_Cart分配的处理器编号不同，造成不同并行度下对粒子的编号方式不同。虽然该编号方式Tag并不是分子动力学模拟所必须的，但是在单原子的模拟中，该Tag对于计算可复现性的影响非常大，因此必须保证该Tag的编号方式要保持一致。针对LAMMPS本身的原子生成方法(create atom)，通过LAMMPS原有的write_data和read_data功能，通过预先Profile生成某模拟规模下的初始原子的位置和速度信息，将这些信息通过write_data写入文件中。在完成不同并行度的模拟中，通过read_data来读取这些信息，这样能够保证在任意并行度下，粒子的编号顺序是一致的。这是针对单原子并行计算可复现的第一个重要步骤。

在Velocity命令中，其中LOOP关键字默认是ALL，但是如果为LOCAL时，每个处理器会使用本地的随机数；因此可以通过固定随机数发生器的方法，可以排除速度生成命令随机数不同带来的影响。

粒子间作用势计算：根据pair命令设定的作用势函数进行粒子间作用势。对于每个粒子，考察位于其作用势截断半径内其他粒子对其的作用势和受力，并计算每个粒子所受合力。

每个粒子的合力计算是针对粒子关系矩阵的一行进行计算。在该行中的粒子通过粒子间作用势函数进行浮点运算，然后累加求和得到每个粒子所受合力。

r为i与j之间的距离

针对基于空间分解划分的分子动力学模拟设计的算法，对于给定模拟系统，首先需要建立该系统中粒子的全邻居列表，如图5所示，对于所有存在作用势关系的粒子作用势对都进行存储，保证参与计算的数据一致。同时，对于得到的邻居列表关系矩阵，我们需要给出确定计算顺序的指标，得到计算索引结构。通过计算索引结构，对关系矩阵中每一行参与计算的数据进行定序合力计算，从而达到在任一并行度下的计算可复现性。

针对上述合力计算过程伪码，计算合力分为对力分量累加定序和内循环定序两个层次：

第一层次：对于力分量累加定序方式最为简单直接，直接在累加操作上，对力的三个分量累加过程，对于参与累加的数据，通过排序方式直接排序。每个时间步都需要对三个分量进行排序，相对耗时。

第二层次：对于内循环的定序方式是结合邻居列表的组织结构，需要建立计算索引结构。基于单粒子计算模拟的特征量决定，通常情况下以粒子对间的距离作为特征量的标识。虽然降低了排序次数，但是索引结构的引入带来额外的内存消耗。

步骤2：依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息。

粒子速度和位置信息更新：采用的积分算法Verlet Velocity积分算法，在一个计算时间步内，将速度的更新分为两个部分，首先更新全部的x和一半的v，继而根据新的位置信息x计算原子的合力信息，进而根据新的合力部分，更新后半部的v。LAMMPS中对于物理属性的维护都是按照x，y，z三个方向分别进行计算和更新，如位置信息(x,y,z)、速度信息(vx,vy,vz)和受力信息(fx,fy,fz)。

步骤3：基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测。

从微观的粒子数据信息进过系综平均，在给定时间步统计宏观数据信息，用于整个模拟系统的演化过程的观测。跟踪模拟问题需要观察模拟的实时性结果，以及通过数据处理观察是否问题达到要求的效果。这种系综平均的核心操作是MPI_Allreduce操作，是一种累加求和的操作。这类的物理属性包括：动能Ke，温度T，势能P。模拟系统的动能和势能的和就是整个模拟系统的能量。

在具有恒定的分子数目，恒定的模拟体积和恒定的能量的微正则系统分子动力学模拟中，外界对模拟系统的能量不产生影响，因此在长时间的观察态情况下，系统的能量值应该与初始的能量值发生很小的波动。但是在不同的并行度模拟中，由于浮点累加操作，造成统计的物理属性值发生了不同。为了解决这个方面的问题，我们在宏观属性统计值的部分，新增加了调整函数，使得无论在何种并行度条件下，统计的宏观属性的计算顺序都是一致的，这也是在第一个重要步骤完成的基础上进行的。调整函数的主要操作步骤如下：

对任意大于1并行度的宏观统计结果与串行宏观统计结果之间的误差，是该计算方法下的并行计算误差，对于一个分子力学模拟计算方法，宏观统计结果为m个，这m个串行和并行计算结果之间差值，求取期望得到并行度对计算结果造成的影响，计算方法如下所示：

以具有周期边界的立方体中充斥着液态Ar流体的模拟为例，其参数为σ＝0.34nm，ε＝1.67×10^-21J和m＝6.63×10^-26kg，m：原子质量。如果采用传统合力计算方法，其并行计算结果与串行计算结果之间的误差约为0.0022；而采用本专利所提出的基于邻居列表特征方法，可以实现任何并行度下的计算可复现性。

由上述实施例可以看出，本发明提供的技术方案采用全邻居列表粒子关系矩阵实现了不同并行度下结果的按位一致，实现单原子模拟情况下并行计算可复现性，同时最小限度的影响原有模拟过程，便于排除模拟算法上的错误。

本发明提供的技术方案在通过对宏观物理量的统计、观测和粒子所受合力的计算过程中采用定序操作的方法，在规避采用高精度数据类型方法引入更大比例的计算耗时的情况下，确保数值计算可复现性，最大程度还原模拟过程，便于程序的移植、调试、验证和确认。

本发明提供的技术方案适用于多种分子动力学软件，实现过程简单、易理解，适用性强。

显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上仅为本发明的实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均包含在申请待批的本发明的权利要求范围之内。

Claims (11)

1.一种并行分子动力学模拟方法，其特征在于，包括：

依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力；

依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息；

基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测。

2.如权利要求1所述的并行分子动力学模拟方法，其特征在于，所述全邻居列表粒子关系矩阵的建立包括：

基于所有粒子坐标，速度、受力关系和作用势建立系统中所有粒子的全邻居列表粒子关系矩阵。

3.如权利要求2所述的并行分子动力学模拟方法，其特征在于，所述全邻居列表粒子关系矩阵的建立还包括确定矩阵中各元素的值；

若两个粒子间存在作用关系，则两个粒子在所述全邻居列表粒子关系矩阵中对应的元素的值为1；

若两个粒子间不存在作用关系，则两个粒子在所述全邻居列表粒子关系矩阵中对应的元素的值为0。

4.如权利要求1所述的并行分子动力学模拟方法，其特征在于，所述依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力包括：

依据所述全邻居列表粒子关系矩阵的顺序，利用所述矩阵中第i行元素计算系统中第i个粒子所受合力，所述i表示正整数，直至计算系统中所有粒子所受的合力。

5.如权利要求4所述的并行分子动力学模拟方法，其特征在于，所述利用所述矩阵中第i行元素计算系统中第i个粒子所受合力包括：

以粒子间距离或列表粒子序作为标识，建立计算索引结构；

依据所述计算索引结构的顺序，利用所述矩阵中第i行元素计算系统中第i个粒子所受合力。

6.如权利要求5所述的并行分子动力学模拟方法，其特征在于，所述第i个粒子所受合力按下式进行计算：

式中，N：系统中与粒子i存在作用势的粒子的个数；h(i,j)：粒子i和粒子j间的关系函数；f(i,j)：粒子i和粒子j间的作用势函数。

7.如权利要求6所述的并行分子动力学模拟方法，其特征在于，所述粒子i和粒子j间的关系函数h(i,j)如下式所示：

所述粒子i和粒子j间的作用势函数f(i,j)如下式所示：

式中：r：粒子i和粒子j间的距离。

8.如权利要求1所述的并行分子动力学模拟方法，其特征在于，所述依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息包括：

依据所述每个粒子所受合力，通过同一个随机数发生器生成相同的速度生成命令随机数，利用Verlet速度积分算法计算和更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息。

9.如权利要求1所述的并行分子动力学模拟方法，其特征在于，所述基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测包括：

利用调整函数对每个粒子的速度和位置信息进行排序，确定累加统计顺序；

基于所述累加统计顺序，进行系统宏观物理量的统计与观测。

10.如权利要求9所述的并行分子动力学模拟方法，其特征在于，所述宏观物理量包括：动能、温度、势能和守恒量。

11.一种并行分子动力学模拟系统，其特征在于，所述系统包括：合力模块、更新模块、统计与观测模块；

合力模块：用于依据预先建立的全邻居列表粒子关系矩阵的顺序计算每个粒子所受合力；

更新模块：用于依据所述每个粒子所受合力计算并更新所述全邻居列表粒子关系矩阵中所有粒子的速度和位置信息；

统计与观测模块：用于基于所有粒子的速度和位置信息对系统宏观物理量进行统计与观测。