oneMD：现代 C++ 和 SYCL 中的分子动力学模拟器
5.00/5 (2投票s)
oneMD 是一款使用 C++ 17 编写的分子动力学模拟器，旨在利用 Intel 的 oneAPI 框架和 SYCL（并行处理的 OpenCL 标准的后继者）在多种设备和架构上利用硬件数据并行性。
引言
在原子和分子尺度上模拟物质是 HPC（高性能计算）和超级计算最迷人的应用之一，并且是许多不同科学研究领域（从材料科学到生物化学，再到药物发现）发现过程的关键部分。牛顿力学理论的发现和发展早在 18 世纪就激发了人们对 N 体牛顿系统时间演化的兴趣，并且 MD 中使用的几种数值方法和技术比使用数字计算机进行高速数值计算和仿真的开发还要早。使用适合电子计算机实现的数值方法模拟数百个相互作用的经典粒子理论的开发始于 20 世纪 50 年代末，而在此之前是蒙特卡罗模拟，蒙特卡罗模拟在二战后不久开发，用于研究核武器，并且是在 ENIAC I 计算机上作为存储程序运行的第一组计算。
oneMD 是一个项目，最初是为了参加 CodeProject 上的“跨架构挑战赛”而创建的。它是一款分子动力学模拟器，其原理与OpenMM 和GROMACS 类似，但使用 C++17 编写，并利用 Intel 的 oneAPI DPC++/C++ 编译器，该编译器实现了 Khronos Group 的SYCL 标准。oneAPI 允许您通过一套 C++ 抽象来利用 CPU 的 SIMD 功能以及 GPU 和 FPGA 加速，以实现加速的并行数据处理。
通常，使用 OpenCL 或类似的并行计算框架（如 C++ 中的 CUDA）意味着您需要在具有 .cl 或 .cu 扩展名的单独源文件中编写内核。（例如，OpenMM 的 OpenCL 平台），使用 C 或 C++ 的变体，如OpenCL C 或CUDA C++。每种并行计算框架都需要自己的语言和概念，并且每种硬件加速设备都有其自身的架构，程序员需要学习才能有效编程。
借助 DPC++ 和 SYCL，您可以使用单一源文件和纯 C++ 17 代码，并将 SYCL 内核封装为常规 C++ 函数或 C++ lambda 表达式。您可以使用相同的 C++ 类型和抽象集来编写通用的数据并行内核，这些内核可以在任何硬件加速设备上运行。您还可以编写针对不同设备的专用内核，包括具有 SIMD 向量单元的通用 CPU、GPGPU 以及具有可编程逻辑单元的 FPGA。
虽然单源模型大大降低了为不同设备编写内核的复杂性，但用 C++ 编写 DPC++ 内核与通用编程有很大不同。主机代码和设备代码之间存在严格的分离，并且设备代码使用受限的设备内存模型，其中所有设备内存缓冲区或分配都必须在使用前声明，并且您必须通过 SYCL 缓冲区显式捕获或共享来自常规 C++ 容器（如 std:vector 和 std:array）的任何数据。
考虑到这些差异，以及我希望从零开始使用现代 C++ 和 CMake 构建项目的愿望，因此，尝试使用 C++17 编写新的模拟器代码库，以 DPC++ 编译器工具链为目标，并在 Intel 免费的DevCloud 云计算环境中进行测试，该环境可访问多种硬件加速设备，这是有道理的。oneMD 也可以使用传统的并行计算库（如 OpenMP）为非 SYCL C++ 编译器（如 GCC 和 MSVC）进行构建，这使我们能够获得数据并行操作和内核的性能基准。
背景
分子动力学模拟器是一种计算机程序，它使用经典力学方程足以准确模拟粒子相互作用的假设，来模拟一组原子或分子在特定区域内如何相互作用。
在分子动力学模拟中，我们有一个由 N 个原子（或分子）组成的系统，其中每个原子都遵循牛顿定律，该定律将作用在单个原子上的力与其质量和加速度联系起来，即 F=ma。该模拟在微观层面计算每个原子的动量和位置，然后使用统计力学将这些信息转换为宏观属性，如压力、能量等。系统的热力学或宏观状态是通过对由位置和动量决定的单个原子的微观状态进行统计平均来计算的。
对于一个立方体区域内 N 个粒子的系统，oneMD 模拟的基本算法是：
-
通过使用随机分布为每个原子分配速度和位置来建立一组初始条件。
-
对于每个原子，跟踪其邻近原子所在的邻居列表。最初距离太远而无法相互作用（根据模型）的原子应从邻居列表中移除。
-
迭代一个固定的时间步数，并在每个时间步执行以下计算：
-
对于每个原子，迭代邻居列表，计算所选原子与其每个邻居之间的距离，并使用类似Lennard-Jones 的方程计算相互作用势能。原子由于邻近原子而受到的力是通过使用 LJ 方程关于距离的微分来找到的。
-
使用 virial 定理，将整个原子系统瞬时动能表示为每个邻近粒子产生的力与它们之间距离的函数。
-
使用动能、体积、温度和玻尔兹曼常数计算瞬时压力。
-
使用Verlet 积分方法计算下一个时间步中每个原子的下一个位置和速度。
开发环境
我们的开发环境将是 Windows 10 和 Visual Studio Code，配备C++ 扩展、CMake 扩展以及Remote SSH 扩展，连接到远程 DevCloud Ubuntu 环境。这种设置使我们能够拥有 C++ IDE 的所有便利，同时在 DevCloud 上构建和测试代码。
SSH
注册 DevCloud 后，您就可以访问关于设置 SSH 访问远程环境的说明。官方 Intel 说明涵盖了为 Windows 上的 Linux、macOS 或 Cygwin 环境设置 SSH 访问，但也可以使用 Windows 10 原生 SSH 客户端。首先，将密钥文件从 Intel 说明页面复制到您 Windows 帐户的 .ssh 文件夹中，路径为 %USERPROFILE%\.ssh。然后，将以下部分添加到 %USERPROFILE%\.ssh\config 中：
################################################################################################
# oneAPI DevCloud SSH config
################################################################################################
Host devcloud
User <user>
IdentityFile <path to key file>
ProxyCommand C:\Windows\System32\OpenSSH\ssh.exe -T -i <path to key file> guest@devcloud.intel.com
################################################################################################
# SSH Tunnel config
################################################################################################
Host *.aidevcloud
User <user>
IdentityFile <path to key file>
ProxyCommand C:\Windows\System32\OpenSSH\ssh.exe -T devcloud nc %h %p
LocalForward 5022 localhost:22
LocalForward 5901 localhost:5901
################################################################################################
################################################################################################
# DevCloud VSCode config
################################################################################################
Host devcloud-vscode
UserKnownHostsFile /dev/null
StrictHostKeyChecking no
Hostname localhost
User <user>
Port 5022
IdentityFile <path to key file>
################################################################################################
请将 <user> 和 <path to key file> 替换为您 DevCloud 用户 ID（uxxx...）和 DevCloud 密钥文件的路径。请注意，在上面的配置片段中，与 Intel 针对 Unix 环境的官方说明相比，有一个不同之处，即我们在 ProxyCommand 配置字段中使用了 Windows SSH 可执行文件的完整路径。
DevCloud
一旦您在 DevCloud 上获得了 shell，您就需要使用 qsub 和 qstat 等命令来设置用于开发或在 Intel 硬件（如 Xeon CPU 和 Iris Xe GPU）上运行程序的环境。这些队列命令继承自老式的集群管理软件，作业被提交到作业队列进行批处理，因此很容易找到命令的文档和示例。Intel 官方的关于提交和管理 DevCloud 作业的参考文档在此处：here。以下是我用于在 DevCloud 上开发 DPC++ 项目的一些基本命令。
为了开发，我们需要创建一个交互式作业，该作业将在我们登录 devcloud 的时间内分配内存和 CPU 资源。要创建运行 24 小时的交互式作业，请使用 qsub 命令，并带有 -I（大写 I）选项指定其为交互式选项，并带有 -l（小写 l）选项设置壁钟时间等参数。
uxxxxx@login-2:~$ qsub -I -l walltime=1:0:0:0
qsub: waiting for job 809147.v-qsvr-1.aidevcloud to start
qsub: job 809147.v-qsvr-1.aidevcloud ready
########################################################################
# Date: Tue Mar 9 11:26:21 PST 2021
# Job ID: 809147.v-qsvr-1.aidevcloud
# User: u61437
# Resources: neednodes=1:batch:ppn=2,nodes=1:batch:ppn=2,walltime=24:00:00
########################################################################
uxxxxx@s001-n045:~$
要查看正在运行的作业，请使用 qstat 命令。
uxxxxx@s001-n045:~/oneMD$ qstat
Job ID Name User Time Use S Queue
------------------------- ---------------- --------------- -------- - -----
809147.v-qsvr-1 STDIN u61437 00:00:04 R batch
要终止正在运行的作业，请使用 qdel 命令，并提供 qstat 显示的作业 ID。
uxxxxx@s001-n045:~$ qstat
Job ID Name User Time Use S Queue
------------------------- ---------------- --------------- -------- - -----
809147.v-qsvr-1 STDIN u61437 00:00:04 R batch
uxxxxx@s001-n045:~$ qdel 809147
Terminated
uxxxx@s001-n045:~$
qsub: job 809147.v-qsvr-1.aidevcloud completed
uxxxxx@login-2:~$
Visual Studio Code
一旦我们能够在 DevCloud 上创建交互式作业，我们就可以使用 Visual Studio Code 的 Remote SSH 扩展连接到交互式环境。请注意，我们需要两次登录 - 一次登录到基础 DevCloud 环境以启动作业，另一次登录到交互式作业环境的主机名（例如，s001-n045.aidevcloud），VS Code 将使用该主机名。我们的 SSH 配置文件已包含上述配置，因此我们只需输入类似 ssh s001-n045.aidevcloud 的命令，然后启动 VS Code 即可。
项目仓库中的 .vscode 文件夹包含用于使 IntelliSense 工作以及在 VS Code 中构建和调试 oneMD 的文件，例如 c_cpp_properties.json 文件包含项目和 DPC++ 库的包含文件夹路径。
有了这个环境，我们就可以像使用本地安装的工具链一样轻松地在 DevCloud 上进行 C++ 开发。
CMake
DevCloud 上当前版本的 CMake 是 3.10。这是一个较旧的版本，您在尝试构建其他库和应用程序时可能会遇到问题，例如，最新支持 SYCL 的 GROMACS 模拟器版本需要 CMake 3.20。DevCloud 不允许使用 sudo，任何您需要的其他程序都必须安装在您的主目录中。
幸运的是，很容易获取一个私有版本的 CMake 并将 VS Code 配置为使用它。您可以从 Kitware 下载最新版本的 CMake，作为 Linux 的二进制安装程序，例如 cmake-3.19.5-Linux-x86_64.sh。然后执行该文件，它将在您方便的位置（例如您的主文件夹）安装 CMake。在 VS Code 中，您可以更改 settings.json 文件中 CMake 可执行文件的路径。请注意，oneMD不需要这样做，因为它与已安装的 CMake 配合良好。
编写内核
我们可以简要看一下 DPC++ 内核的编写方式，与使用OpenMP 相比。oneMD 中 Lennard-Jones 模拟器的第一个迭代版本基于 James Wes Barnett 的 OpenMP LJ 模拟器。模拟器中的第一个内核是最简单的——更新每个原子的邻居列表。对于模拟中的每个原子，我们将其位置存储在一个名为 x 的 std::vector 中，作为 3D 向量。我们希望确定某个半径（称为 rlist）内的邻近原子，并将邻近原子存储在列表中。此操作可以并行化，因为每个原子的距离计算不会改变原子位置，并且它们彼此独立。最初的实现是这样的：
void NeighborList::Update(vector <coordinates> &x, cubicbox &box)
{
#pragma omp parallel for schedule(guided, CHUNKSIZE)
for (unsigned int i = 0; i < this->list.size(); i++)
{
this->list.at(i).resize(0);
}
// Atoms are not double counted in the neighbor list. That is, when atom j
// is on atom i's list, the opposite is not true.
#pragma omp parallel for schedule(guided, CHUNKSIZE)
for (unsigned int i = 0; i < x.size()-1; i++)
{
for (unsigned int j = i+1; j < x.size(); j++)
{
if (distance2(x.at(i), x.at(j), box) < rlist2)
{
this->list.at(i).push_back(j);
}
}
}
return;
}
在将每个原子的邻居列表向量初始化为零大小后，我们迭代每个位置向量，然后再次迭代每个后续向量，计算它们之间的距离，如果距离小于半径阈值，则将第二个向量添加到当前原子的邻居列表中。每个操作都 precedes by OpenMP #pragma 指令，用于使用 CPU 上的线程并行化循环。请注意，此内核最初是使用不同大小的 std::vector 来实现每个原子的邻居列表的。
与原始实现相比，oneMD 版本可以在不同的硬件设备上并行执行。oneMD 中的内核实现如下：
void System::UpdateNeighborList() {
dpc_common::TimeInterval timer;
auto n = static_cast<size_t>(natoms);
auto cut = this->conf.nlist;
int* neighbors = sycl::malloc_device<int>(n * n, q);
sycl::buffer<int, 2> n_buf (neighbors, sycl::range(n, n));
sycl::buffer<Vec3, 1> x_host_buf(&x[0], sycl::range(n));
auto k1 = q.submit([&](sycl::handler &h) {
#include "kernel_logger.hpp"
__kernel_logger("update_neighbor")
auto n_a = n_buf.get_access<sycl::access::mode::write>(h);
auto x_host_a = x_host_buf.get_access<sycl::access::mode::read>(h);
auto b = this->box;
h.parallel_for(sycl::range(n, n), [=](sycl::id<2> idx) {
auto i = idx[0];
auto j = idx[1];
if (i < j)
{
auto d = distance2(x_host_a[i], x_host_a[j], b);
if (d < cut)
{
printijdd2(idx, "Including this atom at distance ", d);
n_a[i][j] = 1;
}
else
{
n_a[i][j] = 0;
}
}
});
});
}
与 OpenMP 内核相比，我们使用固定大小的数组作为邻居列表，这些数组的大小相同，从而创建了一个 2D 网格，使我们能够更有效地并行化使用邻居列表的操作。邻居列表的内存直接在设备上分配。我们创建 SYCL 缓冲区，允许我们使用 SYCL 访问器在内核定义中从位置向量列表 x 读取并写入邻居列表数组。我们使用 C++ lambda 来定义内核，并以 2D 网格的形式并行迭代位置向量 x（作为 x_host_a 访问）。我们计算每对原子之间的距离，如果原子 x[j] 在阈值半径内，则在原子 x[i] 的邻居列表中存储 1。与第一个内核一样，我们不重复计算邻居，并且仅通过考虑当 i > j 时到原子 x[j] 的距离来存储原子 x[i] 的邻居。