CN102411658A - 一种基于cpu和gpu协作的分子动力学加速方法 - Google Patents

Abstract

本发明公开了一种基于CPU和GPU协作的分子动力学加速方法，目的是提出一种通过协同使用CPU和GPU的加速方法，在较低成本的情况下实现分子动力学的加速。技术方案是使用配置有多核CPU和GPU的计算机，对待模拟的分子系统建立“网格-单元”结构存储分子信息，使用CPU和GPU动态申请和处理网格，在处理网格时，CPU以单元为单位处理网格，GPU则直接处理整个网格，实现CPU和GPU高效协作地完成分子动力学模拟。采用本发明可使得CPU核之间负载均衡，避免静态分配计算任务时CPU和GPU可能出现闲置等待的情况，CPU和GPU都能充分发挥自己的计算性能，提高了整个计算机系统的利用率，以低成本实现了分子动力学的加速。

Description

一种基于CPU和GPU协作的分子动力学加速方法

技术领域

本发明涉及一种分子动力学加速方法，尤指基于CPU和GPU协作的分子动力学加速方法。

背景技术

分子动力学方法是一种依靠牛顿力学定律来模拟分子体系的运动，从而计算分子体系宏观性质的方法。它是研究纳米尺度物理现象的重要手段，广泛应用于材料科学、生物物理和药物设计等领域。随着模拟精度的提高和模拟规模的增大，分子动力学模拟时间也随之增长，从而对其实际应用产生了很大影响。目前通常的做法有以下几种：

1、通过组建集群系统来并行的进行模拟。但这种方法硬件成本较高，同时也需要开发人员不断改进并行算法来挖掘集群系统性能。因此在考虑硬件、人力成本时，不适宜采用该方法。

2、通过定制专用计算机来加速模拟。针对分子动力学模拟的不同阶段定制相应的专用硬件进行计算，从而达到很高的模拟速度。但是这种机器硬件成本很高，目前该方法很少采用。

3、利用通用加速单元GPU(Graphics Processing Unit，图形处理单元)来协作普通计算机进行模拟。目前GPU厂商为用户提供了通用编程模型如CUDA(统一计算架构模型)等对GPU进行编程，利用GPU浮点计算能力比CPU强的特点，使用GPU完成分子动力学中计算任务，从而提高模拟速度。但是该方法目前存在一些局限性，主要表现在使用GPU计算时，计算机中CPU处于闲置状态，计算机的总体利用率较低。

第三种方法由于成本较低，是目前前景较大的分子动力学加速方法。但是目前该方法仍存在一些问题：一方面CPU和GPU的计算模型不同，GPU适合一次性处理大量数据，因此它要求分配的计算任务不能太小，而大块的计算任务又不利于CPU进行调度，为了使CPU和GPU都能高效计算，必须设计大小适中的计算任务；另一方面CPU和GPU完成计算任务的时间难以估量，静态地分配计算任务会造成CPU或GPU出现闲置等待的情况，应根据两者的忙闲动态调整计算任务的分配，以保证两者皆不出现闲置等待的情况。但目前还没有公开的技术方案能较好解决这两方面的问题。

发明内容

本发明要解决的技术问题在于：针对当前基于GPU的分子动力学加速方法效率较低的问题，提出一种通过协同使用CPU和GPU的加速方法，使系统中CPU和GPU都能高效计算，从而在保持较低成本的情况下实现分子动力学的加速。

为了解决上述技术问题，本发明提出的技术方案为：使用配置有多核CPU和GPU的计算机，对待模拟的分子系统建立“网格-单元”结构存储分子信息，使用CPU和GPU动态的申请和处理网格，在处理网格时，CPU以单元为单位处理网格，GPU则直接处理整个网格，从而实现CPU和GPU高效协作地完成分子动力学模拟。

具体技术方案为：

第一步、构建配置有多核CPU和GPU的计算机系统。选择包含Q个(Q＞＝4)核的CPU作为计算机系统的处理器，选择支持CUDA的GPU通过PCI-E接口与CPU相连。CPU的Q个核顺序编号为0，...，Q-1。

第二步、从存储于硬盘的文件中读入分子系统相关参数，包括截断半径R_c和所有分子的数据量。其中截断半径R_c表示每个分子与其他分子发生作用的最大距离，如果两个分子之间的距离超过R_c，则它们之间不发生作用。分子数据量包括分子总数目MN以及每个分子的编号、位置向量和速度向量。在CPU的内存中申请二维双精度浮点型数组Pos存储分子的位置向量，申请二维双精度浮点型数组Vel存储分子的速度向量。

第三步、设置模拟参数，包括目标网格权重TW，目标模拟步数T和实际模拟步数t。TW和T的值由模拟需求确定，TW的值为期望网格包含的分子数目，T的值为需要模拟的步数，实际模拟步数t初始化为0。

第四步、使用“网格-单元”结构表示分子系统，“网格-单元”结构的建立过程如下：

4.1使用“单元”结构表示分子系统，步骤为：

4.1.1将分子系统划分为均匀组合的边长为R_c的n个正方体区域，称这些正方体区域为单元，将这些单元顺序编号为0，...，n-1。

4.1.2在CPU的内存中初始化结构体数组SCell保存单元中的分子编号以及单元的邻接单元编号，称该结构体数组为“单元”结构，使用符号SCell[k](0≤k＜n)表示序号为k的单元。每个“单元”拥有一个变量成员和两个数组成员，分别为单元内分子数目AN、单元内分子编号数组AIdx和邻接单元编号数组NB。AIdx顺序存储单元中分子编号。NB顺序存储单元在三维空间中邻接26个单元的编号。用符号SCell[k].AN表示序号为k的单元包含的分子数目；用符号SCell[k].AIdx表示序号为k的单元中分子编号数组；用符号SCell[k].NB表示序号为k的单元邻接单元编号数组。

4.2根据单元内分子数目SCell[k].AN(0≤k＜n)将多个单元组成区域，使每个区域中的分子数目相等，称这些区域为网格。在CPU的内存中用链表保存网格中分子信息以及网格处理标识，链表的每一个节点称为一个“网格”结构。每个“网格”结构拥有一个指针、一个整型变量、一个布尔变量、一个整型数组和四个双精度浮点型数组，分别为后续“网格”指针Next，网格内单元数目CN，网格处理标识R，网格内单元编号数组CIdx，网格内分子位置信息数组Ppos，网格内分子速度信息数组Pvel，网格内分子受力信息数组Pforce，网格外分子位置信息数组Gpos。为了便于描述，用SPatch表示指向链表首节点的指针，称这个节点为SPatch指向的网格；用符号“→”引用节点中的元素。例如，用SPatch→Next表示指向下一个网格的指针；用SPatch→CN表示SPatch指向的网格中单元数目；用SPatch→R表示SPatch指向的网格处理标识，若SPatch→R的值为True则表示SPatch指向的网格已经被处理，若SPatch→R的值为False则表示SPatch指向的网格未被处理；用SPatch→CIdx表示SPatch指向的网格中单元编号数组，用SPatch→CIdx[x](0≤x＜SPatch→CN)存储网格中序号为x的单元编号；用SPatch→Ppos表示SPatch指向的网格内分子位置信息数组；用SPatch→Pvel表示SPatch指向的网格内分子速度信息数组；用SPatch→Pforce表示SPatch指向的网格内分子受力信息数组；用SPatch→Gpos表示SPatch指向的网格外分子位置信息数组。生成网格结构的步骤为：

4.2.1计算每个单元的权重W₀，W₁，...，W_n-1，其中W_k(0≤k＜n)的值为单元k中分子数目SCell[k].AN。初始化当前网格负载PW为0，起始单元编号S_t为0，结束单元编号E_d为0。初始化指针CP，使指针CP＝SPatch。

4.2.2将指向下一个网格的指针CP→Next初始化为NULL，CP指向的网格中单元数目CP→CN初始化为0，CP指向的网格处理标识CP→R初始化为False，CP指向的数组CP→CIdx、CP→Ppos、CP→Pvel、CP→Pforce、CP→Gpos全部初始化为0。

4.2.3判断PW是否小于目标网格权重TW，若是，执行4.2.4，否则转4.2.6。

4.2.4判断E_d是否等于n-1，若是，转4.2.6，否则执行4.2.5。

4.2.5

E_d＝E_d+1，转4.2.3。

4.2.6将编号S_t，S_t+1，...，E_d顺序填入CP →CIdx数组中。根据SCell[a].AIdx(S_t≤a＜E_d)中的分子编号将Pos和Vel数组对应分子的位置向量和速度向量顺序填入CP→Ppos和CP→Pvel中。顺序检查SCell[a].AIdx(S_t≤a＜E_d)数组中的单元编号b，若b＜S_t或b＞E_d，则按SCell[b].AIdx中的分子编号将Pos和Vel数组对应分子的位置向量顺序填入CP→Gpos中。将CP→CN赋值为E_d-S_t+1。判断E_d是否等于n-1，若是，转第五步，否则执行4.2.7。

4.2.7使PW＝0，S_t＝E_d+1，E_d＝E_d+1，CP＝CP→Next。转4.2.2。

第五步、启动两个线程控制GPU和CPU并行的处理网格，称控制GPU的线程为0号线程，控制CPU的线程为1号线程，分别运行在CPU编号为0和1的核上。并行处理网格的方法为：0号线程初始化指针N₀＝SPatch，将N₀→R置为True，执行第六步；1号线程初始化指针N₁＝SPatch→Next，将N₁→R置为True，转第九步。

第六步、0号线程控制GPU对N₀指向的网格进行处理，具体做法为：

6.1将N₀指向的网格的数据拷贝到GPU的内存中。

6.2启动GPU线程采用牛顿力学方程计算N₀中每个分子所受的力，根据力和初始速度计算每个分子当前速度和位置，将计算结果存放在GPU的内存中。

6.3将计算结果从GPU的内存拷贝到位于CPU的内存中的N₀→Ppos和N₀→Pvel数组中。

第七步判断N₀→Next是否为NULL，若是，转第十二步，否则执行第八步。

第八步、使N₀＝N₀→Next，判断N₀→R是否为True，若是转第七步，否则将N₀→R置为True并转第六步。

第九步、1号线程根据CPU的核数Q启动Q-2个子线程T₀，T₁，...，T_Q-3运行于CPU的编号为2，3，...，Q-1的核上，T_c(0≤c＜Q-3)采用牛顿力学方程计算单元

中分子所受的力(其中符号

代表上取整)，T_Q-3采用牛顿力学方程计算单元

...，N₁→CIdx[(N₁→CN)-1]中分子所受的力，T₀，T₁，...，T_Q-3将计算所得的力存放到N₁→Pforce数组中，根据力和初始速度计算每个分子当前速度和位置，将计算结果存放到N₁→Ppos和N₁→Pvel数组中。

第十步、判断N₁→Next是否为NULL，若是，则转第十二步，否则执行第十一步。

第十一步、N₁＝N₁→Next，判断N₁→R是否为True，若是转第十步，否则置N₁→R为True并转第九步。

第十二步、等待0号线程和1号线程都运行到第十二步后，执行第十三步。

第十三步、使t＝t+1，判断t是否小于等于T，是则将并转第五步。否则转第十四步。

第十四步、将计算结果包括分子系统中所有分子的位置向量和速度向量写入位于硬盘的文件中，结束。

与现有技术相比，采用本发明可达到以下技术效果：

1.第五步到第十一步以网格为单位动态地分配计算任务给CPU和GPU进行处理，当CPU或者GPU处理完一个计算任务时，直接申请下一个计算任务，从而避免了静态分配计算任务时CPU和GPU可能出现闲置等待的情况，提高了整个计算机系统的利用率。

2.第四步建立了“网格-单元”结构，并以网格为单位将计算任务分配给CPU和GPU。一个网格的计算任务可以满足GPU一次性处理大量计算任务的需求；而当网格分配给CPU时，CPU在第七步以单元为单位进行计算任务在CPU核上的调度，满足了CPU核之间负载均衡的要求。因此，CPU和GPU都能充分发挥自己的计算性能。

3.由于在整个过程中仅需要用到配置多核CPU和GPU的计算机系统，而这两种部件在目前个人PC上极为常见，因此以低成本实现了分子动力学的加速。

附图说明

图1是本发明的总流程图。

具体实施方案

图1是本发明的总流程图。

步骤1)、构建配置有Q核CPU和GPU的计算机系统；

步骤2)、从硬盘中读入分子系统参数包括截断半径R_c和分子数据量；

步骤3)、设置模拟参数包括目标网格权重TW、目标模拟步数T和实际模拟步数t，初始化实际模拟步数t为0。

步骤4)、对目标分子系统建立“网格-单元”结构，包括结构体数组SCell和链表SPatch；

步骤5)、启动0号线程和1号线程控制GPU和CPU进行计算，将N₀指向SPatch，将N₁指向SPatch→Next，置N₀→R为True，置N₁→R为True，0号线程转步骤6)，1号线程转步骤9)；

步骤6)、0号线程控制GPU根据牛顿力学方程对N₀指向的网格进行计算。

步骤7)、判断N₀→Next是否为NULL，若是，则转步骤12)，否则使N₀＝N₀→Next，并执行步骤8)；

步骤8)、判断N₀→R是否为True，若是则转步骤7)，否则将N₀→R置为True，转步骤6)；

步骤9)、1号线程启动Q-2个子线程以单元为单位根据牛顿力学方程对N₁指向的网格进行计算。

步骤10)、判断N₁→Next是否为NULL，若是，则转步骤12)，否则执行步骤11)；

步骤11)、使N₁＝N₁→Next，判断N₁→R是否为True，若是则转步骤10)，否则将N₁→R置为True，转步骤9)；

步骤12)、等待线程0和线程1都执行到此步；

步骤13)、使t＝t+1，判断t是否小于等于T，若是，则转步骤5)，否则执行步骤14)；

步骤14)、将计算结果写入到硬盘的文件中，结束。

Claims (2)

1.一种基于CPU和GPU协作的分子动力学加速方法，其特征在于包括以下步骤：

第一步、构建配置有Q核CPU和GPU的计算机系统，Q＞＝4，采用的GPU支持CUDA，通过PCI-E接口与CPU相连，CPU的Q个核顺序编号为0，...，Q_-1；

第二步、从存储于硬盘的文件中读入分子系统相关参数，包括截断半径R_c和所有分子的数据量，其中截断半径R_c表示每个分子与其他分子发生作用的最大距离，分子数据量包括分子总数目MN以及每个分子的编号、位置向量和速度向量；在CPU的内存中申请二维双精度浮点型数组Pos存储分子的位置向量，申请二维双精度浮点型数组Vel存储分子的速度向量；

第三步、设置模拟参数，包括目标网格权重TW，目标模拟步数T和实际模拟步数t；

第四步、使用“网格-单元”结构表示分子系统，“网格-单元”结构的建立过程如下：

4.1使用“单元”结构表示分子系统，步骤为：

4.1.1将分子系统划分为均匀组合的边长为R_c的n个正方体区域，称这些正方体区域为单元，将这些单元顺序编号为0，...，n-1；

4.1.2在CPU的内存中初始化结构体数组SCell保存单元中的分子编号以及单元的邻接单元编号，称该结构体数组为“单元”结构，使用符号SCell[k]表示序号为k的单元，0≤k＜n；每个“单元”拥有一个变量成员和两个数组成员，分别为单元内分子数目AN、单元内分子编号数组Aldx和邻接单元编号数组NB；Aldx顺序存储单元中分子编号，NB顺序存储单元在三维空间中邻接26个单元的编号；用符号SCell[k].AN表示序号为k的单元包含的分子数目，用符号SCell[k].AIdx表示序号为k的单元中分子编号数组，用符号SCell[k].NB(表示序号为k的单元邻接单元编号数组；

4.2根据单元内分子数目SCell[k].AN将多个单元组成区域，称这些区域为网格，生成网格结构的步骤为：

4.2.1计算每个单元的权重W₀，W₁，...，W_n+1，其中W_k的值为单元k中分子数目SCell[k].AN，初始化当前网格负载PW为0，起始单元编号S_t为0，结束单元编号E_d为0，初始化指针CP，使指针CP＝SPatch；在CPU的内存中用链表保存网格中分子信息以及网格处理标识，链表的每一个节点称为一个“网格”结构；每个“网格”结构拥有一个指针、一个整型变量、一个布尔变量、一个整型数组和四个双精度浮点型数组，分别为后续“网格”指针Next，网格内单元数目CN，网格处理标识R，网格内单元编号数组CIdx，网格内分子位置信息数组Ppos，网格内分子速度信息数组Pvel，网格内分子受力信息数组Pforce，网格外分子位置信息数组Gpos；SPatch表示指向链表首节点的指针，称这个节点为SPatch指向的网格；用符号“→”引用节点中的元素；

4.2.2将指向下一个网格的指针CP→Next初始化为NULL，CP指向的网格中单元数目CP→CN初始化为0，CP指向的网格处理标识CP→R初始化为False，CP指向的数组CP→CIdx、CP→Ppos、CP→Pvel、CP→Pforce、CP→Gpos全部初始化为0；

4.2.3判断PW是否小于目标网格权重TW，若是，执行4.2.4，否则转4.2.6；

4.2.4判断E_d是否等于n-1，若是，转4.2.6，否则执行4.2.5；

4.2.5

E_d＝E_d+1，转4.2.3；

4.2.6将编号S_t，S_t+1，...，E_d顺序填入CP→CIdx数组中，根据SCell[a].AIdx中的分子编号将Pos和Vel数组对应分子的位置向量和速度向量顺序填入CP→Ppos和CP→Pvel中，S_t≤a＜E_d；顺序检查SCell[a].AIdx数组中的单元编号b，若b＜S_t或b＞E_d，则按SCell[b].AIdx中的分子编号将Pos和Vel数组对应分子的位置向量顺序填入CP→Gpos中，将CP→CN赋值为E_d-S_t+1；判断E_d是否等于n-1，若是，转第五步，否则执行4.2.7；

4.2.7使PW＝0，S_t＝E_d+1，E_d＝E_d+1，CP＝CP→Next，转4.2.2；

第五步、启动两个线程控制GPU和CPU并行的处理网格，称控制GPU的线程为0号线程，控制CPU的线程为1号线程，分别运行在CPU编号为0和1的核上；并行处理网格的方法为：0号线程初始化指针N₀＝SPatch，将N₀→R置为True，执行第六步；1号线程初始化指针N₁＝SPatch→Next，将N₁→R置为True，转第九步；

第六步、0号线程控制GPU对N₀指向的网格进行处理，具体做法为：

6.1将N₀指向的网格的数据拷贝到GPU的内存中；

6.2启动GPU线程采用牛顿力学方程计算N₀中每个分子所受的力，根据力和初始速度计算每个分子当前速度和位置，将计算结果存放在GPU的内存中；

6.3将计算结果从GPU的内存拷贝到位于CPU的内存中的N₀→Ppos和N₀→Pvel数组中；

第七步判断N₀→Next是否为NULL，若是，转第十二步，否则执行第八步；

第八步、使N₀＝N₀→Next，判断N₀→R是否为True，若是转第七步，否则将N₀→R置为True并转第六步；

第九步、1号线程根据CPU的核数Q启动Q-2个子线程T₀，T₁，...，T_Q-3运行于CPU的编号为2，3，...，Q-1的核上，T_c采用牛顿力学方程计算单元

中分子所受的力，0≤c＜Q-3；T_Q-3采用牛顿力学方程计算单元

N₁→CIdx[(N₁→CN)-1]中分子所受的力；T₀，T₁，...，T_Q-3将计算所得的力存放到N₁→Pforce数组中，根据力和初始速度计算每个分子当前速度和位置，将计算结果存放到N₁→Ppos和N₁→Pvel数组中；用符号

代表上取整；用SPatch→CIdx[x]代表网格SPatch中序号为x的单元编号，0≤x＜SPatch→CN；

第十步、判断N₁→Next是否为NULL，若是，则转第十二步，否则执行第十一步；

第十一步、N₁＝N₁→Next，判断N₁→R是否为True，若是转第十步，否则置N₁→R为True并转第九步；

第十二步、等待0号线程和1号线程都运行到第十二步后，执行第十三步；

第十三步、使t＝t+1，判断t是否小于等于T，是则将并转第五步，否则转第十四步；

第十四步、将计算结果包括分子系统中所有分子的位置向量和速度向量写入位于硬盘的文件中，结束。

2.如权利要求1所述的一种基于CPU和GPU协作的分子动力学加速方法，其特征在于所述目标网格权重TW和目标模拟步数T的值由模拟需求确定，TW的值为期望网格包含的分子数目，T的值为需要模拟的步数，实际模拟步数t初始化为0。

Images