GAPD: a GPU-accelerated atom-based polychromatic diffraction simulation code

GAPD, a graphics-processing-unit (GPU)-accelerated atom-based polychromatic diffraction simulation code for direct, kinematics-based, simulations of X-ray/electron diffraction of large-scale atomic systems with mono-/polychromatic beams and arbitrary plane detector geometries, is presented. This code implements GPU parallel computation via both real- and reciprocal-space decompositions. With GAPD, direct simulations are performed of the reciprocal lattice node of ultralarge systems (∼5 billion atoms) and diffraction patterns of single-crystal and polycrystalline configurations with mono- and polychromatic X-ray beams (including synchrotron undulator sources), and validation, benchmark and application cases are presented.

1. Introduction

X-ray/electron diffraction is widely used to characterize microstructure at the lattice level. Diffraction simulations are useful in experimental design and interpretation (Bristowe & Sass, 1980; Budai et al., 1983; Derlet et al., 2005; Hawreliak et al., 2006; Brandstetter et al., 2008; Liu et al., 2014; Wang et al., 2015; E et al., 2015; Huang, 2010; Sun & Fezzaa, 2016). To correlate diffraction signatures with real structure characteristics such as defects, impurities, precipitates and finite grain size, direct atomic-based simulations of diffraction or reciprocal-space mapping are useful. The capability to simulate directly diffraction patterns of large atomic structures (billion atoms and more) with arbitrary configurations is highly desired. For example, size- and strain-induced diffraction peak broadenings of nanocrystalline solids under deformation are always intertwined (Gleiter, 1989, 2000; Revesz et al., 1996; Ungár, 2001; Budrovic et al., 2004; Barabash & Ice, 2014). In order to decouple these two effects, knowledge of the minimum size, above which size-induced broadening can be neglected, would be helpful but the exact value is still under debate.

Given the heavy amount of direct calculations from atomic configurations, the system size is quite limited. To improve computation speed, a method was proposed to implement fast Fourier transformation (FFT) with the assumption that electron density distribution is a sum of Gaussian profiles for each atomic position (Kimminau et al., 2008). One issue is related to the FFT sampling grid, in particular regarding choosing the parameters of Gaussian functions for different configurations. Besides possible artifacts induced by the Gaussian assumption, FFT in reciprocal space is limited to a triperiodic grid, inappropriate for simulations of curved surfaces (Ewald spheres) (Favre-Nicolin et al., 2011). Currently, direct diffraction simulations for systems consisting of millions of atoms are performed by the central processing unit (CPU) MPI (Message Passing Interface) parallel diffraction package (Coleman et al., 2013, 2014) integrated in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) (Plimpton, 1995), or with the stand-alone CPU MPI parallel code SLADS (Chen et al., 2017) and single-node graphics processing unit (GPU) CUDA (Compute Unified Device Architecture) parallel computations (Favre-Nicolin et al., 2011). For the CPU parallel code, calculations involving 10⁸ atoms and 10⁵ points in reciprocal space are still time-consuming (Chen et al., 2017).

For high and ultrahigh strain rate dynamic loading experiments, real-time in situ X-ray diffraction represents a highly promising development at synchrotron radiation (Turneaure et al., 2009; Luo et al., 2012; Fan et al., 2014; Hudspeth et al., 2015; Fan et al., 2016) and X-ray free-electron laser (Briggs et al., 2017) facilities. Unlike traditional quasi-static diffraction experiments performed with monochromatic X-ray beams (bandwidth ≃ 10^-4), higher-bandwidth `pink' beams (bandwidth > 10^-2) (Hudspeth et al., 2015; Fan et al., 2016; LCLS, 2017) are used for dynamic experiments, allowing more photons in a short exposure time (fs to µs).

Simulation codes for traditional polychromatic Laue diffraction have been developed to find the position and peak intensity of a Laue spot (Huang, 2010) and to correlate the spot shape with a given slip system (Tamura, 2014). Different from atom-based diffraction simulations, only the unit-cell structure is normally considered to represent a whole system in these codes. As a result, it is difficult to properly describe diffraction patterns of complex heterogeneous microstructures with varying grain sizes, shapes, strain gradients and phases as in real experiments.

Atom-based diffraction calculations for polychromatic beams require considerable computing resources compared with those for monochromatic beams since more wavelengths need to be considered. For a calculation of 10⁸ atoms and 10⁵ q-points, and a spectrum with 100 wavelengths, the computation time is 100 times that of monochromatic cases, which is largely unacceptable for a CPU cluster. On the other hand, the implementation of GPU may lead to tens to hundreds of times of acceleration from purely CPU-based calculations.

In this work, we present GAPD (https://www.pims.ac.cn/Resources.html), a GPU-accelerated atom-based polychromatic diffraction simulation code for direct, kinematics-based, simulations of X-ray/electron diffraction of large-scale atomic systems with mono-/polychromatic beams and arbitrary plane detector geometries. This code implements GPU parallel computation via both real- and reciprocal-space decompositions, and runs on a multi-node GPU cluster. With GAPD, we perform direct simulations of the reciprocal lattice node of ultra-large systems (∼5 billion atoms), obtain diffraction patterns of single-crystal and polycrystalline configurations with mono- and polychromatic X-ray beams, and present validation, benchmark and application cases.

2. Methodology

2.1. Reciprocal-space representation

Lattice points in reciprocal space as well as the diffraction geometry are shown in Fig. 1(a) for face-centered-cubic single-crystal Cu (blue); only a small portion is displayed for clarity. Here O is the origin of the reciprocal space, and A is the center of the Ewald sphere with a radius of , with λ being the wavelength of the incident beam. = is the incident wave vector, and = is the diffracted wave vector, and the angle between and is . B is a reciprocal lattice point that intersects the Ewald sphere, satisfying Bragg's law (Warren, 1969; Hammond, 2009),

where d_hkl is the crystal interplanar spacing. = is the reciprocal lattice vector or scattering vector, defined as

where and are unit vectors representing the scattered and incident beam directions, respectively. In addition,

An iso- circle on the Ewald sphere, centered at C, defines the azimuthal angles (γ). For lattice point B, ; = 0 () is defined by the user. Here the direction of is referred to as the transverse direction in GAPD.

Figure 1 Diffraction and detection geometries used by GAPD in (a) reciprocal space and (b) real space.

GAPD provides two modes for reciprocal-space construction. For the reciprocal space (RS) mode, one specifies the range and spacing of reciprocal lattice vectors along three reciprocal space axes, q_x, q_y and q_z. In the Ewald sphere (ES) mode, one specifies the direction of the incident beam, the transverse direction and the range and spacing of 2θ and γ, so only the reciprocal space points on an Ewald sphere are sampled.

2.2. Calculation of scattering intensity

The diffraction intensity I of N atoms at scattering vector is the product of structure factor with its complex conjugate, (Warren, 1969),

with

Here, is the position of the jth atom in real space. f_j is the atomic scattering factor to describe the scattering amplitude contributed by atom j at a scattering angle , and is parameterized as

for X-ray scattering and

for electron scattering. Parameters a, b and c have been tabulated for the majority of elements (Fox et al., 1989; Peng et al., 1996). Here the kinematical approximation is used, which assumes full coherent scattering (Vartanyants & Robinson, 2001). Different X-ray sources have different coherence properties. For example, an X-ray free-electron laser has significantly higher coherence than synchrotron radiation (Geloni et al., 2010; Vartanyants et al., 2011; Vartanyants & Singer, 2016), and there exists a coherence effect on diffraction (Vartanyants & Robinson, 2001). This effect will be included in a future version.

The Lorentz-polarization factor (Warren, 1969) for X-ray diffraction,

can be considered as an optional parameter in GAPD. Then the X-ray diffraction intensity follows as

Developed in the C++ language, the parallelization of GAPD is achieved by combining CUDA with MPI for its implementation on GPU high-performance computing clusters. We consider N_a atoms and scattering vectors in reciprocal space. Parallelization is realized in both real space (in terms of atoms) and reciprocal space (in terms of points).

Firstly, atoms are distributed over M MPI CPU threads, each of which contains N_p atoms, i.e.

For each CPU thread, structure factor calculations need to be performed. For each ,

Here, the CPU thread ranked l deals with the atoms labeled from to .

Then, the calculations of for the scattering vectors are distributed over CUDA GPU threads. The results are copied to CPU from GPU when a calculation is finished, and then saved into an array of dimension . When all the threads finish their own calculations, the MPI master thread sums individual in the array as

2.3. Projection from reciprocal space to 2D detector

For the sake of simulating diffraction patterns in synchrotron experiments, GAPD calculates the projection of diffraction wave vectors to a two-dimensional (2D) detector in an arbitrary position. The geometry in real space is shown in Fig. 1(b); xyz denotes the sample coordinate system or laboratory coordinate system (these two systems coincide), while d_xd_y denotes the detector plane. The coordinates of the scattered beam on the detector are expressed as

where P( p_x,p_y) is the point of normal incidence on the detector. and are unit vectors of the detector axes. is the projecting vector of the scattered beam in the sample coordinate system, and

where is parallel to the normal incident direction and is the sample-to-detector distance. The unit vector of the scattered beam can be obtained from , and λ with equation (2). In GAPD, dx is assumed to be always parallel to the xy-plane, and thus and can be derived from as

In this way, user-defined input includes (i) the coordinates of the point of normal incidence on the detector, (ii) the normal incident direction, (iii) the sample-to-detector distance, (iv) the beam incident direction, and (v) the wavelength of the probe beam. The first four items are defined in the sample system. Then, each point constructed in reciprocal space can be projected onto an arbitrarily positioned 2D detector. In GAPD, the beam incident direction and normal incident direction can be input as vectors in the sample coordinate system, or the angles with the xy-plane ( and ) and their orthogonal projections with the y-axis ( and ).

For a polychromatic beam, the intensity at a specific position on a 2D detector, , is the weighted integration over the incident beam wavelength range, [],

Here is the weight factor, e.g. the flux fraction of the incident beam. Each set of () corresponds to a scattering vector .

3. Validation and benchmark

GAPD is validated with electron/X-ray diffraction simulations of single-crystal Cu with various geometries. The detector is set to be perpendicular to the incident beams. The crystal coordinate system coincides with the sample coordinate system. We examine the following cases for the sake of validation: an electron beam (200 keV) with zone axes [001] and , 18.86 keV X-rays with zone axis [100], and 8.91 keV X-rays with zone axis [111]. Their corresponding 2D diffraction patterns are shown in Fig. 2, and are identical to standard indexed diffraction patterns and analytical predictions (Williams & Carter, 1996).

Figure 2 2D diffraction patterns calculated for Cu single crystals with X-rays or electrons along different zone axes: (a) 200 keV electrons, z = [001], (b) 200 keV electrons, z = , (c) 18.86 keV X-rays, z = [100] and (d) 8.91 keV X-rays, z = [111].

With the implementation of GPU acceleration, reciprocal-space mapping and 2D diffraction simulations of atomic systems with several billion atoms become realistic. We use GTX 980 GPUs, and each cluster node contains three GPUs. To evaluate the computing performance of GAPD, several tests with varying number of atoms (N_a), number of points () and number of GPU cluster nodes (N_nodes) are performed (Fig. 3).

Figure 3 Computing performance of GAPD: computation time/speed as a function of number of atoms, number of points and computing nodes as noted. Computation speed: per second.

As shown in Fig. 3(a), the computing time increases linearly as N_a increases at fixed N_nodes (9) and (1.3×10⁵). For instance, the computation time is about 5 h for a system of 5.4 billion atoms, which corresponds to a cube-shaped Cu single crystal with edge length of 400 nm. The computing efficiency increases with increasing [Fig. 3(b)]; the lower efficiency at small is likely to be due to insufficient use of GPU cores. The reciprocal and linear relations in Figs. 3(c) and 3(d), respectively, indicate satifactory parallelization efficiency of GAPD.

4. Application cases

The main features of GAPD, including reciprocal-space visualization, simulation of 2D diffraction patterns on an arbitrarily positioned detector, and considering polychromatic beams, are illustrated with the following three cases: the crystal size effect on reciprocal lattice nodes, the polychromaticity effect on single-crystal diffraction, and the polychromaticity effect on polycrystalline diffraction.

4.1. Crystal size effect on node-broadening in reciprocal space

It is well known that diffraction spot broadening becomes significant for small crystals. But an open question is how small is small? The lower bound above which crystal size-induced broadening can be neglected is still controversial, ranging from 100 nm (Warren, 1969) to 500 nm (Ungár, 2001). A direct simulation of such broadening with GAPD is instrumental for this matter.

A diffraction pattern measured in experiments is essentially a sampling of the Fourier or reciprocal-space representation of a specimen. Therefore, the node broadening in reciprocal space can be used to examine the size effect, without consideration of diffraction geometry and X-ray wavelength. A 3D node in reciprocal space is more appropriate to reflect the size effect than a 2D diffraction pattern, since the latter only samples a slice of the former.

Cube-shaped, defect-free, Cu single crystals with an edge length ranging from 4 nm to 400 nm are constructed, and examined in reciprocal space. The region around the () reciprocal lattice node of a 4 nm crystal is shown in Fig. 4(a), along with the relrods due to the small crystal size. In Fig. 4(b), the perimeter of the node is of a cube shape with round corners, which is more evident in the contour plot of the center cross-section of the node [Fig. 4(c)]. When an Ewald sphere intersects the node, the resulting intersection (a curved surface) depends on the diffraction geometry and X-ray wavelength, and may vary considerably. Consequently, one may obtain diffraction spots of different sizes and shapes, and thus different diffraction broadening for the same reflection plane, so diffraction patterns are less appropriate for such analysis.

Figure 4 (a) 3D visualization of reciprocal lattice node () enclosed in the red cube, along with the relrods due to the small crystal size (4 nm). (b) Enlarged view of the () node in (a). (c) Cross-sectional contour plot of a slice through the node center along the qxqy-plane.

For each crystal size, we obtain three reciprocal lattice nodes, (), () and (002). For each node, we obtain the intensity profile along the q_x-direction on the central q_xq_y-plane which cuts through the node center. The q_x profiles are fitted with a pseudo-Voigt function,

Here x and y refer to and the scattering intensity, respectively; y₀ is offset, m_u is a profile shape factor, x_c is the peak center, A is the peak area and ω is the full width at half-maximum (FWHM).

The FWHM of an intensity profile can be used to characterize node broadening or diffraction peak broadening. The fitted FWHMs are presented in Fig. 5 as a function of crystal size, which can be described with a reciprocal function. The FWHM versus crystal size curves are nearly identical for the three reciprocal-space nodes. The FWHM decreases rapidly with increasing crystal size and approaches zero asymptotically. For a moderately high reciprocal resolution in experiments, ≃ 10^-3 (Lienert et al., 2017). The corresponding values are 4.8×10^-4, 5.5×10^-4 and 9.1×10^-4 Å⁻¹ for (), (002) and (), respectively. For FWHM values below the resolution , the peak broadening is deemed negligible. Then, the minimum crystal size beyond which peak broadening can be neglected ranges from 110 to 208 nm for the three reciprocal lattice nodes investigated. For most metals, q of diffraction index ranges from 0.2 to 1.65 Å⁻¹. Then minimum crystal sizes can be deduced from the reciprocal relation between FWHM and crystal size, ranging from 60 nm to 500 nm.

Figure 5 FWHM of intensity profiles of reciprocal lattice nodes versus crystal size for three reciprocal-space nodes.

4.2. Diffraction of small single crystals with polychromatic X-rays

To evaluate the effect of polychromaticity on the diffraction of small single crystals, we simulate diffraction patterns of Cu single crystals with four different sizes: 4 nm, 14 nm, 36 nm and 100 nm. X-rays are incident perpendicularly on the (112) plane. We choose the first harmonic of undulator U18G13 (period = 18 mm, gap = 13 mm) at the Advanced Photon Source (APS) 32-ID beamline (Fan et al., 2016) as the polychromatic X-ray source [Fig. 6(a)]; its spectral flux peaks at 24.65 keV or = 0.5029 Å. As a reference, we also simulate the diffraction for monochromatic X-rays with a wavelength of = 0.5366 Å satisfying Bragg's law in the same geometry.  is indicated by the red dot in Fig. 6(a). The simulated diffraction patterns are plotted as the γ–2θ plot in Fig. 6(b), and diffraction curves in Fig. 6(c), for diffraction spot ().

Figure 6 (a) X-ray spectrum of the first harmonic of the U18 undulator with a gap of 13 mm at the APS beamline 32-ID (APS U18G13). (b) The γ–2θ plots of diffraction spot from 4 nm single-crystal Cu for monochromatic X-rays ( = 0.5366 Å), and (c) for polychromatic X-rays with APS U18G13. (d) Diffraction profiles along through the center of the diffraction spot of single Cu obtained with the polychromatic source, for various crystal sizes. The black arrow denotes the Bragg peak position.

For the 4 nm single crystal, diffraction spot () is circular on the γ–2θ plot in the case of the monochromatic X-rays [Fig. 6(b)]. However, the corresponding diffraction spot for the polychromatic source is elliptical, elongated along the 2θ direction [Fig. 6(c)], and its center shifts by 0.6° toward lower angle.

The intensity versus profiles through the center of diffraction spot () formed with the polychromatic source are presented in Fig. 6(d) for various crystal sizes. There is an increasing broadening and peak shift toward lower with decreasing grain size. Essentially, the broadening in the direction is induced by the crystal size effect discussed in §4.1, while multiple wavelengths do augment the broadening.

For a monochromatic beam, the Ewald sphere can only sample a small portion of a reciprocal lattice node, while the node intersects with multiple Ewald spheres for a polychromatic beam. A broader bandwidth indicates a larger portion to be sampled by Ewald spheres of different radii, and, thus, more pronounced broadening. As shown above, the size of a reciprocal lattice node increases with decreasing crystal size. The size and intensity of a diffraction spot is the combined result of crystal size (node dimensions) and multiple wavelengths (polychromaticity; multiple Ewald spheres).

The Ewald sphere corresponding to the Bragg wavelength intersects the center of a reciprocal lattice node with the highest intensity for small crystals. Since the intensity at a pixel on the detector is an integration of contributions from different localities of a node weighted by the flux amplitudes of corresponding wavelengths, non-Bragg wavelengths with higher fluxes may still lead to the highest diffraction intensity, and thus a peak shift.

4.3. Diffraction of nanocrystalline solids with polychromatic X-rays

A nanocrystalline Cu structure consisting of 250 randomly oriented grains with a mean grain size of 8 nm is examined. Figs. 7(a)–7(d) present 2D diffraction patterns obtained with monochromatic X-rays (wavelength = 0.5029 Å, 0% bandwidth), polychromatic X-rays with a Gaussian-shaped spectrum (centered at , 4% bandwidth), the APS U18G13 undulator source with single harmonic shown in Fig. 6(a) (spectral flux peak at , ∼8% bandwidth) and the APS U33G25 undulator source with multiple harmonics, respectively. The X-ray spectrum for APS U33G25 is shown in Fig. 7(e). With increasing bandwidth, the diffraction spots on Debye–Scherrer rings become stretched along the direction. The broadening can be simply explained by Bragg's law [equation (1)]. For a polycrystalline solid with random grain orientations and sufficient number of grains, the incident angle of X-rays relative to the same group of crystal planes, θ, can be arbitrary. Thus, for fixed d_hkl, the number of θ angles that satisfy Bragg's law increases with increasing number of wavelengths or bandwidth.

Figure 7 2D diffraction patterns of a 250 grain polycrystalline system for (a) 0.5029 Å monochromatic X-rays, (b) polychromatic X-rays with a Gaussian-shaped spectrum, (c) the first harmonic of APS undulator source U18G13, and (d) the first three harmonics of APS undulator source U33G25. (e) X-ray spectrum of APS U33G25. (f) Corresponding diffraction curves.

We integrate azimuthally the 2D diffraction patterns, and obtain 1D diffraction curves in Fig. 7(f), each normalized by its maximum intensity. Both diffraction curves from the Gaussian and APS U18G13 spectra are broadened along the direction, and the peak shapes are similar to their corresponding X-ray spectra. However, the {111} and {200} peaks shift by 0.14–0.26° toward higher for the APS U18G13 spectrum, while no peak shift is seen for the Gaussian spectrum. The difference is caused by the asymmetry or symmetry of their respective spectra.

The finite grain size leads to diffraction peak broadening in the direction for each wavelength in a polychromatic beam, since the broadened nodes in reciprocal-space representation of the nanocrystalline Cu allow their intersection with Ewald spheres of different radii. As a result, summing those broad­ened peaks (for each wavelength) over the asymmetric APS U18G13 spectrum does not necessarily lead to the maximum intensity at , where the spectral flux peaks.

For multiple harmonics [APS U33G25; the green curve in Fig. 7(f)], both the {111} reflection from the second harmonic and the {220} reflection from the third harmonic contribute to the diffraction peak. The `plateau' following the peak is due to the {200} reflection from the second harmonic and {311} from the third harmonic. Thus, the overlap of diffraction intensities from different crystal planes and different harmonics renders it difficult to analyze diffraction peaks with conventional methods. At present, simulating diffraction patterns with forward simulation codes such as GAPD in comparison with experiments is useful for interpreting multiple-harmonic diffraction data of poly/nanocrystalline solids.

5. Conclusions

We present a GPU-accelerated parallel simulation code, GAPD, for simulating electron/X-ray diffraction with mono-/polychromatic beams directly from atomic configurations. Diffraction simulation on super-large systems (∼5 billion atoms) is demonstrated, and the system size can be scaled up by a factor of 10–100 on more powerful clusters. GAPD is utilized to explore the crystal size effect on node-broadening in reciprocal space, and the influence of polychromaticity on peak-broadening of single-crystal and polycrystalline nanomaterials.

In particular, (i) for a moderately high reciprocal resolution ( ≃ 10^-3), peak-/node-broadening can be neglected for crystal size above 500 nm. Precise minimum crystal sizes depend on (hkl) and q-resolution.

(ii) For small single crystals, polychromatic beams with asymmetric spectra induce peak shift and broadening in the -direction, which diminish with increasing crystal size.

(iii) For polycrystalline solids, a diffraction peak is broaden­ed by a polychromatic beam, and its shape follows that of the beam spectrum. Asymmetric spectra induce both peak broadening and shift.