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

TL;DR

NeuDATool is an open source neutron data analysis software that leverages GPU acceleration and cluster parallelism to significantly speed up structure refinement calculations, facilitating advanced analysis of disordered materials.

Contribution

It introduces a C++ based, GPU-accelerated, and cluster-parallel framework for EPSR, enhancing calculation speed and flexibility over traditional Fortran implementations.

Findings

GPU acceleration increases speed by over 400 times.

Successfully reconstructs microstructures of water and silica.

Supports parallel processing across cluster nodes.

Abstract

Empirical potential structure refinement (EPSR) is a neutron scattering data analysis algorithm and a software package. It was developed by the British spallation neutron source (ISIS) Disordered Materials Group in 1980s, and aims to construct the most-probable atomic structures of disordered liquids. It has been extensively used during the past decades, and has generated reliable results. However, it is programmed in Fortran and implements a shared-memory architecture with OpenMP. With the extensive construction of supercomputer clusters and the widespread use of graphics processing unit (GPU) acceleration technology, it is now necessary to rebuild the EPSR with these techniques in the effort to improve its calculation speed. In this study, an open source framework NeuDATool is proposed. It is programmed in the object-oriented language C++, can be paralleled across nodes within a computer cluster, and supports GPU acceleration. The performance of NeuDATool has been tested with water and amorphous silica neutron scattering data. The test shows that the software could reconstruct the correct microstructure of the samples, and the calculation speed with GPU acceleration could increase by more than 400 times compared with CPU serial algorithm at a simulation box consists about 100 thousand atoms. NeuDATool provides another choice for scientists who are familiar with C++ programming and want to define specific models and algorithms for their analyses.

Tables

Table 1. Table 1: Simulation speeds with GPU acceleration, MPI parallel and CPU serial algorithms [steps/sec.].

Atomic number	$3\times {10}^3$	$3\times {10}^4$	$1\times {10}^5$	$2.5\times {10}^5$	$1.2\times {10}^6$	$3\times {10}^6$	$1\times {10}^7$	$3\times {10}^7$
GPU	1562	9507	12412	10189	3917	1963	757.0	340.3
MPI	2000	334.9	110.5	42.57	—	—	—	—
CPU	1754	189.8	62.67	24.30	—	—	—	—

Equations

$$S(Q)=\sum _{\alpha =1}{c}_{\alpha }{b}_{\alpha }^2+\sum _{\alpha =1,\beta \ge \alpha }(2-{\delta }_{\alpha \beta })c_{\alpha }{c}_{\beta }{b}_{\alpha }{b}_{\beta }\{4\pi \rho {\int }_0^{\infty }{r}^2({\text{g}}_{\alpha \beta }(r)-1)\frac{\sin (Qr)}{Qr}dr\}$$

$$U_j^{EMP}(r)=\frac{1}{4\pi \rho }\sum _{i=1,M}{w}_{ji}^{-1}{\int }_0^{\infty }\Delta S_i(Q)e^{-iQr}dQ$$

Full text

++

††journal: arXiv.org

1 Introduction

Neutron total scattering for disordered materials is a powerful tool to study the most probable atomic structures in an amorphous system. Since the introduction of the first total scattering spectrometer (TSS)[1], which was developed at HELIOS in the 1970s, numerous important scientific problems have been solved with the use of these types of instruments, including the elucidation of the structures of water at high- and low-densities, and the observation of the heterogeneities in mixed alcohol–aqueous solutions.

The success of neutron total scattering for disordered materials is based on deuteration techniques. The scattering pattern $S(Q)$ is a weighted summation of the Fourier transforms of all the pair correlation functions (PDF),

[TABLE]

where $Q$ is the scattering vector, $\sum_{\alpha}c_{\alpha}b_{\alpha}^{2}$ is a flat background, and $(2-\delta_{\alpha\beta})c_{\alpha}c_{\beta}b_{\alpha}b_{\beta}$ is the weighted factors of different partial structural factors, $c_{\alpha/\beta}$ and $b_{\alpha/\beta}$ are the atom ratio and scattering length of atoms with types $\alpha/\beta$, respectively. $g_{\alpha\beta}(r)$ is the PDF between atom types $\alpha$ and $\beta$. If the different atom type number is $M$ in a sample, there are $M(M+1)/2$ different $g_{\alpha\beta}(r)$ functions in Eq. (1), and we must solve all of them first before the most probable atomic structure is obtained. It is thus difficult to reveal the most probable all-atom structure of disordered material, just according to one scattering curve. Fortunately, neutron scattering scientist has deuteration technique. Because deuterated samples have almost the same atomic structure as their hydrogenate forms, each deuterated sample can thus generate a different scattering pattern.

All-atom model simulations, such as EPSR and RMC, are common methods to solve the matrix of **Eq.(1)**s, and reconstruct the atomic structure of the tested samples[2, 3, 4]. Because EPSR can provide a reliable and visualized atomic microstructure, it has been extensively used. However, EPSR still contains constraints which limit its potential applications. First it is programmed in Fortran. Fortran can implement the execution of algorithms with very fast calculation speeds, but it is a procedure-oriented language. Second, EPSR is paralleled in shared memory architecture with OpenMP[5, 6]. It cannot be paralleled across different nodes of a supercomputer cluster. This restricts the calculation speed and the analysis system scale, such as in the case of a macromolecular system. A macromolecule generally contains hundreds of atoms. To cover the entire range of scattering vectors that characterize the entire set of conformation states of macromolecules in a typical total scattering instrument, the simulation box composed of more than half a million atoms should be larger than 10nm. EPSR cannot run in such a large system. Therefore, parallel calculations are necessary.

In light of this, the object-oriented language C

++

and the compute unified device architecture (CUDA)[7, 8] are used to develop a toolkit NeuDATool. In NeuDATool, users can define easily a new simulation box, atoms, molecules, and movement models using the class multiple inheritance mechanism. Graphics processing unit (GPU) hardware acceleration[9, 10, 11, 12] is supported by CUDA C, and this allows the program to take advantage of commonly used GPU computing servers. In addition, with the distributed memory architecture API message passing interface (MPI) mpich2[5, 13], NeuDATool can be paralleled across nodes of a supercomputer cluster. These acceleration techniques make the program have huge speed promotion.

2 NeuDATool simulation method and program Flow

2.1 Algorithmic Principle

NeuDATool is essentially a Monte-Carlo simulation method, but its difference from metropolis MC is based on the fact that the atomic potentials used include neutron scattering data information[2, 14]. In NeuDATool, the atomic potential used in the MC simulation is divided into two categories i.e., “reference potential (RP)” and “empirical potential (EMP)”[2]. RP is similar to that used in molecular dynamics (MD) simulations, so the potential form and parameters can be obtained from all-atom MD force fields, such as the optimized potential for liquid simulations (OPLS)[15, 16], assisted model building with energy refinement (AMBER)[17], chemistry at HARvard molecular mechanics (CHARMM)[18], and others[2, 19]. By contrast, EMP has no fixed form and it is used to reflect the differences of neutron structural factors between experiments and MC simulations ($\Delta S(Q)$). To be exact, EMP is the reverse Fourier transform of $\Delta S(Q)$ as the perturbation to the RP to guide the simulation approach to scattering measurements. In the present form of NeuDATool, EMP is expressed based on a list of Poisson distributions in real space and their corresponding Fourier transforms in Q space[20]. Fig.1 shows a series Poisson distributions with different $\lambda$ values in R space and Q space (top and bottom) respectively. Here, $\lambda$ is the mathematical expectation of the distribution. In NeuDATool, the Poisson series in Q space are used to fit $\Delta S(Q)$. Fig.2 gives examples of the $\Delta S(Q)$ fitting results of H2O samples’ simulation (left) and its corresponding EMP potential (right) respectively. In short, RP is used to assign molecules with reasonable shapes and realize other constraints whose correctness have been approved, while EMP is used as a feedback parameter to lead the MC simulation in a consistently progressive manner with experimental data. In the simulations, only RP is used in the MC simulation at the beginning, and the potential changes of the system ($\Delta U$) are used as the selection criteria of molecules’ or atoms’ random movements. When the simulation reaches equilibrium, EMP is introduced to fit $\Delta S(Q)$ and is added to RP to continue the simulation. When the MC simulation with updated potential reaches equilibrium again, NeuDATool calculates EMP and updates the simulation potential once more. This process is repeated until the EMP becomes equal to zero, i.e., $\Delta S(Q)$ becomes very small.

2.2 NeuDATool Algorithmic Flow

The algorithmic flow of NeuDATool is shown in Fig.3. Details of some processes are described below:

1. I.

Define atoms, molecules, and the simulation box with inheritance in C

++

. Three basic C

++

 classes have been designed to help users to define their special sub-classes or objects:

$\bullet$ Atom: It is used to define an atom type in a simulation. Some basic atomic properties, such as the atom name, element name, isotope name, coordinate position, neutron scattering length, are defined. In the class, the coordinate position of an atom is defined using the Hep3Vector class of the class library for high energy physics (CLHEP)[21, 22] because the Hep3Vector has abundant functions to perform transition, rotation, distance, and angle calculations. Users need to initialize atoms in molecular objects rather than define new atom classes.

$\bullet$ Molecule: It is used to define molecules, intramolecular potentials and their movements, such as translation, rotation, etc. Users need to define new molecular classes through inheritance and can try special inter- and intramolecular movements for their analyses. This class makes the program very flexible and user-friendly.

$\bullet$ SimBox: It is used for generating a simulation box. Users need to define a subclass for their simulation. In the subclass, users only need to edit an initial function with defined molecules. Users can use any suitable algorithm to place the molecules in the model box for generating an initial conformation. 2. II.

Recompile and Run Program. 3. III.

MC: The program performs MC simulations with reference potentials until equilibrium is reached. The program moves molecules or atoms in sequence or randomly. The potential energy variation of the simulation box ($\Delta U=U_{after}-U_{before}$) is used as the movement acceptance criterion. If $\Delta U<0$, the movement is accepted. If $\Delta U>0$, the movement is accepted with a probability $e^{-\Delta U/kT}$. 4. IV.

Calculate $g(r)$, $S(Q)$, and EMP: When MC reaches equilibrium, the program calculates the difference of the neutron structural factor $\Delta S(Q)$ between the simulation $S_{sim}(Q)$ and experiment $S_{exp}(Q)$. EMP is calculated by applying the Fourier transform to $\Delta S(Q)$. The program adds the EMP and RP together as the updated potential to perform the EPMC simulation. 5. V.

EPMC: The empirical potential Monte Carlo (EPMC) algorithm is very similar to MC with the exception that when the simulation reaches equilibrium, the program calculates the EMP again, adds it to the previous potential, and performs the simulation with the updated potential. When the EMP attains a very low value, the program starts to accumulate simulation data. 6. VI.

Accumulate $g(r)$ and $S(Q)$: The program still perform EPMC to accumulate simulation data to improve statistics until smooth $g(r)$ and $S(Q)$ curves are obtained. With the exception of $g(r)$ and $S(Q)$, the program outputs a coordinate file which includes all the atoms in a text format, as used in the GROningen machine for chemical simulations (GROMACS) (with the suffix of .gro)[23], or as used in the large-scale atomic/molecular massively parallel simulator (LAMMPS) (with the suffix of .xyz)[24]. Accordingly, this file can be input to GROMACS or LAMMPS to calculate enthalpy, entropy, etc., and also can be visualized with visual molecular dynamics (VMD)[25, 26]. Users can define, calculate, and output any interesting variables which are related directly with the atomic structure of the sample by adding new output functions.

2.3 Acceleration method

The computation consumptions of $\Delta N(r)$ and $N(r)$ are in approximate proportion to the atom number in the simulation box and its second order, respectively. In general, the amount of atoms in a simulation box is larger than ten thousand. $\Delta N(r)$ needs to be recalculated after every MC simulation step for $\Delta U$ calculation, while $N(r)$ needs to be recalculated after every MC/EPMC equilibrium for $g(r)$ and $S(Q)$ update. Thus, these algorithms represent the highest consumption of the program’s calculation capacity. A graphics processing unit (GPU) has thousands parallel threads, so it is a suitable candidate to accelerate the calculation of $\Delta N(r)$ and $N(r)$. Fig.4 shows the CUDA C kernel function used to calculate $N(r)$ invoking GPU acceleration.

For the implementation of across nodes in a parallel configuration, mpich2 is used in the program[27]. Mpich2 is based on the MPI[13] standard and supports point-to-point and collective data communication among different nodes. Thus, it is highly efficiency in this program. For implementing a shared memory, multithread, parallel configuration within a computer or server node, Open Multi-Processing (OpenMP) is used in the C

++

program[6]. The OpenMP syntax supports the setting of a thread number dynamically, and a thread number cannot be known in advance in most cases. Thus, it is very convenient in programming. NeuDATool uses OpenMP and MPI API to distribute different $\Delta N(r)$/$N(r)$ to different GPU cards belong to a computer node or different nodes in a computer cluster. In addition, CPU serial algorithm version is edited for the programs wide applicability.

3 Performance Test

Neutron scattering data of full hydrogen, full deuterated, half deuterated H2O and amorphous silica (SiO2) samples at ambient temperature and pressure are used to test the performance of NeuDATool regarding the correctness and computational speed.

3.1 Correctness

Experimental neutron data of full hydrogen water, full deuterated water, half deuterated water, and SiO2 distributed with EPSR[28] and GudRun[29, 30] in the ISIS website are used to test the program.

For simulating the sample of water, we define a class to describe molecule of H2O by inheriting form molecule basic class. The molecule structure is maintained with harmonic oscillator potentials between every two atoms. Three different random movements i.e., H2O’s translation, H2O’s rotation, and atom’s (H or O) translation, are implement in the simulation. The EMP ($U^{EMP}$) of the atom type pairs of the sample are calculated with Eq.2.[20]

[TABLE]

Where $U^{EMP}_{j}(r)$ is the EMP of different the atom pairs, i.e., O-O, O-H, and H-H in water; $\Delta S_{i}(Q)$ is the different neutron structure factor between experiment and simulation of full hydrogen, full deuterated, and half deuterated H2O; $w_{ji}^{-1}$ is the inverse of weighted factor matrix in Eq.(1). In water samples, the number of neutron scattering profiles and atom type pairs are the same and equal to 3, So we employed matrix Invert() function provided by CLHEP[21, 22] to calculate the inverse matrix.

In the simulation of SiO2, we define two classes to describe Si and O as single atom molecules. The most difference from water is that there are only one neutron scattering profile but 3 different atom type pairs, i.e., Si-Si, Si-O, and O-O. As do in EPSR[20], we use a $3\times 3$ unit matrix to enlarge the $(1\times 3)$ weighted factor matrix $w_{ij}$ to a $(4\times 3)$ matrix $w_{ij}^{\prime}$, then use a Monte-Carlo method to calculate the pseudo-inverse matrix $w^{\prime-1}$.

Fig.5 shows the NeuDATool simulation results. They are consist with the experiments. It confirms that the distributions of small molecules observed in neutron scattering experiments have been reasonably represented by NeuDATool simulation.

To further verify the reliability of the simulation, we compared the PDF distribution of H2O from NeuDATool with other studies. Fig.6 shows the PDF distribution from NeuDATool simulation. Liquid water is a tetrahedrally random network. On average, 3.5 water molecules form hydrogen bonding with a center one. The results are consist with the studies by Kusalik, Head-Gordona and Soper[31, 32, 33].

All of those prove NeuDAToool can reconstruct the atomic structures of experimental samples correctly based on the neutron diffraction profile.

3.2 Computational Speed

A small computer cluster is used to test the speed performance of different methods. The cluster uses CentOS 7.3 as the operating system and has two nodes. Each node has two Intel Xeon Scalable Gold 6126 CPU (two Skylake–SP architectures, 12 cores, 24 threads 2.6 GHz, Turbo 3.7 GHz and a 19.25 MB L3 Intel smart cache), two Nvidia Tesla V100 calculation GPU card and 128 GB double data rate (DDR4) error correcting code (ECC) registered shared memory. The two nodes are connected with an InfiniBand (IB) connector (data transmission speed can increase up to 56 Gb/s).

A speed comparison of these acceleration methods is shown in Fig.7, and a more detailed quantitative comparison is listed in Tab.1. The water experiment data used in the speed test is the same as in the correctness test. As shown in the figure and table, mpich2 can improve the speed based on the ratio of nodes or thread numbers, while the GPU can provide an excellent acceleration ratio. Most importantly, with the GPU acceleration, the program can simulate a system comprising $>$ 1 million atoms. This is an essential improvement because it allows the program to simulate systems larger than 200 Å, so that it can analyze samples with macromolecules in all atomic models in the future.

4 Conclusions

The neutron scattering data analysis software NeuDATool is programmed with the object-oriented language C

++

. It makes the program flexible and friendly to users who need to define special molecules and MC random movement patterns. Potential functions and the corresponding parameters of the atomic force field can be modified or added by editing the C

++

head and source files. In addition, C

++

is an easy-to-read, high-level computer language, so users can try new algorithms and program flows to improve their analyses, and to calculate and output any important variables.

With the exception of parallel nodes within the server with OpenMP, parallel cross-different nodes of a computer cluster, and GPU hardware acceleration are supported. Specifically, with GPU acceleration, the calculation speed is improved considerably, so the program has the capacity to analyze disordered macromolecular samples and nanoparticles of all atomic models in the future.

Although the program is flexible for users and has a powerful calculation capacity, it was tested with a very limited number of control samples. Accordingly, in its current form, it is not a fully functional software package. The authors aspire to release it as an open-source toolkit framework for public use by interested scientists. In this sense, users will be able to contribute numerous new molecular classes, algorithms, and analyses routines in the future to make the program more powerful.

5 Acknowledgements

This work was supported by the National Key Research and Development Program of China [grant number 2017YFA-0403703]; and the National Natural Science Foundation of China (NSFC) [grant number U1830205, 21674020]. The authors thank the staff of the ISIS Disordered Materials Group for the explanations of the principles of EPSR and for providing experimental neutron diffraction data for testing of the program.