1. 论文基本信息

1.1. 标题

NeuDATool: 一款开源中子数据分析工具，支持图形处理单元硬件加速和跨计算机集群节点并行 (NeuDATool: An Open Source Neutron Data Analysis Tools, Supporting GPU Hardware Acceleration, and Across-computer Cluster Nodes Parallel)

1.2. 作者

• Changli Ma
• He Cheng
• Taisen Zuo
• Guisheng Jiao
• Zehua Han

所属机构：

• 中国科学院高能物理研究所 (Institute of High Energy Physics, Chinese Academy of Sciences (CAS))，北京 100049，中国
• 散裂中子源科学中心 (Spallation Neutron Source Science Center)，东莞 523803，中国
• 中国科学院大学 (University of Chinese Academy of Sciences)，北京 100049，中国

1.3. 发表期刊/会议

预印本平台 arXiv。在相关领域，arXiv 是一个广泛接受且重要的预印本服务器，允许研究人员在同行评审之前分享其研究成果，以促进快速传播和讨论。

1.4. 发表年份

2019年4月12日 (UTC)

1.5. 摘要

经验势结构精修 (Empirical potential structure refinement, EPSR) 是一种中子散射数据分析算法和软件包，由英国散裂中子源 (ISIS) 的无序材料研究小组在1980年代开发，旨在构建无序液体中最可能的原子结构。在过去的几十年中，EPSR 已被广泛使用并产生了可靠的结果。然而，它采用 Fortran 语言编写，并实现了基于 OpenMP 的共享内存架构。随着超级计算机集群的广泛建设和图形处理单元 (Graphics processing unit, GPU) 加速技术的普及，现在有必要利用这些技术重建 EPSR，以提高其计算速度。本研究提出一个名为 NeuDATool 的开源框架。它采用面向对象语言 C++ 编写，可以在计算机集群内部的节点之间并行运行，并支持 GPU 加速。NeuDATool 的性能已通过水和无定形二氧化硅中子散射数据进行了测试。测试结果表明，该软件能够重建样品正确的微观结构，并且在包含约10万个原子的模拟盒中，与 CPU 串行算法相比，GPU 加速的计算速度可以提高400倍以上。NeuDATool 为熟悉 C++ 编程并希望为其分析定义特定模型和算法的科学家提供了另一种选择。

1.6. 原文链接

https://arxiv.org/abs/1904.08756

1.7. PDF 链接

https://arxiv.org/pdf/1904.08756v3.pdf

2. 整体概括

2.1. 研究背景与动机

2.1.1. 核心问题与挑战

中子全散射 (Neutron total scattering) 是研究无序材料（如液体、非晶态固体）原子结构的有力工具。经验势结构精修 (EPSR) 算法是这类研究中广泛使用的方法，它通过蒙特卡洛 (Monte-Carlo, MC) 模拟结合实验中子散射数据来构建最可能的原子结构。然而，EPSR 存在以下局限性：

• 编程语言过时且结构限制： 原始 EPSR 采用 Fortran 语言编写，这是一种过程导向语言，相比现代的面向对象语言，其代码的灵活性、可维护性和扩展性较差。
• 并行计算能力不足： 原始 EPSR 使用 OpenMP 实现共享内存架构并行，这意味着它只能在单个计算机节点内部的多个处理器核心上运行。它无法利用分布式内存架构，从而不能在跨多个节点的超级计算机集群上进行大规模并行计算。这严重限制了其处理更大系统（如包含数十万甚至数百万个原子的宏观分子系统）的能力。
• 未利用现代加速技术： 随着图形处理单元 (GPU) 等硬件加速技术的普及，科学计算的速度得到了显著提升，但原始 EPSR 未能利用这些技术。

2.1.2. 研究重要性

克服上述挑战对于推动无序材料结构研究至关重要。能够处理更大、更复杂的系统（例如具有数十万甚至数百万个原子的宏观分子或纳米粒子系统）将使科学家能够解决目前无法触及的科学问题，例如更精确地理解水在不同条件下的结构、混合溶液中的异质性，以及其他复杂液体的微观结构。提高计算速度也将缩短研究周期，加速科学发现。

2.1.3. 创新思路

为了解决原始 EPSR 的局限性，本文提出了 NeuDATool 框架，其创新思路在于：

• 现代化编程语言和面向对象设计： 采用 C++ 语言重新实现 EPSR，利用其面向对象特性，提高代码的灵活性、可扩展性和用户友好性，使得用户可以更容易地定义新的分子、原子类型和运动模型。
• 引入分布式并行计算： 利用消息传递接口 (Message Passing Interface, MPI) 技术，使程序能够在计算机集群的多个节点上进行并行计算，从而突破共享内存架构的限制，支持更大规模的模拟。
• 整合硬件加速技术： 引入计算统一设备架构 (Compute Unified Device Architecture, CUDA) 和 GPU 硬件加速，大幅提升计算密集型任务（如原子对数和径向分布函数的计算）的速度。

2.2. 核心贡献/主要发现

• 提出 NeuDATool 开源框架： 该框架是原始 EPSR 算法的现代化重构，采用 C++ 编写，支持面向对象编程。
• 实现跨节点并行和 GPU 加速： 整合了 MPI 用于跨节点并行计算，以及 CUDA/GPU 用于硬件加速，显著提高了计算效率和可扩展性。
• 验证算法的正确性： 通过水和无定形二氧化硅的中子散射数据测试，NeuDATool 成功重建了样品的正确微观结构，并与实验数据和现有研究结果一致。
• 展示显著的计算速度提升： 在包含约10万个原子的模拟盒中，GPU 加速的计算速度比 CPU 串行算法提高了400倍以上，使得模拟百万级原子系统成为可能。
• 提供灵活和用户友好的工具： C++ 的面向对象特性使得科学家能够方便地定义和修改特定的模型和算法，为中子散射数据分析提供了新的选择和更大的灵活性。

3. 预备知识与相关工作

3.1. 基础概念

3.1.1. 中子全散射 (Neutron Total Scattering)

中子全散射是一种利用中子束照射样品，通过探测散射中子的动量和能量变化来研究材料原子结构的技术。它特别适用于研究无序材料（如液体、非晶态固体、玻璃等），因为中子与原子核相互作用，对轻元素（如氢、氘）和同位素（如氢与氘）的散射能力差异显著，这使得它能够提供关于原子对关联和整体微观结构的详细信息。

3.1.2. 经验势结构精修 (Empirical Potential Structure Refinement, EPSR)

EPSR 是一种基于蒙特卡洛模拟的算法，用于从实验中子散射数据中推断无序材料的原子结构。其核心思想是通过调整原子间相互作用势（包括一个参考势和一个经验势），使模拟产生的结构因子 (Structure Factor) 与实验测量的中子散射数据（结构因子）尽可能匹配。通过这种迭代优化过程，EPSR 能够构建出与实验数据一致的最可能的原子配置。

3.1.3. 蒙特卡洛模拟 (Monte-Carlo, MC)

蒙特卡洛模拟是一大类使用随机数来解决问题的方法，尤其适用于物理、化学中的多体系统模拟。在原子分子模拟中，MC 模拟通过随机地对原子或分子进行位移、旋转等操作，并根据能量变化（通常使用 Metropolis 准则）来决定是否接受这些操作，从而探索系统的构象空间，最终达到统计平衡。

3.1.4. 径向分布函数 (Pair Correlation Function, PDF, $g_{\alpha\beta}(r)$)

径向分布函数 $g_{\alpha\beta}(r)$ 描述了在距离 $r$ 处找到一个 $\beta$ 类型原子，给定原点处有一个 $\alpha$ 类型原子的概率密度。它是描述无序材料微观结构的关键函数，能够揭示原子间的平均距离、配位数和短程有序结构。

3.1.5. 结构因子 (Scattering Pattern / Structure Factor, S(Q))

结构因子 S(Q) 是中子散射实验的直接测量结果，它反映了材料在倒空间 (Q 空间) 中的结构信息。$Q$ 是散射矢量，其大小 $Q = \frac{4\pi}{\lambda}\sin\theta$ 与中子波长 $\lambda$ 和散射角 $2\theta$ 有关。S(Q) 与径向分布函数 $g_{\alpha\beta}(r)$ 之间通过傅里叶变换 (Fourier Transform) 相关联。

3.1.6. 氘化技术 (Deuteration Techniques)

氘化技术是指通过化学方法用氘（氢的同位素）替换样品中的氢原子。由于氢和氘具有非常不同的中子散射长度，对样品进行不同程度的氘化可以显著改变其散射模式 S(Q)，而对原子结构本身影响很小（因为它们的化学性质非常相似）。这使得通过分析多个氘化样品的散射数据，可以解耦出不同原子对之间的径向分布函数，从而更全面地解析复杂系统的结构。

3.1.7. 图形处理单元 (Graphics Processing Unit, GPU) 加速

GPU 是一种专门为并行处理大量数据而设计的处理器。与传统的中央处理器 (CPU) 相比，GPU 拥有数千个小的处理核心，能够同时执行大量简单的计算任务。在科学计算中，GPU 加速通常用于计算密集型且高度并行的任务，例如本文中涉及的原子对距离计算。

3.1.8. 计算统一设备架构 (Compute Unified Device Architecture, CUDA)

CUDA 是 Nvidia 公司推出的一种并行计算平台和编程模型，允许开发者使用 C、C++ 等高级语言直接在 Nvidia GPU 上进行编程。它提供了一套完整的开发工具，包括编译器、库和调试器，使得 GPU 加速的应用程序开发变得更加容易。

3.1.9. 消息传递接口 (Message Passing Interface, MPI)

MPI 是一种用于分布式内存并行计算的标准和库。它允许在不同处理器（通常位于不同计算机节点上）之间通过发送和接收消息进行通信和数据交换。MPI 是构建大规模并行应用程序的基石，尤其适用于超级计算机集群。

3.1.10. Open Multi-Processing (OpenMP)

OpenMP 是一种用于共享内存多处理器编程的应用程序接口 (API)。它允许程序员通过编译器指令（#pragma）在 C、C++ 和 Fortran 代码中指定并行区域，从而在单个计算机节点内部的多个 CPU 核心上实现并行计算。

3.1.11. C++ (面向对象语言)

C++ 是一种通用的、静态类型的、支持多编程范式（包括面向对象、泛型和过程式编程）的编程语言。面向对象编程 (Object-Oriented Programming, OOP) 允许通过类 (class) 和对象 (object) 来组织代码，具有封装 (encapsulation)、继承 (inheritance) 和多态 (polymorphism) 等特性，提高了代码的模块化、可重用性和可扩展性。

3.2. 前人工作

3.2.1. 原始 EPSR 软件

EPSR 算法最初由英国 ISIS 散裂中子源的无序材料研究小组在1980年代开发，并以 Fortran 语言实现。该软件包在过去几十年中被广泛应用于无序液体结构的分析，并产生了大量可靠的科学结果。然而，其 Fortran 语言和基于 OpenMP 的共享内存架构限制了其在现代超级计算机集群上的性能和可扩展性。

3.2.2. 反向蒙特卡洛 (Reverse Monte Carlo, RMC)

RMC 是另一种与 EPSR 类似的，用于从实验数据（如中子或X射线散射）中重建原子结构的全局优化方法。与 EPSR 不同的是，RMC 通常不使用显式势能函数，而是直接通过随机移动原子并接受那些能改善模拟数据与实验数据匹配程度的移动来优化结构。EPSR 结合了 RMC 的思想和物理势能函数的概念，使其在物理合理性上通常更强。

3.2.3. 分子动力学 (Molecular Dynamics, MD) 力场

分子动力学 (MD) 模拟是另一种重要的原子模拟方法，它通过牛顿运动方程模拟原子和分子的运动。MD 模拟需要精确的原子间相互作用势，这些势被称为力场 (Force Fields)。本文提到的 reference potential (RP) 与 MD 模拟中使用的力场类似，其参数可以从现有的 MD 力场中获取，例如：

• OPLS (Optimized Potential for Liquid Simulations): 一系列针对液体模拟优化的力场。
• AMBER (Assisted Model Building with Energy Refinement): 一种广泛用于生物大分子（如蛋白质和核酸）模拟的力场。
• CHARMM (Chemistry at HARvard Molecular Mechanics): 另一种用于生物分子和化学系统模拟的通用力场。 这些力场为原子提供了合理的形状和相互作用，构成了 EPSR 算法中 reference potential 的基础。

3.3. 技术演进

中子散射数据分析工具的技术演进反映了计算科学和硬件技术的发展：

• 早期阶段 (1970-1980s): 随着第一代中子全散射谱仪的出现，对数据分析方法的需求增长。EPSR 在这一时期诞生，以 Fortran 语言实现，并使用 OpenMP 进行单节点共享内存并行。这在当时是先进的，但受限于当时的计算资源。
• 中期阶段 (1990s-2000s): 超级计算机集群开始普及，分布式计算（如 MPI）成为处理大规模科学问题的标准。然而，许多遗留代码（如原始 EPSR）未能及时更新以利用这些新的并行范式。
• 现代阶段 (2000s至今): GPU 加速技术异军突起，为计算密集型任务提供了前所未有的加速能力。同时，面向对象语言（如 C++）在软件工程中的优势日益凸显，使得代码更易于开发、维护和扩展。NeuDATool 正是这一技术演进的产物，旨在将 EPSR 算法现代化，使其能够充分利用最新的计算硬件和软件范式。

3.4. 差异化分析

本文提出的 NeuDATool 与原始 EPSR 及其相关工作的主要差异和创新点体现在以下几个方面：

NeuDATool 的核心创新在于将现代计算技术（面向对象编程、分布式并行、GPU 加速）引入到经典的 EPSR 算法中，从而解决了原始 EPSR 在处理复杂系统和追求更高计算效率方面的瓶颈。

4. 方法论

NeuDATool 本质上是一种蒙特卡洛 (Monte-Carlo) 模拟方法，但与传统的 Metropolis MC 模拟不同，它将中子散射实验数据信息融入到原子相互作用势中，以指导模拟过程。

4.1. 方法原理

在 NeuDATool 中，蒙特卡洛模拟使用的原子间相互作用势分为两类：参考势 (reference potential, RP) 和 经验势 (empirical potential, EMP)。

• 参考势 (RP): 类似于分子动力学 (Molecular Dynamics, MD) 模拟中使用的力场。它用于赋予分子合理的形状，并实现其他已被验证的物理化学约束。RP 的形式和参数可以从现有的全原子 MD 力场中获取，如 OPLS、AMBER 和 CHARMM 等。
• 经验势 (EMP): 没有固定的形式，其作用是反映模拟产生的结构因子 $S_{sim}(Q)$ 与实验测量的结构因子 $S_{exp}(Q)$ 之间的差异 $\Delta S(Q)$。具体来说，EMP 是 $\Delta S(Q)$ 的逆傅里叶变换，作为对 RP 的扰动，引导模拟结果趋近于实验测量值。

4.1.1. 经验势的表达形式

在 NeuDATool 的当前实现中，经验势 (EMP) 是基于一系列在实空间和 Q 空间中的泊松分布 (Poisson distributions) 及其相应的傅里叶变换来表示的。

该图像是一个示意图，展示了在R空间（上）和Q空间（下）中不同au值下的泊松分布。上半部分的分布用于拟合riangle S(Q)，图中的峰值经归一化处理，所有分布的峰值均为1.0。

图 1（原文 Figure 1）展示了一系列在R空间（上图）和Q空间（下图）中具有不同 $\lambda$ 值（数学期望）的泊松分布。这些Q空间中的分布列表被用于拟合 $\Delta S(Q)$。每个分布的峰值都被归一化到1.0。

4.1.2. 模拟过程中的势能更新

模拟开始时，MC 模拟只使用 参考势 (RP)。系统势能的变化 ($\Delta U = U_{after} - U_{before}$) 用于判断分子或原子的随机移动是否被接受。

• 如果 $\Delta U < 0$，移动被接受。

• 如果 $\Delta U > 0$，移动以 $e^{-\Delta U/kT}$ 的概率被接受，其中 $k$ 是玻尔兹曼常数， $T$ 是温度。 当 MC 模拟达到平衡后，程序计算 经验势 (EMP)。EMP 通过拟合模拟和实验之间的结构因子差异 $\Delta S(Q)$ 得到。然后，这个计算得到的 EMP 会被添加到 RP 中，形成更新后的总势能，继续进行模拟。这个过程会重复进行，直到 EMP 变得非常小（即 $\Delta S(Q)$ 趋近于零），表示模拟结果已与实验数据充分吻合。

  该图像是图表，展示了不同_ ext{H}_2 ext{O}样品的riangle S(Q)拟合结果（左图）及对应的经验势（右图）。其中，HD ext{H}_2 ext{O}、FH ext{H}_2 ext{O}和FD ext{H}_2 ext{O}分别表示半氘化、完全氢化和完全氘化的水样品，riangle S(Q)的值及势的变化随$Q$和$R$的变化而变化。

图 2（原文 Figure 2）展示了 $\mathrm{H_2O}$ 样品（HD $\mathrm{H_2O}$ 表示半氘化 HDO，FH $\mathrm{H_2O}$ 是完全氢化的 $\mathrm{H_2O}$，FD $\mathrm{H_2O}$ 是完全氘化的 $\mathrm{D_2O}$）的 $\Delta S(Q)$ 拟合结果（左图）以及相应的经验势（右图）。左图显示了使用泊松分布拟合 $\Delta S(Q)$ 的效果，右图则展示了 O-H 和 O-O 原子对的经验势。

4.2. NeuDATool 算法流程

NeuDATool 的算法流程如下图（原文 Figure 3）所示，以下将详细描述各个步骤。

该图像是示意图，展示了NeuDATool的算法流程。图中描绘了从定义原子、分子和模拟盒类开始，通过重编译程序和运行程序，进行蒙特卡洛（MC）模拟，判断是否达到平衡，进而计算g(r)、S(Q)和经验势的步骤。程序继续验证经验势是否接近于0，进行EPMC步骤。最后，当达到平衡后，会累积g(r)和S(Q)，并终止程序。

图 3（原文 Figure 3）展示了 NeuDATool 的算术流程图。

4.2.1. 步骤 I：定义原子、分子和模拟盒

这一步通过 C++ 的继承机制来定义模拟所需的物理实体。设计了三个基本的 C++ 类来帮助用户定义其特定的子类或对象：

• Atom 类：用于定义模拟中的原子类型。它定义了基本的原子属性，如原子名称、元素名称、同位素名称、坐标位置、中子散射长度等。原子坐标位置使用高能物理类库 (CLHEP) 的 Hep3Vector 类来定义，因为 Hep3Vector 提供了丰富的函数来执行平移、旋转、距离和角度计算。用户通常在分子对象中初始化原子，而不是定义新的原子类。
• Molecule 类：用于定义分子、分子内势及其运动，例如平移、旋转等。用户需要通过继承该类来定义新的分子类，并可以尝试针对其分析的特殊分子间和分子内运动。这个类使得程序非常灵活和用户友好。
• SimBox 类：用于生成模拟盒。用户需要为其模拟定义一个子类。在该子类中，用户只需编辑一个初始函数，其中包含已定义的分子。用户可以使用任何合适的算法将分子放置在模型盒中，以生成初始构象。

4.2.2. 步骤 II：重新编译并运行程序

在定义了所需的原子、分子和模拟盒类后，用户需要重新编译 C++ 程序，然后运行模拟。

4.2.3. 步骤 III：蒙特卡洛 (MC) 模拟

程序使用 参考势 (RP) 进行 MC 模拟，直到系统达到平衡。程序会按顺序或随机移动分子或原子。模拟盒的总势能变化 ($\Delta U = U_{after} - U_{before}$) 被用作接受移动的判据。

• 如果 $\Delta U < 0$，移动被接受。
• 如果 $\Delta U > 0$，移动以概率 $e^{-\Delta U/kT}$ 被接受。

4.2.4. 步骤 IV：计算 g(r)、S(Q) 和 EMP

当 MC 模拟达到平衡时，程序会执行以下计算：

• 计算中子结构因子差异 $\Delta S(Q)$： 这是模拟得到的结构因子 $S_{sim}(Q)$ 与实验测量得到的结构因子 $S_{exp}(Q)$ 之间的差值。
• 计算经验势 (EMP)： EMP 是通过对 $\Delta S(Q)$ 进行傅里叶逆变换得到的。
• 更新势能： 程序将计算得到的 EMP 与 RP 相加，形成更新后的总势能，用于接下来的经验势蒙特卡洛 (EPMC) 模拟。

4.2.5. 步骤 V：经验势蒙特卡洛 (EPMC) 模拟

EPMC 算法与 MC 算法非常相似，不同之处在于：

• 当模拟达到平衡时，程序会再次计算 EMP。
• 将新的 EMP 添加到之前的势能中，然后使用更新后的势能继续进行模拟。
• 这个迭代过程会重复进行，直到 EMP 的值变得非常低，表明模拟结构与实验数据高度一致。

4.2.6. 步骤 VI：累积 g(r) 和 S(Q)

当 EMP 达到非常低的值后，程序开始累积模拟数据以提高统计精度，直到获得平滑的径向分布函数 g(r) 和结构因子 S(Q) 曲线。除了 g(r) 和 S(Q)，程序还会输出包含所有原子坐标的文件，格式兼容 GROMACS (.gro) 或 LAMMPS (.xyz)。这些文件可以输入到 GROMACS 或 LAMMPS 中进行焓、熵等计算，也可以使用 Visual Molecular Dynamics (VMD) 进行可视化。用户还可以通过添加新的输出函数，定义、计算和输出任何与样品原子结构直接相关的感兴趣变量。

4.3. 加速方法

程序中计算量最大的部分是原子对计数函数 $\Delta N(r)$ 和 N(r) 的计算，它们的计算复杂度大致与模拟盒中原子数量的平方成正比。鉴于模拟盒中的原子数量通常超过一万个，这些计算在每个 MC 模拟步骤（用于 $\Delta U$ 计算）或每个 MC/EPMC 平衡点（用于 g(r) 和 S(Q) 更新）都需要重新计算。因此，对这些算法进行加速是提高程序整体性能的关键。

4.3.1. GPU 硬件加速

图形处理单元 (GPU) 具有数千个并行线程，非常适合加速 $\Delta N(r)$ 和 N(r) 的计算。NeuDATool 利用 CUDA C 支持 GPU 加速。

该图像是CUDA C代码，用于计算N(r)，涉及原子之间的距离计算和数据存储。代码中定义了全局函数Fill_NR()，处理原子坐标数组并根据距离范围填充配对数组。

图 4（原文 Figure 4）展示了用于计算 N(r) 的 CUDA C 内核函数代码。该代码定义了一个全局函数 Fill_NR()，它接收原子坐标数组 d_pos、原子数量 atom_num、模拟盒大小 box_size、最大距离 R_max、距离步长 dR、原子对计数数组 d_NR 和原子类型对数 atom_pair_num 作为输入。该内核通过并行线程计算每个原子对在不同距离下的贡献，并填充 d_NR 数组，从而利用 GPU 的大规模并行能力加速了 N(r) 的计算。

4.3.2. 跨节点并行

为了实现跨计算机集群节点的并行配置，程序使用了 mpich2。mpich2 是基于消息传递接口 (MPI) 标准的实现，支持不同节点之间的点对点 (point-to-point) 和集体 (collective) 数据通信，因此在该程序中具有很高的效率。

4.3.3. 共享内存多线程并行

为了在单个计算机或服务器节点内部实现共享内存、多线程并行配置，程序使用了 OpenMP。OpenMP 语法支持动态设置线程数量，这在大多数情况下线程数量无法预先确定的场景中非常方便。

4.3.4. 综合并行策略

NeuDATool 结合使用 OpenMP 和 MPI API，将不同的 $\Delta N(r)/N(r)$ 计算任务分配给属于同一计算机节点的不同 GPU 卡，或分配给计算机集群中的不同节点。此外，为了程序的广泛适用性，还编辑了 CPU 串行算法版本。

5. 实验设置

为了测试 NeuDATool 的正确性和计算速度，研究人员使用了完全氢化水、完全氘化水、半氘化水和无定形二氧化硅 ($\mathrm{SiO_2}$) 样品在环境温度和压力下的中子散射数据。

5.1. 数据集

实验使用了来自 ISIS 网站上随 EPSR 和 GudRun 发布的完全氢化水、完全氘化水、半氘化水和 $\mathrm{SiO_2}$ 样品的实验中子数据。

5.1.1. 水样品模拟设置

• 分子定义： 定义一个继承自 Molecule 基类的 C++ 类来描述 $\mathrm{H_2O}$ 分子。
• 分子内势： 使用谐振子势 (harmonic oscillator potentials) 维持每两个原子之间的分子结构。
• 随机运动： 采用三种不同的随机运动模式：
  1. $\mathrm{H_2O}$ 分子的平移 (translation)。
  2. $\mathrm{H_2O}$ 分子的旋转 (rotation)。
  3. 原子（H 或 O）的平移。
• 经验势计算： 水样品中原子类型对（O-O、O-H 和 H-H）的经验势 $U_j^{EMP}(r)$ 使用以下公式计算： $$U _ { j } ^ { E M P } ( r ) = \frac { 1 } { 4 \pi \rho } \sum _ { i = 1 , M } w _ { j i } ^ { - 1 } \int _ { 0 } ^ { \infty } \Delta S _ { i } ( Q ) e ^ { - i Q r } d Q$$ 其中：
  • $U_j^{EMP}(r)$ 是不同原子对的经验势，如水中的 O-O、O-H 和 H-H。
  • $\rho$ 是原子数密度。
  • $M$ 是中子散射图谱的数量。
  • $w_{ji}^{-1}$ 是权重因子矩阵的逆矩阵。在水样品中，中子散射图谱的数量和原子类型对的数量都是3，因此可以直接使用 CLHEP 提供的 Invert() 函数计算逆矩阵。
  • $\Delta S_i(Q)$ 是第 $i$ 个样品模拟与实验之间中子结构因子的差异。
  • $Q$ 是散射矢量，$r$ 是实空间距离。
  • $\int_{0}^{\infty} \Delta S_i(Q) e^{-iQr} dQ$ 表示傅里叶逆变换。

5.1.2. $\mathrm{SiO_2}$ 样品模拟设置

• 分子定义： 定义两个 C++ 类来描述 Si 和 O 作为单原子分子。
• 经验势计算差异： 与水不同，$\mathrm{SiO_2}$ 只有一个中子散射图谱，但有三种不同的原子类型对（Si-Si、Si-O 和 O-O）。
• 权重因子矩阵处理： 借鉴 EPSR 的方法，使用一个 $3 \times 3$ 单位矩阵将 $(1 \times 3)$ 的权重因子矩阵 $w_{ij}$ 扩展为一个 $(4 \times 3)$ 的矩阵 w'_{ij}。然后，使用蒙特卡洛方法计算伪逆矩阵 $w'^{-1}$。这种处理方式允许在只有单个实验散射谱的情况下，通过引入额外的约束来解耦出不同的原子对势。

5.2. 评估指标

5.2.1. 正确性 (Correctness)

• 中子散射谱比较： 将 NeuDATool 模拟得到的结构因子 S(Q) 与实验测量的 S(Q) 进行直接比较。通过目视检查曲线的吻合程度来评估模型是否能准确复现实验数据。
• 径向分布函数 (PDF) 比较： 将 NeuDATool 模拟得到的 PDF 曲线与已发表的，经过验证的液体结构研究结果（如水在 Kusalik、Head-Gordon 和 Soper 等人研究中的 PDF 分布）进行比较，以验证重建微观结构的物理合理性。

5.2.2. 计算速度 (Computational Speed)

• 每秒步数 (steps/sec)： 这是一个直接衡量模拟吞吐量的指标，表示在给定时间内程序可以执行多少个模拟步骤。更高的 steps/sec 值表示更快的计算速度。通过比较在不同原子数量下，CPU 串行、MPI 并行和 GPU 加速算法的 steps/sec，来量化不同加速方法的性能提升。

5.3. 对比基线

5.3.1. 正确性验证

• 实验数据： 将模拟结果与实际中子散射实验数据进行比较。
• 现有研究： 将模拟得到的径向分布函数与已发表的、广泛接受的关于水结构的研究结果进行对比。

5.3.2. 计算速度验证

• CPU 串行算法： 作为性能的基准，衡量在单个 CPU 核心上不使用任何并行或加速技术时的计算速度。
• MPI 并行算法： 衡量在多个 CPU 核心或跨多个计算机节点上使用 MPI 进行并行计算时的速度提升。
• GPU 加速算法： 衡量在配备 GPU 的硬件上，通过 CUDA 优化后的算法所能达到的最高计算速度。

5.3.3. 硬件环境

用于测试速度性能的实验平台是一个小型计算机集群，其配置如下：

• 操作系统： CentOS 7.3。
• 节点数量： 2 个节点。
• 每个节点配置：
  • CPU： 两颗 Intel Xeon Scalable Gold 6126 CPU（SkylakeSP 架构，12 核，24 线程，2.6 GHz，睿频 3.7 GHz），配备 19.25 MB L3 Intel 智能缓存。
  • GPU： 两块 Nvidia Tesla V100 计算 GPU 卡。
  • 内存： 128 GB 双倍数据速率 (DDR4) 错误纠正码 (ECC) 注册共享内存。
• 节点互联： 两个节点通过 InfiniBand (IB) 连接器互联，数据传输速度可达 56 Gb/s，确保了节点间高效的数据通信。

6. 实验结果与分析

6.1. 核心结果分析

6.1.1. 正确性验证

论文首先通过水和无定形二氧化硅样品的中子散射数据测试了 NeuDATool 的正确性。

该图像是图表，展示了NeuDATool与实验样品的中子散射谱比较。左侧为_ ext{H}_2 ext{O}样品，右侧为无定形ext{SiO}_2样品。图中点表示实验数据，实线表示NeuDATool模拟结果，虚线为随机初始模拟框。

图 5（原文 Figure 5）展示了 NeuDATool 模拟的中子散射谱与实验样品的比较。左侧是 $\mathrm{H_2O}$ 样品，右侧是无定形 $\mathrm{SiO_2}$ 样品。图中的点表示实验数据，实线表示 NeuDATool 的模拟结果，虚线则表示随机初始模拟盒的结果。从图中可以看出，NeuDATool 的模拟结果（实线）与实验数据（点）吻合良好，无论是对于水样品还是二氧化硅样品，都能够准确地复现实验观测到的散射谱，而初始的随机模拟盒（虚线）则与实验数据存在显著差异。这有力地证明了 NeuDATool 能够基于中子衍射谱正确地重建实验样品的原子结构。

该图像是图表，展示了水的PDF分布，包括O-H、H-H和O-O的相应曲线，并且插图展示了液态水的分子结构。图中标注的数值代表了不同的原子间距。

图 6（原文 Figure 6）展示了 NeuDATool 模拟的水的 O-H、H-H 和 O-O 径向分布函数 (PDF) 曲线，插图为液态水结构的卡通示意图。模拟结果显示，液态水呈现出四面体随机网络结构，平均每个水分子与中心水分子形成约 3.5 个氢键。这些 PDF 分布与 Kusalik、Head-Gordon 和 Soper 等人的研究结果一致。例如，O-H 峰值在约 1.0 Å 和 1.8 Å，H-H 峰值在约 1.5 Å 和 2.3 Å，O-O 峰值在约 2.8 Å 和 4.5 Å 附近，这些都符合液态水的典型结构特征。这进一步证实了 NeuDATool 不仅能重现宏观散射谱，还能在微观层面重建出准确的原子间关联。

综上所述，所有这些结果都证明了 NeuDATool 能够根据中子衍射谱正确地重建实验样品的原子结构。

6.1.2. 计算速度验证

论文还对 NeuDATool 在不同加速方法下的计算速度进行了测试。

该图像是图表，展示了不同原子数量下，GPU加速、MPI并行和CPU串行算法的计算速度比较。横坐标表示模拟盒中的原子数量，纵坐标表示每秒计算步骤数，体现了GPU在较大原子数量时的显著性能优势。

图 7（原文 Figure 7）展示了在不同原子数量的模拟盒下，串行算法、MPI 并行算法和 GPU 加速算法的计算速度比较。X 坐标表示模拟盒中的原子数量，Y 坐标表示每秒的模拟步骤数。从图中可以看出，随着原子数量的增加，GPU 加速算法的性能优势显著提升，其每秒步数远高于 MPI 并行和串行算法。

以下是原文 Table 1 的结果：

表 1（原文 Table 1）详细列出了在不同原子数量下，GPU 加速、MPI 并行和 CPU 串行算法的模拟速度（单位：步/秒）。

分析：

• 小规模系统 (3 × 10^3 原子)： 在原子数量较少时，CPU 串行、MPI 和 GPU 的速度差异不显著，甚至 MPI 略快于 GPU (2000 vs 1562)。这可能是因为 GPU 在初始化和数据传输上存在开销，对于小任务无法完全发挥其并行优势。
• 中等规模系统 (3 × 10^4 至 2.5 × 10^5 原子)： 随着原子数量增加，GPU 的优势开始凸显。
  • 在 3 × 10^4 原子时，GPU (9507) 远超 CPU (189.8) 约 50 倍，也远超 MPI (334.9) 约 28 倍。
  • 在 1 × 10^5 原子时，GPU 达到了最高速度 12412 步/秒，而 CPU 仅为 62.67 步/秒，MPI 为 110.5 步/秒。此时 GPU 比 CPU 串行快约 198 倍，比 MPI 快约 112 倍。
  • 在 2.5 × 10^5 原子时，GPU 速度为 10189 步/秒，CPU 串行速度为 24.30 步/秒。此时 GPU 比 CPU 串行快约 419 倍。这个数据与摘要中“超过 400 倍”的声明相符。
• 大规模系统 (1.2 × 10^6 至 3 × 10^7 原子)：
  • MPI 和 CPU 串行算法在原子数量超过 2.5 × 10^5 后未提供数据，这可能意味着其计算时间过长或内存不足以进行如此大规模的模拟。
  • GPU 算法在百万级原子数量时仍然能够运行，并保持相对较高的速度。例如，在 1.2 × 10^6 原子时，速度为 3917 步/秒；在 3 × 10^7 原子时，速度为 340.3 步/秒。这表明 GPU 能够处理远超 CPU/MPI 独立运行所能处理的系统规模。

结论：

• MPI 确实能够根据节点或线程数量提高速度。
• GPU 提供了“卓越的加速比”，并且是实现超大规模模拟的关键。最重要的是，通过 GPU 加速，程序可以模拟包含超过 100 万个原子的系统。这是一个实质性的改进，因为它允许程序模拟大于 200 Å 的系统，从而有望在未来分析全原子模型中的大分子样品。

7. 总结与思考

7.1. 结论总结

本研究成功地开发并提出了 NeuDATool，一个用于中子散射数据分析的开源框架。该工具使用面向对象语言 C++ 编写，极大地增强了程序的灵活性，方便用户定义特殊的分子和蒙特卡洛随机运动模式。原子力场的势函数和相应参数可以通过修改 C++ 头文件和源文件进行轻松修改或添加。此外，C++ 作为一种易读的高级计算机语言，使用户能够尝试新的算法和程序流程，以改进分析，并计算和输出任何重要变量。

NeuDATool 不仅支持 OpenMP 在服务器内部节点进行并行，还支持跨计算机集群不同节点之间的并行（通过 MPI），并集成了 GPU 硬件加速技术。特别是 GPU 加速，显著提高了计算速度，使得程序未来有能力分析全原子模型中的无序大分子样品和纳米粒子。通过对水和无定形二氧化硅的测试，NeuDATool 证明了其能够准确重建样品微观结构的能力，并在约10万个原子的模拟盒中，相比 CPU 串行算法实现了超过400倍的加速。

7.2. 局限性与未来工作

尽管 NeuDATool 提供了强大的计算能力和用户灵活性，但作者也坦诚指出了其目前的局限性，并展望了未来的工作：

• 测试样本有限： 目前程序仅用有限数量的控制样品进行了测试。这表明其在更广泛、更复杂材料系统上的通用性和鲁棒性仍需进一步验证。
• 非完全功能软件包： 在当前形式下，NeuDATool 并非一个完全功能的软件软件包，更像是一个工具包框架。这意味着用户可能需要进行额外的开发才能满足其特定需求。
• 未来展望： 作者希望将其作为开源工具包框架发布供感兴趣的科学家使用。他们设想用户能够贡献众多新的分子类、算法和分析例程，以使程序在未来变得更加强大和完善。这体现了开源社区协作开发的潜力。

7.3. 个人启发与批判

7.3.1. 个人启发

这篇论文提供了一些重要的启发：

• 旧有科学软件现代化的重要性： 许多经典且经过验证的科学算法（如 EPSR）是用旧语言编写的，限制了它们的性能和可扩展性。NeuDATool 的成功案例表明，通过采用现代编程语言（C++）、面向对象设计和最新硬件加速技术（GPU/MPI），可以显著提升这些算法的价值和应用范围。这对于其他领域的遗留代码现代化具有借鉴意义。
• 面向对象设计的优势： 在科学计算中，面向对象编程的灵活性和可扩展性是巨大的优势。通过 Atom、Molecule、SimBox 等类的设计，NeuDATool 使得用户能够更容易地定制和扩展模型，这对于探索新的物理化学系统和开发新算法至关重要。
• 异构计算的强大潜力： 结合 CPU（用于控制和串行任务）、GPU（用于计算密集型并行任务）和多节点并行（MPI），能够突破单一计算模式的瓶颈，实现前所未有的计算能力，为处理多尺度、大规模的科学问题提供了有效的解决方案。
• 开源在科学发展中的作用： 将工具作为开源框架发布，鼓励社区贡献，是加速科学软件发展和促进跨学科合作的有效途径。

7.3.2. 批判与潜在改进

尽管 NeuDATool 取得了显著进展，仍有一些潜在问题或可以改进的地方：

• 初始功能集的吸引力： 论文提到它“不是一个完全功能的软件包”，这可能在初期劝退一些不熟悉 C++ 编程或不愿投入大量二次开发时间的潜在用户。一个更健壮、涵盖更广泛材料和功能的初始版本，可能会加速其被社区采纳和贡献的速度。
• 大规模系统验证的深度： 尽管论文展示了在百万级原子系统中的速度提升，但对于这些超大规模系统，除了速度，其模拟结果的物理准确性、稳定性以及可能出现的新的计算挑战（例如，长程相互作用的处理、边界效应、数据存储和传输）的详细分析较少。在未来工作中，可以对大规模系统的结构特征进行更深入的分析，以进一步验证其可靠性。
• 与其他现代 RMC/EPSR 实现的比较： 论文主要将 NeuDATool 与原始 Fortran EPSR 的性能进行了比较。如果在当时有其他基于 C++ 或 GPU 加速的 RMC/EPSR 算法实现，进行更广泛的基准测试可以更全面地评估 NeuDATool 在同类工具中的竞争力。
• 用户友好性与文档： 虽然 C++ 提供了灵活性，但对于非专业的程序员科学家来说，直接修改 C++ 头文件和源文件可能仍然具有挑战性。未来可以考虑提供更高级别的脚本接口（如 Python 绑定）或更详细的用户文档和教程，以进一步降低使用门槛。
• 误差分析和收敛性： 论文详细介绍了速度提升和结果正确性，但对 EMP 收敛标准的具体量化、收敛速度、以及不同参数设置（如温度、势能更新频率）对模拟结果和效率的影响讨论较少。更严谨的误差分析和收敛性研究可以增强方法的可靠性。