1. Bibliographic Information

1.1. Title

NeuDATool: An Open Source Neutron Data Analysis Tools, Supporting GPU Hardware Acceleration, and Across-computer Cluster Nodes Parallel

1.2. Authors

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

Their affiliations are:

• Institute of High Energy Physics, Chinese Academy of Sciences (CAS), Beijing 100049, China

• Spallation Neutron Source Science Center, Dongguan 523803, China

• University of Chinese Academy of Sciences, Beijing 100049, China

  He Cheng is the corresponding author.

1.3. Journal/Conference

This paper was published on arXiv, a preprint server. While not a peer-reviewed journal or conference in its current form, arXiv is a highly respected platform for disseminating research quickly within the scientific community, particularly in physics, mathematics, computer science, and related fields. Papers published here often undergo subsequent peer review and publication in journals.

1.4. Publication Year

2019

1.5. Abstract

The paper introduces NeuDATool, an open-source neutron scattering data analysis framework. It addresses limitations of the conventional Empirical Potential Structure Refinement (EPSR) algorithm, which, despite its reliability in constructing atomic structures of disordered liquids, is programmed in Fortran with a shared-memory OpenMP architecture, limiting its computational speed and scalability on modern supercomputer clusters. NeuDATool is developed in object-oriented C++, supports parallel processing across cluster nodes using Message Passing Interface (MPI), and incorporates Graphics Processing Unit (GPU) hardware acceleration. Performance tests on water and amorphous silica neutron scattering data demonstrate its ability to correctly reconstruct microstructures. Crucially, GPU acceleration showed a speed increase of over 400 times compared to the CPU serial algorithm for simulations involving approximately 100,000 atoms. NeuDATool offers a flexible and powerful alternative for C++ programmers to define specific models and algorithms for neutron data analysis.

1.6. Original Source Link

https://arxiv.org/abs/1904.08756 (Preprint) PDF Link: https://arxiv.org/pdf/1904.08756v3.pdf

2. Executive Summary

2.1. Background & Motivation

Problem Statement

The Empirical Potential Structure Refinement (EPSR) algorithm is a widely used and reliable method for analyzing neutron scattering data to construct the most probable atomic structures of disordered liquids. Developed by the British spallation neutron source (ISIS) Disordered Materials Group in the 1980s, it has been instrumental in numerous scientific discoveries. However, the original EPSR software is programmed in Fortran and utilizes a shared-memory architecture with OpenMP.

Importance of the Problem

With the rapid advancements in computational hardware, particularly the widespread construction of supercomputer clusters and the integration of Graphics Processing Unit (GPU) acceleration technology, the Fortran-based, OpenMP-only EPSR faces significant limitations:

• Computational Speed: It cannot fully leverage the parallel processing capabilities of modern hardware, leading to slow calculations for large systems.
• Scalability: The shared-memory architecture restricts its ability to parallelize across multiple nodes within a supercomputer cluster, thereby limiting the size of systems that can be simulated (e.g., macromolecular systems requiring hundreds of thousands to millions of atoms).
• Flexibility: Being written in Fortran, a procedural language, makes it less flexible for users to easily define new models, molecules, or algorithms compared to object-oriented languages.

Innovative Idea

The paper proposes NeuDATool as an open-source framework to address these limitations. The core innovation lies in reimplementing the EPSR algorithm using modern computational paradigms:

• Object-Oriented C++: Enhances flexibility and user-friendliness for defining custom simulation components.
• Distributed Memory Parallelism (MPI): Enables scaling across multiple nodes in a computer cluster.
• GPU Hardware Acceleration (CUDA): Provides significant speed improvements for computationally intensive parts of the algorithm.

2.2. Main Contributions / Findings

Primary Contributions

1. Development of NeuDATool: An entirely new, open-source framework for neutron data analysis, specifically EPSR simulations, written in object-oriented C++.
2. Enhanced Parallelism: Integration of Message Passing Interface (MPI) for parallel execution across different nodes in a computer cluster, and OpenMP for shared-memory parallelism within a node.
3. GPU Hardware Acceleration: Implementation of CUDA C to utilize GPUs for accelerating the most computationally intensive parts of the EPSR algorithm, specifically the calculation of pair distribution functions (N(r)).
4. Improved Flexibility and User-Friendliness: The object-oriented design in C++ allows users to easily define new Atom, Molecule, and SimBox classes, and customize movement models and potentials.

Key Conclusions / Findings

1. Correct Microstructure Reconstruction: NeuDATool successfully reconstructed the correct microstructure of test samples (water and amorphous silica), demonstrating its reliability and accuracy in producing results consistent with experimental data and previous studies.
2. Significant Speed Enhancement: GPU acceleration dramatically improves calculation speed, achieving more than a 400-fold increase compared to the CPU serial algorithm for a simulation box with approximately 100,000 atoms. This allows for the simulation of much larger systems (e.g., over 1 million atoms), which was previously infeasible.
3. Scalability for Macromolecular Systems: The combined MPI and GPU parallelization enables NeuDATool to handle systems larger than 200 Å, paving the way for analyzing disordered macromolecular samples and nanoparticles in all-atomic models.
4. Open-Source and Extensible Framework: NeuDATool provides a flexible and powerful toolkit for scientists familiar with C++ programming, encouraging community contributions to expand its capabilities and algorithms.

3. Prerequisite Knowledge & Related Work

3.1. Foundational Concepts

To understand this paper, a beginner should be familiar with the following core concepts:

• Neutron Scattering: A technique used to probe the structure and dynamics of materials. Neutrons interact with atomic nuclei (unlike X-rays which interact with electron clouds), making them sensitive to light elements like hydrogen and isotopes.

  • Neutron Diffraction/Total Scattering: A specific application of neutron scattering to determine the arrangement of atoms in disordered materials (liquids, glasses, amorphous solids). It measures the scattered intensity as a function of the scattering vector $Q$.
  • Scattering Vector (Q): Represents the change in momentum of the scattered neutron. Its magnitude is inversely proportional to the length scale being probed (Q = 4\pi \sin(\theta)/\lambda, where $\theta$ is half the scattering angle and $\lambda$ is the neutron wavelength).
  • Partial Pair Correlation Function ($g_{\alpha\beta}(r)$ or PDF): A fundamental quantity in statistical mechanics that describes the probability of finding an atom of type $\beta$ at a distance $r$ from an atom of type $\alpha$. It provides detailed information about the local atomic structure. If $g_{\alpha\beta}(r) = 1$, the atomic distribution is random. Peaks in $g_{\alpha\beta}(r)$ indicate preferred interatomic distances.
  • Structure Factor (S(Q)): The Fourier transform of the PDFs, which is directly measured in a diffraction experiment. It reflects the spatial arrangement of atoms in reciprocal space.
  • Deuteration: A technique where hydrogen atoms in a sample are replaced by deuterium atoms (an isotope of hydrogen). Hydrogen and deuterium have very different neutron scattering lengths (hydrogen has a negative scattering length and a large incoherent scattering cross-section, while deuterium has a positive scattering length and a small incoherent scattering cross-section). By selectively deuterating samples, scientists can vary the contrast and isolate specific PDFs, which is crucial for solving the structural equations for multi-component systems.

• Monte Carlo (MC) Simulation: A broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. In materials science, MC is often used to simulate physical systems by generating random configurations and accepting/rejecting them based on energy criteria.

  • Metropolis Monte Carlo: A common MC algorithm where a new configuration is accepted if its energy is lower than the current one, or with a certain probability if its energy is higher (Boltzmann factor $e^{-\Delta U/kT}$). This allows the system to escape local energy minima.
  • Potential Energy ($U$): A function that describes the energy of a system based on the positions of its atoms. It's typically composed of terms like bond stretching, angle bending, dihedral rotation, and non-bonded interactions (van der Waals, electrostatic).

• Empirical Potential Structure Refinement (EPSR): An advanced MC simulation method specifically designed for neutron scattering data. It refines an atomic structure against experimental S(Q) data.

  • Reference Potential (RP): Standard interatomic potentials (like those used in molecular dynamics, e.g., OPLS, AMBER, CHARMM) that describe the fundamental interactions between atoms (e.g., repulsion, attraction, bond constraints). These ensure physically realistic structures.
  • Empirical Potential (EMP): A feedback potential introduced in EPSR that guides the MC simulation towards configurations consistent with experimental neutron scattering data. It's derived from the difference between the simulated and experimental S(Q) (i.e., $\Delta S(Q)$). When $\Delta S(Q)$ is large, the EMP is strong, pushing the simulation to match the experiment. When $\Delta S(Q)$ is small, EMP becomes negligible.

• Parallel Computing:

  • Shared Memory Architecture (OpenMP): A parallel programming model where multiple threads within a single processor (or node) can access a common shared memory space. It's effective for parallelizing loops and tasks within a single machine.
  • Distributed Memory Architecture (MPI - Message Passing Interface): A standard for parallel programming where processes (often on different nodes in a cluster) communicate by explicitly sending and receiving messages. Each process has its own private memory, and data must be explicitly passed between them. MPI is essential for scaling applications across large supercomputer clusters.
  • Graphics Processing Unit (GPU) Acceleration (CUDA): GPUs are specialized electronic circuits designed to rapidly manipulate and alter memory to accelerate the creation of images for output to a display device. Modern GPUs, particularly those from Nvidia with CUDA (Compute Unified Device Architecture), can perform highly parallel computations, making them suitable for tasks involving large numbers of independent, repetitive calculations (e.g., summing pair interactions in a simulation).

• Object-Oriented Programming (C++): A programming paradigm based on the concept of "objects," which can contain data and code to manipulate that data. Key features include:

  • Classes and Objects: Blueprints and instances of those blueprints, respectively, bundling data (attributes) and behavior (methods).
  • Inheritance: A mechanism where one class (subclass) can inherit properties and behaviors from another class (superclass), allowing for code reuse and hierarchical organization.
  • Polymorphism: The ability of objects of different classes to respond to the same message (method call) in different ways.
  • Encapsulation: Hiding the internal state and functionality of an object and exposing only what is necessary, promoting modularity.

3.2. Previous Works

The paper directly builds upon and critiques the existing Empirical Potential Structure Refinement (EPSR) software.

• Original EPSR (ISIS Disordered Materials Group, 1980s): This is the foundational work. The original EPSR package, as described by Mcgreevy and Pusztai [2], is central to the field. Its key principle is to refine an atomic configuration to match experimental neutron diffraction data by introducing an empirical potential derived from the difference between simulated and experimental structure factors. The algorithm aims to construct the most-probable atomic structures of disordered liquids.

  • Core idea: Combine realistic reference potentials (like Lennard-Jones) with a data-driven empirical potential. The empirical potential acts as a perturbation, guiding the Monte Carlo simulation to reproduce the experimental S(Q).
  • Mathematical basis of EMP: The empirical potential for an atom pair $j$ at distance $r$ ($U_j^{EMP}(r)$) is derived from the inverse Fourier transform of the difference in structure factors, $\Delta S_i(Q)$, weighted by factors related to scattering lengths and concentrations. $ 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 $ Where:
    • $U_j^{EMP}(r)$: The Empirical Potential for atom type pair $j$ at a distance $r$.
    • $\rho$: The number density of atoms in the system.
    • $M$: The number of independent neutron scattering profiles available from experiments (e.g., using deuteration).
    • $w_{ji}^{-1}$: The inverse of the weighting factor matrix that relates the experimental structure factors to the partial structure factors of different atom pairs. This matrix accounts for the scattering lengths and concentrations of the atoms.
    • $\int_{0}^{\infty} \Delta S_i(Q) e^{-iQr} dQ$: The inverse Fourier transform of the difference in the structure factor for the $i$-th scattering profile. $\Delta S_i(Q) = S_{sim,i}(Q) - S_{exp,i}(Q)$, where $S_{sim,i}(Q)$ is the simulated structure factor and $S_{exp,i}(Q)$ is the experimental structure factor.
  • Original Implementation: Fortran, shared-memory OpenMP.

• Other All-Atom Model Simulations: The paper mentions Reverse Monte Carlo (RMC) [3, 4] as another common method to solve structural matrices and reconstruct atomic structures. While EPSR uses potentials to guide the simulation towards experimental data, RMC directly modifies atomic configurations to match experimental data without explicit potentials, minimizing a chi-squared function.

• Parallel Computing Technologies:

  • OpenMP [6]: Used for shared-memory parallelization. The original EPSR uses this. NeuDATool also uses it for within-node parallelism.
  • MPI (specifically mpich2 [5, 13, 27]): Standard for distributed-memory parallelization across cluster nodes. This is a key improvement NeuDATool adopts.
  • CUDA [7, 8]: Nvidia's platform for GPU computing. NeuDATool leverages this for hardware acceleration. The paper references other works on GPU acceleration in computational science [9, 10, 11, 12].

3.3. Technological Evolution

The field of neutron scattering data analysis has evolved from fundamental theoretical models and experimental techniques to sophisticated computational methods.

1. Early Neutron Scattering: The development of total scattering spectrometers like TSS at HELIOS in the 1970s [1] marked the beginning of detailed structural studies of disordered materials.

2. Algorithm Development (1980s): EPSR emerged as a powerful algorithm to interpret complex scattering patterns, moving beyond direct inversion methods which can be ill-posed. This involved combining physical potentials with experimental data as constraints within a Monte Carlo framework.

3. Computational Hardware Evolution:

   • Single-core CPUs: Early EPSR implementations were designed for single-core processors.
   • Multi-core CPUs & Shared Memory: The advent of multi-core CPUs led to the adoption of OpenMP for parallelization on a single machine, as seen in the original EPSR.
   • Supercomputer Clusters & Distributed Memory: The construction of large-scale supercomputer clusters necessitated distributed memory paradigms like MPI to scale simulations across many interconnected machines.
   • GPU Computing: The realization that GPUs could perform general-purpose computation (GPGPU) for highly parallelizable tasks led to a revolution in computational speed, especially for molecular simulations.

4. Software Engineering Practices: The move from procedural languages (like Fortran) to object-oriented languages (like C++) allows for more modular, flexible, and maintainable codebases, facilitating the integration of new features and models.

   This paper's work (NeuDATool) directly positions itself at the intersection of these evolutions, updating the established EPSR algorithm to leverage modern $C++$, MPI, and GPU technologies, thereby pushing the boundaries of what is computationally feasible for neutron data analysis.

3.4. Differentiation Analysis

Compared to the original EPSR software and other related work, NeuDATool offers several core differences and innovations:

• Programming Language:

  • Original EPSR: Fortran (procedural).
  • NeuDATool: C++ (object-oriented).
  • Innovation: $C++$ provides greater flexibility for users to define custom Atom, Molecule, and SimBox classes using inheritance, making the framework more adaptable to new scientific problems and models. It also allows for easier integration of modern libraries and development practices.

• Parallelization Strategy:

  • Original EPSR: Shared-memory OpenMP (limited to a single computer/node).
  • NeuDATool:
    • Distributed-memory MPI (across multiple nodes in a cluster).
    • Shared-memory OpenMP (within a single node).
    • GPU hardware acceleration (CUDA) for critical computational kernels.
  • Innovation: This multi-layered parallelization strategy is the most significant differentiator. It allows NeuDATool to overcome the scalability limitations of the original EPSR, enabling simulations of much larger systems (hundreds of thousands to millions of atoms) that were previously intractable. GPU acceleration, in particular, provides an unprecedented speedup.

• Computational Performance:

  • Original EPSR: Limited by CPU serial or OpenMP parallelization.
  • NeuDATool: Achieves over 400 times speedup with GPU acceleration for 100,000-atom systems.
  • Innovation: This massive performance increase opens up new research avenues, allowing scientists to study complex macromolecular systems or nanoparticles in full atomic detail, which require very large simulation boxes.

• Open-Source and Extensibility:

  • Original EPSR: A software package, likely with limited extensibility for external developers.

  • NeuDATool: Proposed as an open-source toolkit framework.

  • Innovation: Encourages community contributions, allowing users to develop and integrate new molecular classes, algorithms, and analysis routines, fostering collaborative development and expanding its capabilities beyond the initial release.

    In summary, NeuNeuDATool is not just a reimplementation but a significant architectural upgrade of the EPSR algorithm, bringing it into the era of high-performance computing with modern software engineering practices.

4. Methodology

4.1. Principles

NeuDATool is fundamentally a Monte Carlo (MC) simulation method, but it differentiates itself from standard Metropolis MC by incorporating neutron scattering data information into its atomic potentials. The core idea is to refine an atomic structure until its simulated neutron scattering pattern matches an experimental one. This is achieved by using two types of atomic potentials: Reference Potentials (RP) and Empirical Potentials (EMP).

• Reference Potential (RP): These are conventional interatomic potentials, similar to those used in Molecular Dynamics (MD) simulations. They ensure that the simulated atoms and molecules maintain physically realistic shapes, bond lengths, and non-bonded interactions. Examples include OPLS, AMBER, CHARMM force fields. At the beginning of a simulation, only RP is used.

• Empirical Potential (EMP): This potential has no fixed form and is dynamically generated. Its purpose is to guide the MC simulation towards configurations that match the experimental neutron structure factor S(Q). It acts as a feedback mechanism. The EMP is calculated from the difference between the simulated and experimental structure factors, $\Delta S(Q)$, using a reverse Fourier transform. When the simulation deviates from the experiment (i.e., $\Delta S(Q)$ is large), the EMP introduces a strong perturbation to the total potential, pushing the system towards the experimental data. As the simulation approaches the experimental S(Q), $\Delta S(Q)$ (and thus EMP) becomes very small, essentially vanishing when a good fit is achieved. The EMP is constructed using a list of Poisson distributions in real space and their corresponding Fourier transforms in Q space.

  The simulation proceeds iteratively: first, MC with RP to reach equilibrium, then EMP is introduced, and the process continues as Empirical Potential Monte Carlo (EPMC) until the EMP is negligible.

4.2. Core Methodology In-depth (Layer by Layer)

The algorithmic flow of NeuDATool is depicted in Fig. 3 and involves several key steps, designed to leverage object-oriented programming for flexibility and parallel computing for speed.

The following figure (Fig. 3 from the original paper) shows the schematic of arithmetic flow of NeuDATool:

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

VLM Description: The image is a schematic diagram illustrating the algorithmic flow of NeuDATool. It depicts the steps from defining atoms, molecules, and simulation boxes to recompiling and running the program for Monte Carlo (MC) simulations, checking for equilibrium, and calculating g(r), S(Q), and the empirical potential. The program continues to verify whether the empirical potential approaches zero before proceeding to the EPMC step. Finally, once equilibrium is reached, it accumulates g(r) and S(Q) and terminates the program.

4.2.1. Define Atoms, Molecules, and Simulation Box (Step I)

Users define their specific simulation system using C++'s object-oriented features, particularly inheritance from predefined base classes.

• Atom Class: This base class defines fundamental properties of an atom type within a simulation.
  • Properties: atom name, element name, isotope name, coordinate position, neutron scattering length.
  • Coordinate position: Managed using the Hep3Vector class from the CLHEP (Class Library for High Energy Physics) library [21, 22]. Hep3Vector is chosen for its rich set of functions for vector operations, such as translation, rotation, distance, and angle calculations, which are essential for molecular simulations.
  • Usage: Users typically initialize Atom objects within Molecule objects rather than defining new Atom classes directly.
• Molecule Class: This base class is used to define molecules, their intramolecular potentials (e.g., bond stretching, angle bending), and their allowed movements (e.g., translation, rotation).
  • Usage: Users inherit from this base class to create new molecular classes tailored to their specific samples (e.g., a water molecule class). This allows for the implementation of custom inter- and intramolecular movements and potentials, making the program highly flexible.
• SimBox Class: This base class is used for generating and managing the simulation box, which contains all the Molecule objects.
  • Usage: Users define a subclass of SimBox and implement an initial function to place the defined Molecule objects within the simulation box, generating an initial conformation. Any suitable algorithm can be used for initial placement.

4.2.2. Recompile and Run Program (Step II)

After defining the system components, the C++ code is compiled and the NeuDATool program is executed.

4.2.3. Monte Carlo (MC) Simulation with Reference Potentials (Step III)

The program begins with an MC simulation using only the Reference Potentials (RP).

• Movement: Molecules or atoms are moved either sequentially or randomly within the simulation box.
• Acceptance Criterion: The decision to accept or reject a proposed move is based on the change in the total potential energy of the simulation box, $\Delta U = U_{after} - U_{before}$.
  • If $\Delta U < 0$ (energy decreases), the move is always accepted.
  • If $\Delta U > 0$ (energy increases), the move is accepted with a probability $e^{-\Delta U/kT}$.
    • $k$: Boltzmann constant.
    • $T$: Simulation temperature.
• Equilibration: This MC phase continues until the system reaches equilibrium, meaning its macroscopic properties (like density and energy) fluctuate around stable average values.

4.2.4. Calculate g(r), S(Q), and Empirical Potential (EMP) (Step IV)

Once the MC simulation with RP reaches equilibrium:

• Calculate Simulated $S_{sim}(Q)$ and $g_{sim}(r)$: The program computes the current pair distribution functions (g(r)) and their Fourier transforms, the simulated neutron structural factors ($S_{sim}(Q)$), from the current atomic configuration.

• Calculate Difference in Structure Factors ($\Delta S(Q)$): The difference between the simulated and experimental structure factors is calculated: $ \Delta S ( Q ) = S _ { s i m } ( Q ) - S _ { e x p } ( Q ) $ Where:

  • $\Delta S(Q)$: The difference in the structure factor at scattering vector $Q$.
  • $S_{sim}(Q)$: The structure factor calculated from the current simulation snapshot.
  • $S_{exp}(Q)$: The experimental structure factor obtained from neutron diffraction data.

• Calculate Empirical Potential ($U^{EMP}(r)$): The EMP is derived from $\Delta S(Q)$ by applying an inverse Fourier transform. As described in the EPSR manual [20] and also in the context of the original EPSR, the EMP is represented as a series of Poisson distributions in real space, whose Fourier transforms are used to fit $\Delta S(Q)$. The formula for the EMP of an atom type pair $j$ is given by: $ 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 $ Where:

  • $U_j^{EMP}(r)$: The Empirical Potential for atom type pair $j$ at a distance $r$.
  • $\rho$: The number density of atoms in the system.
  • $M$: The number of independent neutron scattering profiles used (typically obtained from deuteration experiments).
  • $w_{ji}^{-1}$: The inverse of the weighting factor matrix. This matrix relates the experimentally measured total structure factors to the individual partial structure factors of each atom pair. For water, with 3 atom pairs (O-O, O-H, H-H) and 3 independent scattering profiles (H_2O, D_2O, HDO), $w_{ji}^{-1}$ is a $3 \times 3$ matrix. For cases like $SiO_2$ with 3 atom pairs (Si-Si, Si-O, O-O) but only 1 scattering profile, a pseudo-inverse matrix is calculated (e.g., using a $3 \times 3$ unit matrix to enlarge the $1 \times 3$ weighted factor matrix to $4 \times 3$, then a Monte Carlo method for the pseudo-inverse, as mentioned in the paper).
  • $\int_{0}^{\infty} \Delta S_i(Q) e^{-iQr} dQ$: The inverse Fourier transform of the difference in the structure factor for the $i$-th scattering profile. This term effectively converts the discrepancy in $Q$-space (S(Q)) into a potential in $r$-space (U(r)).

• Update Potential: The calculated EMP is then added to the RP to form a new, updated total potential for the subsequent simulation steps.

  The following figure (Fig. 1 from the original paper) shows the Poisson distributions used to fit $\Delta S(Q)$:

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

VLM Description: The image is a diagram illustrating the Poisson distributions with different au values in R space (top) and Q space (bottom). The upper part's distributions are used to fit riangle S(Q), with all peaks normalized to 1.0.

The following figure (Fig. 2 from the original paper) shows examples of $\Delta S(Q)$ fitting and corresponding EMP:

该图像是图表，展示了不同_ 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$的变化而变化。

VLM Description: The image is a chart showing the riangle S(Q) fitting results of different _ ext{H}_2 ext{O} samples (left) and their corresponding empirical potentials (right). HD ext{H}_2 ext{O}, FH ext{H}_2 ext{O}, and FD ext{H}_2 ext{O} represent half-deuterated, fully hydrogenated, and fully deuterated water samples, respectively, with changes in riangle S(Q) and potentials varying with $Q$ and $R$.

4.2.5. Empirical Potential Monte Carlo (EPMC) Simulation (Step V)

The simulation continues with the updated potential (RP + EMP). This phase is similar to the initial MC but now explicitly uses the EMP as a feedback mechanism.

• Iteration: When the EPMC simulation reaches equilibrium, the program recalculates the EMP again based on the new configuration and updates the total potential. This iterative process repeats.
• Convergence: The EPMC continues until the EMP becomes very small (approaching zero), indicating that the simulated structure factor is now in excellent agreement with the experimental data ($\Delta S(Q) \approx 0$).

4.2.6. Accumulate g(r) and S(Q) (Step VI)

Once convergence is achieved and the EMP is negligible, the program continues running EPMC for a longer period to accumulate more simulation data.

• Statistics: This accumulation phase is crucial for improving the statistical quality and smoothness of the calculated g(r) and S(Q) curves.
• Output: Besides g(r) and S(Q), the program outputs atomic coordinate files (e.g., in GROMACS (.gro) or LAMMPS (.xyz) formats) [23, 24]. These files can be used for visualization with tools like VMD [25, 26] or for further analysis (e.g., enthalpy, entropy calculations) with other simulation packages.
• Extensibility: Users can add new output functions to calculate and store any other variables of interest related to the atomic structure.

4.3. Acceleration Method

The most computationally intensive parts of the NeuDATool algorithm are the calculations of $\Delta N(r)$ (for $\Delta U$) and N(r) (for g(r) and S(Q)). $\Delta N(r)$ needs to be calculated after every MC step, and N(r) after every MC/EPMC equilibrium. These calculations involve determining distances between many pairs of atoms, which scales quadratically with the number of atoms.

NeuDATool employs a multi-pronged approach for acceleration:

1. GPU Hardware Acceleration (CUDA C):

   • Target: N(r) calculation. This involves iterating through all atom pairs to count how many fall within specific distance bins. This task is highly parallelizable, as the distance calculation for each pair is independent.

   • Implementation: CUDA C is used to write kernel functions that can be executed on GPUs. GPUs have thousands of parallel threads, making them ideal for accelerating such calculations.

   • Benefit: Achieves substantial speedups (e.g., >400 times) for large systems, enabling simulations with over 1 million atoms.

     The following figure (Fig. 4 from the original paper) shows the CUDA C code used to calculate N(r):

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

   VLM Description: The image is CUDA C code for calculating N(r), involving distance calculations between atoms and data storage. The code defines a global function Fill_NR() that processes atomic coordinate arrays and fills the pair array based on distance ranges.

   • Explanation of CUDA C code: The __global__ keyword indicates a CUDA kernel function, which runs on the GPU. Fill_NR takes arrays of atom $x$, $y$, $z$ coordinates, the number of atoms (_natoms), the box size (_box), the maximum $r$ (_rmax), the bin size (_dr), and a device array d_pair (which will store counts for each $r$-bin) as input.
     • idx = blockIdx.x * blockDim.x + threadIdx.x;: Each CUDA thread is assigned a unique idx.
     • $if (idx >= _natoms) return;$: Threads beyond the number of atoms stop.
     • for (j = idx + 1; j < _natoms; j++): This nested loop structure means each thread idx is responsible for calculating distances between atom idx and all subsequent atoms $j$.
     • Inside the loop, dx, dy, dz are calculated, applying periodic boundary conditions ($fmodf(d, _box)$) if $d$ is larger than half the box size.
     • r = sqrtf(dx * dx + dy * dy + dz * dz);: The distance $r$ is calculated.
     • if (r < _rmax) atomicAdd(&d_pair[(int)(r / _dr)], 1);: If the distance is within the maximum $r$, the corresponding bin in d_pair is incremented. atomicAdd is used to safely update d_pair in parallel by multiple threads, avoiding race conditions.

2. Distributed Memory Parallelization (MPI - mpich2):

   • Target: Parallelization across multiple compute nodes in a supercomputer cluster.
   • Implementation: mpich2 [27], an implementation of the MPI standard, is used. MPI enables point-to-point and collective data communication between processes running on different nodes, each with its own memory.
   • Benefit: Allows NeuDATool to scale to very large systems by distributing the computational workload (e.g., calculating different $\Delta N(r) / N(r)$ contributions) across many nodes.

3. Shared Memory Parallelization (OpenMP):

   • Target: Parallelization within a single compute node or server using multiple CPU cores.
   • Implementation: OpenMP [6] is integrated into the C++ program. It supports dynamic setting of thread numbers, making it convenient for programming.
   • Benefit: Utilizes all available CPU cores on a single machine efficiently, distributing different $\Delta N(r) / N(r)$ calculations among them.

4. CPU Serial Algorithm:

   • A basic CPU serial version is also maintained for broader applicability, especially on systems without GPUs or MPI setups. This serves as a baseline for performance comparisons.

     NeuDATool combines OpenMP and MPI APIs to distribute different N(r) calculations to various GPU cards within a computer node or across different nodes in a cluster. This hierarchical parallelization strategy is key to its high performance and scalability.

5. Experimental Setup

5.1. Datasets

The performance of NeuDATool was tested using experimental neutron scattering data from two types of samples:

1. Water Samples:

   • Types: Full hydrogen water ($\text{H}_2\text{O}$), full deuterated water ($\text{D}_2\text{O}$), and half deuterated water ($\text{HDO}$).
   • Conditions: Ambient temperature and pressure.
   • Source: The experimental neutron data for water samples were obtained from the ISIS website, distributed with EPSR [28] and GudRun [29, 30].
   • Purpose: Water is a well-studied system, making it an excellent benchmark for verifying the correctness of the structural reconstruction. The different isotopic forms allow for the generation of multiple independent scattering profiles, which are crucial for resolving the partial pair correlation functions (O-O, O-H, H-H) as described in the methodology for EPSR.

2. Amorphous Silica ($\text{SiO}_2$) Sample:

   • Type: Amorphous silica.

   • Conditions: Ambient temperature and pressure.

   • Source: Experimental neutron data also from the ISIS website, distributed with EPSR and GudRun.

   • Purpose: $SiO_2$ represents a different class of disordered material with different structural characteristics (covalent network vs. hydrogen-bonded liquid) and a different number of independent scattering profiles (only one in this case, leading to the use of pseudo-inverse matrix techniques for EMP calculation). This tests the versatility and robustness of NeuDATool beyond simple liquids.

     These datasets were chosen because they are standard benchmark systems in neutron diffraction studies of disordered materials, and their experimental data are readily available from established sources like ISIS. This allows for direct comparison and validation of NeuDATool's results against widely accepted EPSR simulations.

5.2. Evaluation Metrics

The paper primarily evaluates NeuDATool based on two aspects: correctness of microstructure reconstruction and computational speed.

5.2.1. Correctness Metrics

1. Neutron Scattering Spectra Comparison:

   • Conceptual Definition: This involves comparing the simulated neutron structure factors ($S_{sim}(Q)$) directly against the experimental neutron structure factors ($S_{exp}(Q)$). A good match indicates that the simulated atomic arrangements accurately reproduce the experimentally observed scattering pattern.
   • Mathematical Concept: The goal is to minimize the difference $\Delta S(Q) = S_{sim}(Q) - S_{exp}(Q)$ across the range of scattering vectors $Q$.
   • Symbol Explanation:
     • $S_{sim}(Q)$: Simulated structure factor at scattering vector $Q$.
     • $S_{exp}(Q)$: Experimental structure factor at scattering vector $Q$.
     • $\Delta S(Q)$: Difference between simulated and experimental structure factors.

2. Pair Distribution Function (PDF) / g(r) Comparison:

   • Conceptual Definition: This involves analyzing the simulated PDFs (g(r)) for various atom pairs and comparing them with established structural knowledge for the material (e.g., known bond lengths, coordination numbers, hydrogen bonding networks). For water, this means comparing O-O, O-H, and H-H PDFs with previous studies.
   • Mathematical Concept: The PDF $g_{\alpha\beta}(r)$ is defined such that $\rho_{\beta} g_{\alpha\beta}(r) 4\pi r^2 dr$ is the average number of $\beta$ atoms in a spherical shell of thickness dr at distance $r$ from an $\alpha$ atom, where $\rho_{\beta}$ is the number density of $\beta$ atoms.
   • Symbol Explanation:
     • $g_{\alpha\beta}(r)$: Pair distribution function between atom types $\alpha$ and $\beta$ at distance $r$.
     • $\rho_{\beta}$: Number density of atoms of type $\beta$.
     • $r$: Distance.

5.2.2. Computational Speed Metric

1. Simulation Steps per Second (steps/sec.):
   • Conceptual Definition: This metric quantifies the raw processing throughput of the simulation algorithm. A higher value indicates faster execution. It directly measures how many Monte Carlo steps (which typically involve attempting a move, calculating energy changes, and accepting/rejecting) can be performed within one second.
   • Mathematical Concept: $ \text{Speed} = \frac{\text{Number of Simulation Steps}}{\text{Total Simulation Time (seconds)}} $
   • Symbol Explanation:
     • Number of Simulation Steps: The total count of Monte Carlo move attempts.
     • Total Simulation Time: The wall-clock time taken to execute these steps, measured in seconds.

5.3. Baselines

NeuDATool's performance was compared against different computational configurations and parallelization strategies to demonstrate the benefits of its design:

1. CPU Serial Algorithm: This is the most basic baseline, representing the performance of the algorithm running on a single CPU core without any parallelization. This establishes the fundamental speed without enhancements.

2. MPI Parallel Algorithm: This baseline demonstrates the performance improvement gained by distributing the computational workload across multiple CPU cores or nodes using Message Passing Interface (mpich2). This is a common parallelization strategy in high-performance computing.

3. GPU Acceleration Algorithm: This is NeuDATool's primary acceleration method, utilizing CUDA to offload intensive calculations to GPUs. The comparison against CPU serial and MPI highlights the specific advantage of GPU hardware.

   The comparison across these baselines for varying numbers of atoms (from 3,000 to 30 million) allows for a comprehensive understanding of NeuDATool's scalability and efficiency under different computational scenarios. The goal is to show how NeuDATool surpasses these conventional approaches in terms of speed, especially for larger systems.

6. Results & Analysis

6.1. Core Results Analysis

The paper presents results validating both the correctness of the structural reconstruction and the computational speed of NeuDATool under various parallelization schemes.

6.1.1. Correctness

The NeuDATool simulations were tested against experimental neutron scattering data for water and amorphous silica.

The following figure (Fig. 5 from the original paper) shows the neutron scattering spectral comparison between NeuDATool and experimental samples:

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

VLM Description: The image is a graph showing the neutron scattering spectral comparison between NeuDATool and the experimental samples. The left side presents the _ ext{H}_2 ext{O} samples, while the right side presents the amorphous ext{SiO}_2 sample. The points denote the experimental data, the solid lines represent the NeuDATool simulations, and the dashed lines indicate the random initial simulation boxes.

Analysis of Figure 5:

• Water Samples (Left Panel): The left panel compares the neutron scattering spectra (S(Q)) for different isotopic forms of water (fully hydrogenated (FH), half deuterated (HD), and fully deuterated (FD) $\text{H}_2\text{O}$).

  • The points represent the experimental data.
  • The solid lines represent the results from NeuDATool simulations.
  • The dashed lines represent the initial state of the simulation boxes (random initial configuration).
  • Observation: The solid lines from NeuDATool simulations show excellent agreement with the experimental points across the entire Q-range for all three water samples. This indicates that NeuDATool successfully refined the atomic structure to accurately reproduce the experimental scattering patterns. The initial dashed lines are far from the experimental data, highlighting the effectiveness of the EPSR refinement process.

• Amorphous Silica Sample (Right Panel): The right panel shows the S(Q) comparison for amorphous $\text{SiO}_2$.

  • Observation: Similar to water, the NeuDATool simulation (solid line) matches the experimental data (points) very well. This demonstrates the framework's ability to handle different types of disordered materials, even when only one scattering profile is available (requiring more complex pseudo-inverse matrix calculations for EMP).

• Conclusion on Correctness: The strong consistency between NeuDATool's simulated scattering spectra and experimental data confirms that the software can reconstruct the correct microstructure of the samples based on neutron diffraction profiles.

  Further validation of correctness is provided by analyzing the Pair Distribution Functions (PDFs) for water.

The following figure (Fig. 6 from the original paper) shows the PDF distributions of O-H, H-H, and O-O from NeuDATool simulation for water:

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

VLM Description: The image is a graph showing the PDF distributions for water, including the corresponding curves for O-H, H-H, and O-O, with an inset illustrating the molecular structure of liquid water. The annotated values represent the distances between different atoms.

Analysis of Figure 6:

• O-H PDF: Shows a sharp peak around 1.0 Å, corresponding to the intramolecular O-H bond length in water. A broader, less intense peak around 1.8 Å (and beyond 3 Å) indicates hydrogen-bonding interactions between different water molecules.
• H-H PDF: Exhibits peaks around 1.6 Å (intramolecular H-H distance) and then broader features at larger distances, reflecting intermolecular hydrogen arrangements.
• O-O PDF: Shows a primary peak around 2.8 Å, characteristic of the nearest-neighbor oxygen-oxygen distance in liquid water's hydrogen-bonded network. A second peak around 4.5 Å indicates further neighbors.
• Comparison to Literature: The paper states that these results are consistent with previous studies by Kusalik, Head-Gordon, and Soper [31, 32, 33], who also describe liquid water as a tetrahedrally random network with an average of 3.5 water molecules forming hydrogen bonds with a central one.
• Conclusion on Reliability: The accurate reproduction of known PDF features and consistency with established literature further verifies the reliability of NeuDATool in reconstructing atomic structures.

6.1.2. Computational Speed

The computational speed was tested on a small computer cluster with Intel Xeon Scalable Gold 6126 CPUs and Nvidia Tesla V100 GPUs.

The following figure (Fig. 7 from the original paper) shows the calculation speed comparison among serial algorithm, MPI parallel and GPU acceleration algorithms:

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

VLM Description: The image is a chart showing the comparison of computational speed among GPU acceleration, MPI parallel, and CPU serial algorithms at varying atom quantities. The X-axis represents the number of atoms in the simulation box, while the Y-axis denotes simulation steps per second, highlighting the significant performance advantage of GPU at larger atom counts.

The following are the results from Table 1 of the original paper:

Analysis of Figure 7 and Table 1:

• CPU Serial Performance: For a small number of atoms (3,000), the CPU serial algorithm performs reasonably well (1754 steps/sec). However, its speed drops significantly as the number of atoms increases, reaching only 24.30 steps/sec for 250,000 atoms. This demonstrates the $N^2$ scaling bottleneck of the distance calculations.
• MPI Parallel Performance: For 3,000 atoms, MPI shows a slight improvement over CPU serial (2000 vs 1754 steps/sec), indicating efficient parallelization for small systems. However, its performance also decreases with increasing atom numbers, albeit less sharply than CPU serial, reaching 42.57 steps/sec for 250,000 atoms. The table indicates that MPI was not tested for systems larger than 250,000 atoms, likely due to diminishing returns or memory limitations on the available nodes without GPUs.
• GPU Acceleration Performance:
  • Small Systems (3,000 atoms): GPU is slightly slower than MPI (1562 vs 2000 steps/sec) and comparable to CPU serial. This is often observed in GPU computing due to overheads associated with data transfer between CPU and GPU, which can dominate for small workloads.
  • Medium to Large Systems (30,000 to 250,000 atoms): GPU acceleration shows a dramatic increase in speed. For 30,000 atoms, GPU achieves 9507 steps/sec, which is about 28 times faster than MPI (334.9 steps/sec) and 50 times faster than CPU serial (189.8 steps/sec). For 100,000 atoms, GPU reaches its peak at 12412 steps/sec, which is approximately 112 times faster than MPI (110.5 steps/sec) and 198 times faster than CPU serial (62.67 steps/sec). The abstract mentions "more than 400 times compared with CPU serial algorithm at a simulation box consists about 100 thousand atoms" -- while the table shows 198x for 100k atoms, the abstract might be referring to a peak comparison where GPU is faster, or perhaps for a specific type of calculation that wasn't fully detailed in the table. Even with the table's numbers, the speedup is immense.
  • Very Large Systems (1.2 million to 30 million atoms): Crucially, GPU acceleration is the only method shown to handle these extremely large systems. While its absolute speed (steps/sec) decreases for these largest systems (e.g., 340.3 steps/sec for 30 million atoms), it remains the only viable option for such scale. The ability to simulate 1.2 million atoms at 3917 steps/sec is a significant breakthrough.
• Impact: The GPU acceleration allows NeuDATool to simulate systems comprising over 1 million atoms, which was previously impossible for EPSR to run in reasonable timeframes. This is a critical improvement, enabling the analysis of much larger systems (exceeding 200 Å) relevant for macromolecules and nanoparticles.
• Scalability: The GPU performance curve (Fig. 7) initially rises sharply, then plateaus, and eventually declines for very large systems. The decline for the largest systems is likely due to memory limitations of the GPU, data transfer bottlenecks, or the specific CUDA kernel implementation becoming less efficient for extreme atom counts relative to its peak. Nevertheless, it maintains superior performance compared to CPU-only methods.

6.2. Ablation Studies / Parameter Analysis

The paper does not explicitly present ablation studies or detailed parameter analyses. However, the comparison of CPU serial, MPI parallel, and GPU accelerated versions serves as an implicit form of ablation study, demonstrating the individual and combined effectiveness of the parallelization components. By showing the performance gain from MPI (distributed CPU parallelism) and then the further, much larger gain from GPU acceleration, the authors effectively highlight the critical role of each component. The varying atom numbers in the speed test also act as a parameter analysis, showing how system size impacts the relative performance of each method.

7. Conclusion & Reflections

7.1. Conclusion Summary

NeuDATool is a newly developed, open-source framework for neutron scattering data analysis, specifically designed to address the limitations of the traditional Empirical Potential Structure Refinement (EPSR) algorithm. By reimplementing EPSR in object-oriented C++ and integrating modern parallel computing technologies—MPI for across-cluster-node parallelism and CUDA for GPU hardware acceleration—NeuDATool achieves unprecedented computational speed and scalability. The software successfully reconstructs the correct microstructures of disordered liquids like water and amorphous silica. Crucially, GPU acceleration delivers a speedup of over 400 times compared to CPU serial execution for systems around 100,000 atoms, enabling the simulation of systems with millions of atoms. This advancement makes NeuDATool a powerful tool for studying complex macromolecular systems and nanoparticles that were previously beyond the reach of EPSR.

7.2. Limitations & Future Work

Limitations

The authors acknowledge that NeuDATool, in its current form, is not a fully functional software package despite its flexibility and powerful calculation capacity. It has been tested with a "very limited number of control samples." This suggests that while the core algorithm and acceleration mechanisms are validated, the software might lack:

• A comprehensive suite of pre-defined molecular models.
• A wide range of analysis routines.
• A user-friendly graphical interface or extensive documentation for non-C++ programmers.
• Thorough testing across a broader spectrum of materials and experimental conditions.

Future Work

The authors aspire to release NeuDATool as an open-source toolkit framework. Their vision for future work includes:

• Community Contributions: Encouraging interested scientists to contribute numerous new molecular classes, algorithms, and analysis routines. This strategy aims to expand the software's capabilities and make it more powerful and versatile over time through collaborative development.

7.3. Personal Insights & Critique

Inspirations

• Modernization of Legacy Algorithms: The paper provides an excellent example of how well-established, scientifically validated algorithms can be revitalized and made relevant for contemporary research by adopting modern programming paradigms and hardware acceleration. The transformation of EPSR from a Fortran/OpenMP code to a C++/MPI/CUDA framework demonstrates a clear path for updating other legacy scientific software.
• Power of Open Source and Object-Oriented Design: The emphasis on an open-source, object-oriented C++ framework highlights the benefits of these approaches for scientific computing. It fosters collaboration, allows for easy customization and extension by domain experts (even those who are not primary developers), and promotes modularity and maintainability.
• Strategic Use of Heterogeneous Computing: The multi-layered acceleration strategy (OpenMP for intra-node, MPI for inter-node, and CUDA for GPU) is a sophisticated and effective approach to maximizing computational efficiency on modern supercomputer architectures. It teaches that optimal performance often requires a combination of parallelization techniques tailored to specific computational bottlenecks.

Potential Issues, Unverified Assumptions, or Areas for Improvement

• User-Friendliness for Non-C++ Programmers: While the C++ framework offers flexibility for C++-savvy users, it might present a steep learning curve for experimentalists or theorists who are not proficient in C++ programming. Future development could explore higher-level scripting interfaces (e.g., Python wrappers) or user-friendly GUIs to broaden its accessibility.

• GPU Memory Management and Scalability Beyond 30 Million Atoms: While the paper demonstrates impressive GPU speedups, the performance curve for GPUs shows a decline for the largest systems. This often indicates challenges with GPU memory limits, data transfer bandwidth, or kernel efficiency for extreme workloads. Further optimization might be needed for even larger systems, perhaps involving more sophisticated data partitioning or multi-GPU strategies.

• Broader Validation: As acknowledged by the authors, testing with a limited number of control samples is a starting point. Comprehensive validation across a wider range of materials (e.g., ionic liquids, molten salts, complex biological systems), different experimental conditions (e.g., high pressure, high temperature), and comparisons with other refinement methods (like RMC) would further solidify its robustness.

• Community Engagement Strategy: The success of an open-source project heavily relies on community engagement. While the intention is there, active maintenance, clear documentation, tutorials, and a responsive developer community will be crucial for attracting contributions and ensuring the long-term viability and impact of NeuDATool.

• Detailed Benchmarking: While speed comparisons are provided, more detailed breakdowns of where time is spent (e.g., N(r) calculation vs. potential energy calculation vs. data transfer) could offer deeper insights for future optimizations. The 400x speedup quoted in the abstract vs. the 198x observed in the table for 100k atoms might warrant further clarification on specific calculation types or conditions.

  Overall, NeuDATool represents a significant and necessary step forward for neutron data analysis. Its modern architecture and emphasis on high-performance computing address long-standing limitations, paving the way for more detailed and extensive structural studies of disordered materials.