基于FPGA的分子动力学模拟交互控制系统

王 鑫,吴俊辉

(1.江南大学教育部轻工过程先进控制重点实验室,江苏 无锡 214122; 2.江苏省未来网络创新研究院,江苏 南京 211111)

0 引 言

分子动力学(Molecular Dynamics， MD)模拟一直是高性能计算领域的热门话题，在化学、材料、生命科学等多个学科领域具有广泛的应用[1-5]。分子动力学模拟主要基于牛顿第二定律来实现对分子运动状态的追踪，这需要花费大量时间去计算粒子间的作用力，特别是短程力[6-8]。在进行大时间尺度的模拟时，需要在极短时间内对粒子的特征参数进行高次迭代来计算短程力。对于以CPU为主导的分子动力学模拟系统，这会产生高额的CPU设备数据传输成本[9]。

为了降低数据通信过程中的高额成本以及提升短程力计算效率，本文放弃传统的以CPU为主导的设计，采用基于FPGA加速卡的系统设计[10-13]。通过引入FPGA，提出一种交互控制系统。该系统可以使用片上存储缓存数据，减少数据通信频率[14]；同时可以进行并行计算[15-17]来提高计算效率。

FPGA加速卡的核心是短程力计算模块。在分子动力学模拟过程中，短程力计算模块不断地进行粒子间的短程力运算，而粒子信息的传输速度影响着短程力计算的效率[18]。短程力计算模块支持多组粒子间短程力的并行计算[19]。计算的粒子组数越多，相应的数据传输压力就越大。同时，并行计算多组粒子之间的短程力时，还可能出现同时读取同一粒子数据的情况，这又会引发访存的冲突问题。如何利用FPGA上有限的硬件资源来实现粒子信息从片上存储到短程力计算模块的快速可靠的传输就显得至关重要。针对该问题，本文提出一种交互控制系统(Interactive Control System, ICS)，用较少的FPGA硬件资源，实现粒子信息的高效传递，大大减轻了计算模块在并行化时所面临的数据传输压力，同时避免了访存的冲突问题。

1 总体设计

本文FPGA加速卡的总体架构设计如图1所示。上位机将二进制的粒子数据通过PCIE总线传递给FPGA加速卡上的片上存储BRAM(Block RAM)，接着发送开始信号到ICS。ICS收到开始信号后，取出BRAM中的数据，并根据粒子配置对数据进行解析，之后将完整的粒子信息发送给短程力计算模块进行运算。短程力计算模块在完成运算后，向ICS发送新一轮的数据请求，并将力的计算结果发送给力回写模块，令其将数据回写到BRAM。当完成所有数据的运算后，ICS反馈给上位机一个结束信号。上位机收到结束信号后，将BRAM中的力信息取出，用于粒子运动信息的更新。

图1 基于FPGA的分子动力学模拟系统总体设计框图

1.1 短程力计算的NT算法

本文的分子动力学模拟系统中，短程力计算模块的设计思路与波士顿大学提出的设计[12,20]类似，是基于NT算法(Neutral Territory method)[21-25]来实现粒子间配对和计算的。这是Shaw[25]在2005年开发的一种并行化短程力计算多体问题的方法。该方法中，一块物理空间被划分为多个tower区域(固定宽度的长条形区域)和plate区域(以tower区域的中间段为中心形成的半圆盘形区域)[25]。在进行分子动力学模拟时，需要同时计算各个tower区域中粒子与对应plate区域中粒子间的交互作用。

为了便于计算，短程力计算模块将tower区域与plate区域分别划分为固定大小的单元。tower区域单元中的粒子静态缓存在短程力计算模块内，plate区域单元中的粒子依次动态流入短程力计算模块，与缓存的tower区域的粒子进行短程力的计算。本文把上述2类单元形象地称为静态盒子和动态盒子，并统称为盒子。鉴于粒子在物理空间中随机分布，粒子在盒子中的分布也是随机的，且每个盒子内的粒子数目不定。

从总体设计与NT算法的原理可以看出，ICS工作逻辑的设计直接影响着短程力计算模块的工作效率，而BRAM内的数据结构[26]又影响着ICS的数据解析效率。

1.2 BRAM数据结构设计

本文设计的ICS硬件系统在逻辑上将BRAM按盒子进行划分。这些盒子有各自的编号，如图2所示。为了便于管理，每一块盒子所对应的存储器空间内，可存储的粒子数据最大容量保持一致，而需要使用的粒子数据则是存储在不同盒子对应的逻辑存储空间的不同位置(对应粒子在物理空间中随机分布的状态)，所以，BRAM内缓存的粒子数据就存在有效(需要使用的数据段)与无效(不关心的数据段)的区别。如果一个盒子中的所有粒子数据都是无效的，该盒子就被称为空盒子。粒子数据包括用于寻址的粒子盒子编号和盒内编号，以及用于计算短程力的粒子的三维坐标、电荷数、sigma(σ)和epsilon()[27]。

图2 片上BRAM的数据分布图

为了提高系统的工作效率，避免对无效数据段进行解析，本文引入位图来获得有效粒子在BRAM中的具体位置，进而完成数据的高效传递。具体来说，采用位图的每一位依次与粒子的位置空间对应，1表示该位置的数据是有效粒子数据，0表示无效的。这样就无需取出整个盒子的所有粒子数据，而是由位图得到每一个有效粒子的确切位置，同时也避免了解析全是无效数据的空盒子。这样的数据结构有助于提高粒子数据的传输速度和后续短程力计算模块的计算速度。为简化描述，后文提及的粒子在不加说明的情况下，指的都是有效粒子。

除了粒子数据，分子动力学模拟涉及的另一类重要数据是力数据。注意到，BRAM的访问端口至多2个，而其中一个端口需要用于加载二进制数据文件。为了提高粒子数据解析模块和力回写模块的并行化程度，避免粒子数据读取和力数据读写在访存时出现冲突，本文将力数据存储于另一块BRAM中[28]，划分方式与粒子数据一致，数据内容为粒子所受的力在3个方向的分量为f_x、f_y、f_z，如图2所示。

1.3 ICS工作逻辑设计

基于BRAM数据结构的设计方案，本文将ICS分为2个模块，分别是取数控制模块和粒子数据解析模块(见图1)。整个ICS的工作流程如下：ICS一开始处于待机状态，上位机先将二进制的位图数据和粒子数据通过PCIE发送到FPGA加速卡，并按照1.2节设计的数据结构存储在FPGA加速卡的BRAM中，同时将BRAM中对应的力存储空间的力清零；上位机完成数据写入后，再发送盒子对号(用于指定计算哪2个盒子内粒子间的短程力)给ICS的取数控制模块，取数控制模块根据盒子对号与图2所示的BRAM地址映射关系得到粒子数据的存储地址，并发送给粒子数据解析模块；粒子数据解析模块在收到地址后，通过BRAM_CTRL(自定义的与BRAM交互的控制块)从BRAM中依次读取出粒子数据，然后根据预定义的粒子配置(包含短程力计算模块所需静态粒子个数、粒子配置信息各参数及对应位宽等)对二进制数据进行解析，并输出给短程力计算模块；当完成一组盒子号的所有粒子的力计算后，短程力计算模块将结果(包含力信息和粒子位置信息)发送给力回写模块，然后，力回写模块更新BRAM中对应的力数据，同时，ICS开始下一轮的数据解析与传输；在完成所有盒子对中粒子的数据解析并收到短程力计算模块的运算完成信号后，取数控制模块将结束信号反馈给上位机，上位机再将力数据从BRAM中取回，用于粒子运动信息的更新。

2 ICS详细设计

2.1 取数控制模块

为了避免从BRAM中读取粒子数据时引发冲突问题，本文使用如下设计：从BRAM中依次取出相应盒子内的粒子数据，解析处理后按照先静态粒子后动态粒子的顺序传输给短程力计算模块。为了进一步提高传输效率，先取一个静态盒子内的粒子数据进行缓存，然后将所有需要计算力的动态盒子内的粒子数据依次取出，从而避免对同一个静态盒子内粒子数据进行多次解析。该设计在节省传输带宽的同时，又与短程力计算模块的计算速度相匹配，为实现分子动力学模拟系统的流水线设计提供支持。

取数控制模块的设计如图3所示。图中的FIFO_S与FIFO_D是2个先入先出的双口缓冲器。该模块工作流程如下：

图3 交互控制系统的取数控制模块的设计

1)在接收到上位机发来的盒子对组号时，分别将其存入FIFO_S和FIFO_D中。

2)在接收到上位机发来的开始指令后，前述的2个FIFO同时将顶端的盒子号分别输出到寄存器S和寄存器D中进行暂存。

3)将寄存器S中的盒子号与目前顶端的盒子号进行比较，判断下一轮的静态盒子是否产生变化。同时根据盒号与地址空间的映射关系求出2个寄存器中盒号所对应的位图地址与盒子地址，并将地址与比较器结果一起发送给粒子数据解析模块。

4)根据粒子数据解析模块反馈的信号，返回至步骤2开始FIFO中下一组盒号对的输出；或根据短程力计算模块的反馈信号结束整个取数过程。

2.2 粒子数据解析模块

粒子数据解析模块是ICS的核心部分，也是与片上BRAM交互最为频繁的部分。其内部通过精准的时延和逻辑控制，实现了对粒子数据高效快速的解析与输出。该模块的基本处理原则是，先对静态盒子进行操作，再对动态盒子进行操作。该模块的设计如图4所示。

图4 交互控制系统的粒子数据解析模块的设计

1)该模块将取数控制模块输出的比较器结果用寄存器暂存，并对输入进来的2组位图地址与盒子地址进行选择。地址选择器先对盒子地址进行暂存，用于后面计算该盒内的粒子数据的地址。根据位图地址，通过BRAM_CTRL与BRAM进行交互得到位图数据，并输出一个盒子标志位到寄存器T以记录当前解析的盒子类型。对位图数据进行分析后，计算盒子中粒子的个数并用寄存器Q暂存，并将粒子数据所处的盒内位置送入FIFO。根据FIFO中的索引得到粒子数据在盒内的偏移地址，加上之前暂存的盒子地址，就可以计算出粒子数据的具体地址。

2)根据粒子数据的具体地址，取出粒子数据进行解析。为了减少数据占用的空间大小，粒子的各个信息要尽可能连续存放，所以解析的时候需要根据软硬件约定的数据结构对数据进行位的拆分与重新赋值。

3)将已赋值的粒子信息输出到短程力计算模块，记录已解析完成的粒子个数，并与前面记录的单个盒子的粒子个数进行比较。若未完成一个盒子内的所有粒子的数据解析，则从FIFO中取出下一个索引，重复前述操作。

在实际的工作流程中，每个静态盒子内的粒子个数是不确定的，而短程力计算模块在更换缓存在模块内的静态粒子时，会出现输入减少、模块无法正常工作的问题。为了保证短程力计算模块每次都有等量的静态粒子输入，本文引入空粒子(所有参数都为0，不参与实际计算)来补充输入。当完成单个静态盒子内所有有效粒子的数据解析后，根据短程力计算模块的静态粒子配置(默认配置为单个盒子的最大粒子容量)要求，结合当前的解析的粒子类型来进行空粒子的补充输出。

在完成动态盒子中粒子数据的解析之后，根据FIFO_A(用于暂存输入到短程力计算模块的粒子)的存储空间来决定是否进行下一轮的数据解析。若下一轮的静态盒子保持不变，在进行新的一轮的数据解析时，则跳过静态盒子内的粒子数据解析，因为短程力计算模块已对其进行了暂存。

图5、图6所示的状态转移图对动静态盒子的整个数据解析过程做了更为细致的描述。

图5 粒子数据解析(静态盒子)的状态转移图

图6 粒子数据解析(动态盒子)的状态转移图

3 系统的测试与分析

本文使用Verilog硬件描述语言实现了ICS的设计，并在之后进行了波形仿真和板级验证。在测试过程中使用的粒子数据集(7051个有效粒子数据)是来自蛋白质晶体结构资料数据库(Protein Data Bank， PDB)的多聚泛素酶基因(UBQ)数据集。该数据集中的浮点数据已经进行了定点转换，统一采用了16进制的定点数表示。作为系统测试的准备工作，根据2.2节的数据结构设计，将所有的粒子数据随机划分到1372个盒子中，并生成对应的位图数据文件。

3.1 仿真测试

首先，在Vivado中对整个系统进行小数据集的仿真。将位图数据和粒子数据的二进制文件通过DMA加载到Vivado封装好的BRAM上。这里对(0，0)号盒子对所对应的粒子数据进行解析，粒子数据解析模块的结果如图7和图8所示。

图7最左边第1行的clk为工作时钟，最后一行的valid表示该粒子信息是否有效，倒数第2行至第4行为用于短程力计算模块的2个控制位，依次为is_real(空粒子标志)和is_static(动静态粒子标志)，其余行为粒子所带的属性。从图7中可以看到一共输出了28个粒子，其中有16个静态粒子(4个为空粒子)，12个动态粒子。作为代表，图8展示了11号粒子和空粒子的波形细节图。

图7 0号盒子的粒子数据解析结果

图8 典型粒子与空粒子的波形细节图

图9则是数据集用于短程力计算的0号盒子数据，从左到右依次为number(盒子编号)、idx(盒内编号)、loc(粒子的三维地址)、ec(电荷数)、用于计算短程力的2个参数sigma和epsilon。

图9 0号盒子的粒子信息(定点表示)

对比图8与图9，可以发现2幅图中关于11号粒子的信息是完全一致的。更多的实验结果显示，其余粒子的信息也是一致的。限于篇幅，这里不再详细展示。需要注意的是，粒子的信息并未实现连续的解析输出。这是因为，目前BRAM的数据宽度设定的是128位，而一个粒子的完整信息需要解析512位的数据才能得到，所以粒子信息是每4个时钟周期输出一个。

3.2 硬件测试

下面在硬件上对本文设计的系统进行测试。为了验证实用性与有效性，本文将设计好的ICS与短程力计算模块以及力回写模块进行联合测试。测试系统由Intel Xeon Silver 4216处理器以及Xilinx Virtex UltraScale+ HBM VCU128 FPGA组成。上位机首先通过PCIE将所有盒子的位图数据和粒子数据加载到片上BRAM中，然后发出41931组盒子对号(这些盒子对组合内的粒子间可能存在短程力)。ICS在收到盒子对号后开始工作，辅助整个系统完成粒子数据的解析、短程力的计算以及力信息到BRAM的回写。最后，上位机从BRAM中取出最终的结果，并与基于NAMD软件模拟的结果进行对比。为便于比较，已将数据转换浮点类型。

对比结果表明，短程力在3个方向上的数值在小数点后6位内完全一致，个别不一致的数据已在表1列出。可以看到，这些数据在小数点后5位前都是一致的，完全在允许的误差范围内。这说明本文设计的ICS可靠地实现了与分子动力学模拟系统各模块的交互及数据解析。

表1 ICS硬件测试的力数据对比表

关于ICS的硬件资源使用情况如表2所示。可以看到，该部分的资源占用百分比很低，这为其它模块的设计开发留存了丰富的资源。

表2 片上资源消耗

4 结束语

本文设计的用于分子动力学模拟的交互控制系统以经典的NT算法为基础，实现了粒子数据的有序传输。它具有以下几个特点：

1)考虑了粒子在空间中较为分散的特点。

2)加入位图数据来减少对无效数据的解析。

3)采用依次解析盒子内粒子数据的方式来减轻带宽压力。

4)通过监测静态盒子的变化避免静态盒子内粒子数据的多次重复解析。

5)可以根据短程力计算模块的计算速度来自动调节粒子数据的解析速度。

该系统很好地协调了各个模块的工作，而且能实现向多个计算模块同时供数，有效提高了计算效率。然而，它也有一些需要完善的地方，比如在取数控制模块，使用2个FIFO来存储盒子对号消耗的资源较多。在未来的工作中，拟根据短程力计算模块的需求来对该系统进一步优化，以实现更高的性能。