用于原子-离子动力学的GPU并行化MATLAB软件

《Computer Physics Communications》：GPU-parallelized MATLAB software for atom-ion dynamics

GPU加速的离子-原子经典力学模拟及MATLAB优化实现。摘要：本文开发了基于GPU并行化的MATLAB脚本atomiongpu.m，通过优化ode45算法实现ode45gpu，在CPU上速度提升7倍，GPU上百万级模拟加速达600-3500倍，支持多种原子/离子参数、陷阱势及初始条件，自动计算复合概率、散射角分布及轨迹可视化。

　　

Saajid Chowdhury | Jesús Pérez-Ríos

物理与天文系，石溪大学，石溪，11794，纽约，美国

摘要

我们提供了一个MATLAB脚本atomiongpu.m，该脚本可以利用GPU并行化技术每天运行数百万次独立的模拟实验，模拟一个被俘离子与低密度原子云的相互作用，计算被俘离子以及从远处开始的原子的经典轨迹。该脚本使用了ode45gpu，这是我们对MATLAB的ODE求解器ode45中使用的Runge-Kutta算法的优化和专用实现。首先，我们讨论了物理系统，并展示了如何使用ode45gpu在CPU上以大约7倍于ode45的速度解决问题；当在GPU上并行运行一百万条轨迹时，与在CPU上使用ode45相比，速度可以提高600到3500倍。接着，我们展示了如何轻松修改atomiongpu.m的输入参数，以适应不同类型的原子、离子、原子-离子相互作用、陷阱势能、模拟参数和计算硬件，从而使atomiongpu.m能够自动计算复合体形成的概率、可观测量（如散射角和复合体寿命）的分布以及特定轨迹的图谱。

程序概述

程序名称： atomiongpu.m

CPC库链接到程序文件： https://doi.org/10.17632/sjw4hzw9jx.1

开发者仓库链接： https://github.com/saajidchowdhury/supplementGPU

许可条款： CC0 1.0

编程语言： MATLAB R2023a，已安装并行计算工具箱

补充材料： https://github.com/saajidchowdhury/supplementGPU

问题性质： 模拟离子和原子的经典动力学（牛顿定律），使用多达数百万组不同的初始条件，存储最终状态和其他一些标量可观测量及其分布，并绘制特定轨迹。

解决方法： 实现了MATLAB的第四/五阶自适应时间步长Runge-Kutta算法ode45来求解常微分方程，我们编写了一个独立的函数ode45gpu。然后，我们使用MATLAB的arrayfun在多个CPU或GPU上对其进行并行化。最后，我们编写了包装脚本atomiongpu.m，以便快速方便地使用ode45gpu。

其他说明：atomiongpu.m、ode45gpu及其图表的源代码可以在仓库中找到。

引言

在许多领域，如化学物理、生物物理和材料科学中，分子动力学（MD）是一种常用的预测技术。通常，MD模拟系统的原子和分子的运动，遵循牛顿定律，持续一段时间。MD模拟的结果揭示了系统的性质，如热分布和反应速率[1]、[2]。

在原子和分子物理学中，MD的一个具体应用是研究被俘离子与自由原子的相互作用，这对于冷化学和超冷化学至关重要[3]。在原子-离子相互作用的实验研究中，Paul陷阱被用来在旋转的四极电场中捕获离子，而缓冲气体被用来冷却离子，以隔离原子-离子散射的量子力学方面[4]、[5]、[6]。此外，观察到被俘离子会暂时与原子弱结合，形成所谓的复合体。除了这些实验进展外，MD模拟还被广泛用于理解这些结果并进行理论预测，例如复合体形成的概率和热化速率[1]、[2]、[7]、[8]、[9]、[10]、[11]、[12]、[13]、[14]、[15]、[16]。

由于缓冲气体的密度通常较低，可以假设离子一次只与一个原子相互作用，因此只需要模拟两个粒子。然而，精确的初始条件是未知的，所以仅运行一次原子-被俘离子系统的MD模拟是不够的。通常需要运行数千次MD模拟，每次模拟都使用不同的初始条件，以获得全面的统计描述。然而，为了观察到高分辨率的分形动力学或计算出复合体形成的概率（精度达到0.1%），需要运行数百万次MD模拟。

因为每个原子-离子轨迹都是独立的，所以可以并行化这些模拟。此外，由于每个轨迹只包含两个粒子，并且通常不需要存储完整的轨迹，因此所需的内存非常少。因此，原子-被俘离子系统是使用GPU进行大规模并行化MD模拟的理想候选者。GPU包含数十到数百个微处理器，它们可以同时工作，每个微处理器都能进行低内存的基本计算。GPU加速的MD已经很常见[17]、[18]、[19]、[20]、[21]、[22]，但现有的软件并不适合这种冷化学问题。例如，一些MD软件包只能模拟被俘离子，但不包括与原子的相互作用[23]、[24]，而一些更通用的软件则针对非常大的分子进行了优化[25]，但没有实现陷阱势能或可观测量（如反弹次数）。此外，从未有过简单的开源MATLAB实现来并行化经典MD。

要运行MD模拟，必须数值求解常微分方程（ODE）。MATLAB有多种求解ODE的方法，包括基于Dormand-Prince对的第四/五阶自适应时间步长Runge-Kutta方法ode45[26]、[27]、[28]。要在GPU上运行计算，通常需要重新设计和实现计算算法，以适应CUDA接口[29]、[30]。幸运的是，一些MATLAB函数也可以在GPU上运行，利用MATLAB将基本操作内部转换为CUDA代码的能力。这使我们能够编写一个MATLAB函数ode45gpu，该函数既可以在CPU上运行，也可以在GPU上运行。

本文的结构如下。首先，我们详细描述了所模拟的物理系统，最终得到一个形式为ode45来解决这个系统，以及如何通过以最佳方式编写来显著提高性能。接下来，我们介绍了我们在MATLAB中实现的第四/五阶自适应时间步长Runge-Kutta算法ode45gpu，它在单个CPU核心以及多个CPU核心并行运行时都优于ode45。然后我们讨论了如何使用gpuArray和arrayfun在GPU上使用ode45gpu，包括我们如何编写。我们还讨论了ode45gpu在不同类型GPU上的性能。最后，我们描述了atomiongpu.m的输入和输出，这是一个MATLAB脚本，它可以自动化并简化ode45gpu的使用，以便研究原子-离子相互作用。

章节片段

物理系统

MATLAB脚本atomiongpu.m能够模拟由原子和离子组成的任何系统，这些离子可能是未被俘获的，或者处于谐振陷阱或Paul陷阱中。它忽略了内部自由度，并假设原子-离子相互作用的形式为$V\left(r\right)=C_n/r^n+C_m/r^m$。对于并行化，原子的初始速度可以是朝向离子的固定速度，也可以是在随机方向上的热分布；代码还可以模拟特定的自定义情况

在CPU上使用ode45

通常，要数值求解，可以使用ODE求解器，例如MATLAB的ode45。需要向求解器提供函数、所需的模拟时间t_max以及诸如绝对和相对容差之类的选项。然后，求解器以离散步骤dt进行积分，模拟一段时间，在我们的案例中是$t_{\text{max}}=1$毫秒。在每个时间步长，如果估计的误差超过由容差确定的某个阈值，则

Atomiongpu.m：ode45gpu的用户界面

在上一小节中，我们描述了如何从零开始在GPU上使用ode45gpu。然而，并不是每个人都有时间详细编辑代码，我们可以为被俘离子与低密度原子云的相互作用这一特定用例自动化代码。因此，为了使ode45gpu更易于使用，我们开发了一个名为atomiongpu.m的“用户界面”脚本。用户只需更改图10中显示的几个参数，其余工作由代码完成。

结论

在本文中，我们介绍了我们的MATLAB函数ode45gpu，它可以在CPU或GPU上求解描述原子-被俘离子动力学的常微分方程，在CPU上的运行速度比MATLAB的ode45快7倍，在GPU上并行运行时速度可以提高600到3500倍。然后，我们提供了一个名为atomiongpu.m的脚本，用户可以轻松自定义原子、离子、陷阱、ODE传播、初始条件和并行化的所有参数，让代码完成其余工作

CRediT作者贡献声明

Saajid Chowdhury：撰写 – 审阅与编辑，撰写 – 原始草稿，可视化，软件开发。Jesús Pérez-Ríos：撰写 – 审阅与编辑，监督，方法论，资金获取，概念化。

利益冲突声明

作者声明以下可能被视为潜在利益冲突的财务利益/个人关系：Saajid Chowdhury报告称获得了美国空军科学研究办公室的财务支持。如果有其他作者，他们声明没有已知的可能会影响本文所述工作的财务利益或个人关系。