The original manuscript, published in the traditional manner, can be found in Microscopy and Microanalysis.

#Abstract

Crystalline materials used in technological applications are often complex assemblies composed of multiple phases and differently oriented grains. Robust identification of the phases and orientation relationships from these samples is crucial, but the information extracted from the diffraction condition probed by an electron beam is often incomplete. We therefore have developed an automated crystal orientation mapping (ACOM) procedure which uses a converged electron probe to collect diffraction patterns from multiple locations across a complex sample. We provide an algorithm to determine the orientation of each diffraction pattern based on a fast sparse correlation method. We test the speed and accuracy of our method by indexing diffraction patterns generated using both kinematical and dynamical simulations. We have also measured orientation maps from an experimental dataset consisting of a complex polycrystalline twisted helical AuAgPd nanowire. From these maps we identify twin planes between adjacent grains, which may be responsible for the twisted helical structure. All of our methods are made freely available as open source code, including tutorials which can be easily adapted to perform ACOM measurements on diffraction pattern datasets.

#Introduction

Polycrystalline materials are ubiquitous in technological applications. An ideal crystal structure can be fully defined with a small number of parameters: the 3 vectors defining its unit cell, and the position and species of each atom inside the unit cell (Borchardt 2011). To fully describe crystalline materials in the real world however, we require a description of both the crystal lattice, and all defects present in a given material. These include point defects such as dopants, vacancies, or interstitials (Dederichs, 1978), line defects such as dislocations (LeSar, 2014), planar defects including internal boundaries and surfaces (Tang, 2006), and volume defects such as precipitates (Kleiven, 2020). Strain fields in the surrounding material can be induced by each of these defects, or generated by the boundary or growth conditions of the material such as in thin film stresses (Janssen, 2007). One large subset of crystalline materials are polycrystalline phases, which consist of many small crystalline grains, arranged in either a random or organized fashion. Many material properties such as mechanical strength (Thompson, 2000), optical response (Park, 2019, Londono, 2021), or thermal or electrical conductivity (Castro-Méndez, 2019) are strongly modulated by the density and orientation of the boundaries between crystalline grains (Thompson, 1995). Thus characterizing the orientation of polycrystalline grains is essential to understanding these materials.

The two primary tools used to study the orientation of polycrystalline materials are electron backscatter diffraction (EBSD) in scanning electron microscopy (SEM), and transmission electron microscopy (TEM). EBSD can measure the orientation of crystalline grains with very high accuracy, but has limited resolution and is primarily sensitive to the surface of materials (Humphreys, 2001, Wright, 2011, Wright, 2015). Alternatively, we can directly measure the atomic-scale structure and therefore the orientation of polycrystalline grains, either by using plane wave imaging in TEM (Li, 2020), or by focusing the probe down to sub-atomic dimensions and scanning over the sample surface in scanning TEM (STEM) (Peter, 2018). This is possible due to the widespread deployment of aberration correction for both TEM and STEM instruments (Linck, 2016, Ramasse, 2017). Atomic resolution imaging, however, strictly limits the achievable field-of-view, and requires relatively thin samples, and thus is primarily suited for measuring polycrystalline grain orientations of 2D materials (Ophus, 2015, Qi., 2020).

Another approach to orientation mapping in TEM is to use diffraction space measurements. For crystalline materials, diffraction patterns will contain Bragg spots with spacing inversely proportional to the spacing of atomic planes which are approximately perpendicular to the beam direction (described by both the Laue condition and Bragg equations (Fultz, 2012)). To generate a spatially-resolved orientation map, we can focus a STEM probe down to dimensions of 0.5 to 50 nm, scan it over the sample surface, and record the diffraction pattern for each probe position, a technique referred to as nanobeam electron diffraction (NBED) (Ozdol, 2015), scanning electron nanobeam diffraction (SEND) (Tao, 2009]), or four dimensional-scanning transmission electron microscopy (4D-STEM) (we choose this nomenclature for this text) due to the 4D shape of the collected data (Bustillo, 2021). 4D-STEM experiments are increasingly enabled by fast direct electron detectors, as these cameras allow for much faster recording and much larger fields of view (Ophus, 2019, Nord, 2020, Paterson, 2020).

By performing template matching of diffraction pattern libraries on 4D-STEM datasets, we can map the orientation of all crystalline grains with sufficient diffraction signal. This method is usually named automated crystal orientation mapping (ACOM), and has been used by many authors in materials science studies (Zaefferer, 1994, Rauch, 2005, Wu, 2009, Kobler, 2013, MacLaren, 2020, Londoño-Calderon, 2020, Jeong, 2021, Zuo, 2021). ACOM experiments in 4D-STEM are highly flexible; two recent examples include (Lang, 2021) implementing ACOM measurements in liquid cell experiments, and (Wu, 2021) adapting the ACOM method to a scanning confocal electron diffraction (SCED) experimental configuration. ACOM is also routinely combined with precession electron diffraction, where the STEM beam is continually rotated around a cone incident onto the sample, in order to excite more diffraction spots and thus produce more interpretable diffraction patterns (Rauch, 2010, Brunetti, 2011, Moeck, 2011, Eggeman, 2015, Brunetti, 2011). Recently, (Mehta, 2020) have combined simulations with machine learning segmentation to map orientations of 2D materials, and (Yuan, 2021) have used machine learning methods to improve the resolution and sensitivity of orientation maps by training on simulated data. For more information, (Zaefferer, 2011) has provided a review of ACOM methods in SEM and TEM.

In this study, we introduce a new sparse correlation framework for fast calculation of orientation maps from 4D-STEM datasets. Our method is based on template matching of diffraction patterns along only the populated radial bands of a reference crystal's reciprocal lattice, and uses direct sampling of the first two Euler angles (which, in the convention we have adopted, correspond to the zone axis), and a fast Fourier transform correlation step to solve for the final Euler angle (in our convention, the in-plane rotation of the pattern). We test our method on both kinematical calculations, and simulated diffraction experiments incorporating dynamical diffraction. Finally, we generate orientation maps of polycrystalline AuAgPd helically twisted nanowires, and use clustering to segment the polycrystalline structure, and map the shared (111) twin planes of adjacent grains. Future parts of this study will improve our approach using machine learning methods (Munshi, 2021), and will extend our ACOM methods to multibeam electron diffraction experiments (Hong, 2021).

#Methods

#Overview

The problem we are solving is to identify the relative orientation between a given diffraction pattern measurement and a parent reference crystal. We solve this problem with three steps:

• First, we generate a diffraction pattern library which covers all unique crystal orientations using kinematical simulation. This library, stored in a sparse polar coordinate representation $P$, is called an ``orientation plan.''

• We find all diffracted spots/disks in each diffraction pattern, and convert them into the same sparse polar coordinate representation $X$.

• We determine the best fit orientation(s) by finding the maximum value(s) of the correlation $C$ between the diffraction patterns and the orientation plan.

All of the previously discussed ACOM implementations work in essentially the same way, i.e. by precomputing the diffraction library in some form, and then comparing each diffraction pattern to this library using a cost function based on some form of correlation. Performing template matching directly on diffraction patterns, which may contain millions of pixels, against a library of similarly sized patterns, is computationally expensive. However, the underlying information we are interested in, i.e.~the projected lattice in the pattern, is typically composed of at most a few dozen non-zero points. Our sparse correlation method involves reducing the diffraction patterns to a simpler representation where the correlation can be evaluated rapidly, by first detecting the positions of the Bragg disks in the pattern, then segmenting the data into radial bands, and only evaluating the correlation in the populated bands.

The primary advances of this paper are:

1. We use Fourier transforms along the annular direction in polar coordinates for both the diffraction library and diffraction patterns to efficiently solve for the in-plane image rotation. For a full polar coordinate transform, only a small number of radial bins will contain reciprocal lattice points, and thus the output is sparse along the radial direction. We utilize this sparsity by only evaluating the polar coordinate correlations on radial shells that contain reciprocal lattice points of the reference structure, making the calculations much faster.

2. We use Fourier transforms along the annular direction in polar coordinates for both the diffraction library and diffraction patterns to efficiently solve for the in-plane image rotation. For a full polar coordinate transform, only a small number of radial bins will contain reciprocal lattice points, and thus the output is sparse along the radial direction. We utilize this sparsity by only evaluating the polar coordinate correlations on radial shells that contain reciprocal lattice points of the reference structure, making the calculations much faster.

3. We give users fine-grained control over the relative weighting of diffraction peak radii and intensities in the correlation calculation, as well defining a kernel size which can be increased to allow more pattern distortion, or decreased to reduce the chance of false positive signals from grains with close orientations.

4. We provide all methods and codes as an open source implementation for the community to freely use and modify. Below we detail each of the steps for our orientation matching algorithm, and their required input calculations.

#Structure Factor Calculations

The structure factors of a given crystalline material are defined as the complex coefficients of the Fourier transform of an infinite crystal (Spence, 1993). We require these coefficients in order to simulate kinematical diffraction patterns, and thus we briefly outline their calculation procedure here.

First, we define the reference crystal structure. This structure consists of two components, the first being its unit cell defined by its lattice vectors $\vec{a}$, $\vec{b}$, and $\vec{c}$ composed of positions in $\vec{r}=(x,y,z)$, the 3D real space coordinate system. The second component of a crystal structure is an array with dimensions $[N, 4]$ containing the fractional atomic positions $\vec{p}_n = (p_{{a}}, p_{{b}}, p_{{c}})_n$ and atomic number $Z_n$, for the $n$th index of $N$ total atoms in the unit cell. Together these positions and atomic numbers are referred to as the atomic basis. Because $\vec{p}_n$ is given in terms of the lattice vectors, all fractional positions have values inside the range $[0, 1)$. The unit cell and real space Cartesian coordinates of the fcc Au structure are plotted in Figure 1a.

All subsequent calculations are performed in reciprocal space (also known as Fourier space or diffraction space). Thus the next step is to compute the reciprocal lattice vectors, defined by (Gibbs, 1884)

where $\times$ represents the vector cross product and $\Omega$ is the unit cell volume in real space. Note that this definition does not include factors of $2\pi$, and therefore all reciprocal coordinates have spatial frequency units.

Figure 1 - Crystal data. (a) Unit cell in real space for fcc Au. (b) Single atom scattering for Au atoms. (c) Structure factors of fcc Au in reciprocal space.

Next, we calculate the position of all reciprocal lattice points required for our kinematical diffraction calculation, given by

where $h$, $k$, and $l$ are integers representing the reciprocal lattice index points corresponding to the Miller indices $(h,k,l)$. We include only points where $|\vec{q}_{hkl}| < k_{\rm{max}}$, where $\vec{q}=(q_x,q_y,q_z)$ are the 3D coordinates in reciprocal space, i.e.\ those which fall inside a sphere given by the maximum scattering vector $k_{\rm{max}}$. To find all reciprocal lattice coordinates, we first determine the shortest vector given by linear combinations of $(\vec{a}^*,\vec{b}^*,\vec{c}^*)$, and divide $k_{\rm{max}}$ by this vector length to give the range for $(h,k,l)$. We then tile $(h,k,l)$ in both the positive and negative directions up to this value, and then remove all points with vector lengths larger than $k_{\rm{max}}$.

The reciprocal lattice defined above represents all possible coordinates where the structure factor coefficients $V_g(\vec{q})$ could be non-zero. The structure factor coefficients depend only the atomic basis and are given by

where $f_n$ are the the single-atom scattering factors for the $n$th atom, which describe the scattering amplitude for a single atom isolated in space. There are multiple ways to parameterize $f_n$, but here we have chosen to use the factors defined by (Lobato, 2014) which are implemented in py4DSTEM. Figure 1b shows the atomic scattering factor for an Au atom, while Figure 1c shows the structure factor amplitudes for fcc Au in 3D.

#Calculation of Kinematical Diffraction Patterns

Here we briefly review the theory of kinematical diffraction of finite crystals, following De Graef, 2003. We can fully describe an electron plane wave by its wavevector $\vec{k}$, which points in the direction of the electron beam and has a length given by $|\vec{k}| = 1/\lambda$, where $\lambda$ is the (relativistically-corrected) electron wavelength. Bragg diffraction of the electron wave along a direction $\vec{k}'$ occurs when electrons scatter from equally spaced planes in the crystal, described in reciprocal space as

For elastic scattering, $\vec{k}'$ has the same length as $\vec{k}$, and so scattering can only occur along the spherical surface known as the Ewald sphere construction (Ewald, 1921). For a perfect infinite crystal, scattering will seemingly almost never happen since it requires intersection of the Ewald sphere with the infinitesimally small points of the reciprocal lattice vectors. However, real samples have finite dimensions, and thus in reciprocal space their lattice points will be convolved by a \emph{shape factor} $D(\vec{q})$. Therefore diffraction can still occur, as long as Eq.~\ref{eq:diffraction_condition} is approximately satisfied.

If the sample foil is tilted an angle $\alpha$ away from the beam direction, the vector between a reciprocal lattice point $\vec{g}$ and its closest point on the Ewald sphere has a length equal to

The $s_{\vec{g}}$ term is known as the excitation error of a given reciprocal lattice point $\vec{g}$. When the excitation error $s_{\vec{g}}=0$, the Bragg condition is exactly satisfied. When the length of $s_{\vec{g}}$ is on the same scale as the extent of the shape factor, the Bragg condition is approximately satisfied.

A typical TEM sample can be approximately described as a slab or foil which is infinite in two dimensions, and with some thickness $t$ along the normal direction. The shape function of such a sample is equal to

Because this expression is convolved with each reciprocal lattice point, we can replace $q_z$ with the distance between the Ewald sphere and the reciprocal lattice point. For the orientation mapping application considered in this paper, we assume that $\alpha = 0$, and that the sample thickness $t$ is unknown. Instead, we replace Eq.~\ref{Eq:shape_foil} with the approximation

where $\sigma$ represents the excitation error tolerance for a given diffraction spot to be included. We chose this expression for the shape function because it decreases monotonically with increasing distance between the diffraction spot and the Ewald sphere $q_z$, and produces smooth output correlograms.

To calculate a kinematic diffraction pattern for a given orientation $\vec{w}$, we loop through all reciprocal lattice points and use Eq.\ref{eq:excitation error} to calculate the excitation errors. The intensity of each diffraction spot is given by the intensity of the structure factor $|F_{hkl}|^2$, reduced by a factor defined by either Eq.\ref{Eq:shape_foil} or Eq.\ref{Eq:shape_gaussian}. We define the position of the diffraction spots in the imaging plane by finding two vectors perpendicular to the beam direction, and projecting the diffraction vectors $q$ into this plane. The result is the intensity of each spot $I_m$, and its two spatial coordinates $(q_{m_x},q_{m_y})$, or alternatively their polar coordinates $q_m = \sqrt{{q_{m_x}}^2 + {q_{m_y}}^2}$ and $\gamma_m =$ arctan2$(q_{m_y}, q_{m_x})$. Note that the in-plane rotation angle is arbitrarily defined for kinematical calculations in the forward direction. The resulting diffraction patterns are defined by the list of $M$ Bragg peaks, each defined by a triplet $(q_{m_x},q_{m_y},I_m)$ in Cartesian or $(q_m,\gamma_m ,I_m)$ in polar coordinates.

Figure 2b shows diffraction patterns for fcc Au, along five different zone axes (orientation directions). Each pattern includes Bragg spots out to a maximum scattering angle of $k_{\rm{max}} = 1.5$ Å$^{-1}$, and each spot is labeled by the $(hkl)$ indices. The marker size shown for each spot scales with the amplitude of each spot's structure factor, decreased by Eq.\ref{Eq:shape_gaussian} using $\sigma = 0.02$ Å$^{-1}$.

Figure 2 - Diffraction pattern and orientation plan calculations. (a) Zone axes included in an orienation plan for fcc Au. (b) Diffraction patterns for all zones axes / projection directions included in the orientation plan in (a). (c) Slices of the orientation plan, where radial shells are shown in the vertical axis and the in-plane angle is shown on the horizontal axis.

#Generation of an Orientation Plan

The orientation of a crystal can be uniquely defined by a $[3 \times 3]$-size matrix $\overline{m}$, which rotates vectors $\vec{d}_0$ in the sample coordinate system to vectors $\vec{d}$ in the parent crystal coordinate system

where the first two columns of $\overline{m}$ given by $\vec{u}$ and $\vec{v}$ represent the orientation of the in-plane $x$ and $y$ axis directions of the parent crystal coordinate system, respectively, and the third column $\vec{w}$ defines the zone axis or out-plane-direction. The orientation matrix can be defined in many different ways, but we have chosen to use a $Z-X-Z$ Euler angle scheme (Rowenhorst, 2015), defined as

where $C_1 = \cos(\phi_1)$, $S_1 = \sin(\phi_1)$, $C_2 = \cos(\theta_2)$, $S_2 = \sin(\theta_2)$, $C_3 = \cos(\phi_3)$, and $S_3 = \sin(\phi_3)$. The Euler angles $(\phi_1, \theta_2, \phi_3)$ chosen are fairly arbitrarily, as are the signs of rotation matrices given above. %This convention was chosen to preserve internal consistency and to produce orientation matrices with sorted directions.

In order to determine the orientation $\overline{m}$ of a given diffraction pattern, we use a two-step procedure. The first step is to calculate an orientation plan $P((\phi_1, \theta_2), \phi_3, q_s)$ for a given reference crystal. The second step, which is defined in the following section, is to generate a correlogram from each reference crystal, from which we directly determine the correct orientation.

The first two Euler angles $\phi_1$ and $\theta_2$ represent points on the unit sphere which will become the zone axis of a given orientation. The first step in generating an orientation plan is to select 3 vectors delimiting the extrema of the unique, symmetry-reduced zone axes possible for a given crystal. Fig. Figure 2a shows these boundary vectors for fcc Au, which are given by the directions $[001]$, $[011]$, and $[111]$. We next choose a sampling rate or angular step size, and generate a grid of zone axes to test. We define a 2D grid of vectors on the unit sphere which span the boundary vectors by using the spherical linear interpolation (SLERP) formula defined by Shoemake, 1985. These points with a step size of $2^\circ$ are shown in Fig. \ref{Fig:ACOM_corr_match}d. The rotation matrices which transform the zone axis vector (along the z axis) are given by the matrix inverse of the first two terms in Eq. \ref{Eq:rotation_matrices}.

We then examine the vector lengths of all non-zero reciprocal lattice points $\vec{g}_{hkl}$ and find all unique spherical shell radii $q_s$. These radii will become the first dimension of our orientation correlogram, where each radius is assigned one index $s$. We loop through all included zone axes, and calculate a polar coordinate representation of the kinematical diffraction patterns. (Wu, 2009) pointed out that a polar transformation can make the in-plane rotation matching step more efficient, as the pattern rotation becomes a simple translation. We will further speed up the in-plane matching by using Fourier correlation along the angular dimension after the polar transformation (De Castro, 1987).

For each zone axis, the first step to compute the plan is to rotate all structure factor coordinates by the matrix inverse of the first two terms in Eq. \ref{Eq:rotation_matrices}. Next, we compute the excitation errors $s_g$ for all peaks assuming a $[0,0,1]$ projection direction, and the in-plane rotation angle of all peaks $\gamma_{\vec{q}}$. The intensity values of the orientation plan for all $q_s$ shells and in-plane rotation values $\phi_3$ are defined using the expression

where $\delta$ is the correlation kernel size, $\gamma$ and $\omega$ represent the power law scaling for the radial and peak amplitude terms respectively, $\rm{max}(...)$ is the maximum function, which returns the maximum of its two arguments, $\rm{mod}(...)$ is the modulo operator, and the summation includes only those peaks $\vec{g}$ which belong to a given radial value $q_s$. We have used the combined indexing notation for $(\phi_1, \theta_2)$ to indicate that in practice, this dimension of the correlation plan contains all zone axes, and thus the entire array has only 3 dimensions. The correlation kernel size $\delta$ defines the azimuthal extent of the correlation signal for each reciprocal lattice point. Note that Eqs. \ref{Eq:shape_gaussian} and \ref{Eq:shape_foil} are not used for the calculation of orientation plans.

We normalize each zone axis projection using the function

yielding the final normalized orientation plan

By default, we have weighted each term in the orientation plan with the prefactor $q_s |V_g|$, i.e.\ setting $\gamma=\omega=1$. The $q_s$ term gives slightly more weight to higher scattering angles, while the $|V_g|$ term is used to weight the correlation in favour of peaks with higher structure factor amplitudes, which was found to be more reliable than weighting the orientation plan by $|V_g|^2$, which weights each peak by its structure factor intensity.

Figure 2c shows 2D slices of the 3D orientation plan, for the diffraction patterns shown in Figure 2b. The in-plane rotational symmetry of each radial band is obvious for the low index zone axes, e.g. for the $[001]$ orientated crystal, the first row of the corresponding orientation plan consists of four spots which maintains the 4-fold symmetry of the diffraction pattern and can be indexed as $[020]$, $[200]$, $[0\overline{2}0]$ and $[\overline{2}00]$. The final step is to take the 1D Fourier transform along the $\phi_3$ axis in preparation for the Fourier correlation step defined in the next section.

#Correlation Pattern Matching

For each diffraction pattern measurement, we first measure the location and intensity of each Bragg disk by using the template matching procedure outlined by Savitzky, 2021. The result is a set of $M$ experimental diffraction peaks defined by the triplets $(q_m,\gamma_m ,I_m)$ in polar coordinates. Note that while all ACOM approaches we are aware of store the diffraction libraries in vector format (Rauch, 2005), here we also reduce the experimental data to a list of peak position and intensity vectors. This has the effect of deconvolving the probe shape from the diffracted disk, and thus improving the resolution. From the experimental peaks, we calculate the sparse polar diffraction image $X(\phi_3, q_s)$ using the expression

Note that the polar coordinates $q_s$ and $\phi_3$ used in this expression are identical to those used in the orientation plan calculation. \hl{The measured diffraction intensity is not normalized, as realistic sample thicknesses we expect the intensity to vary significantly from the kinematically predicted values.}

By default, we again use prefactors weighted by the peak radius and estimated peak amplitude given by the square root of the measured disk intensities. However, if the dataset being analyzed contains a large number of different sample thicknesses, multiple scattering can cause strong oscillations in the peak amplitude values. The intensity weighting factor $\omega$ provides a similar effect as the ``gamma correction'' used in many diffraction template matching routines (Cautaerts, 2021), but acts on the measured disk intensities rather than the original diffraction pattern. As we will see in the simulations below, in these situations the best results may be achieved by setting $\omega=0$, i.e.\ ignoring peak intensity and weighting only by the peak radii. Note that in the diffraction image, the correlation kernel size $\delta$ again gives the azimuthal extent of the correlation signal. However, in Eq. \ref{eq:image_polar} it also sets the range over which peaks are included in a given radial bin, and the fraction of the intensity assigned to each radial bin. To prevent experimental disk position errors from causing peaks to be assigned erroneously when the radial bins are near to one another (such as due to different reflections with nearly-similar spacing), experimental peaks can be included in multiple radial bins if they fall within the correlation kernel size of multiple bins. The kernel size $\delta$ can be optimized for each type of sample: if the sample contains crystals with large lattice distortions, a larger kernel size can be used to increase the tolerance. Alternatively, if a sample consists of many overlapping grains then the kernel size can be decreased to lower the probability of false positive matches for nearby orientations.

Finally, we calculate the correlation $C((\phi_1, \theta_2), \phi_3)$ of this image with the orientation plan using the expression

where $\mathcal{F}$ and $\mathcal{F}^{-1}$ are 1D forward and inverse fast Fourier transforms (FFTs) respectively along the $\phi_3$ direction, and the ${}^*$ operator represents taking the complex conjugate. We use this correlation over $\phi_3$ to efficiently calculate the in-plane rotation of the diffraction patterns. The maximum value in the correlogram will ideally correspond to the most probable orientation of the crystal. In order to account for mirror symmetry of the 2D diffraction patterns, we can also compute the correlation

where the mirror operation is accomplished by taking the complex conjugate of $\mathcal{F}\left\{X(\phi_3, q_s)\right\}$. For each zone axis $(\phi_1, \theta_2)$, we take the maximum value of $C$ and $C_{\rm{mirror}}$ in order to account for this symmetry. Figure 3b and c show output correlograms, for the corresponding diffraction patterns shown in Figure 3a. For each zone axis $(\phi_1, \theta_2)$, we have computed the maximum correlation value, which are plotted as a 2D array in Figure 3b. In each case, the highest value corresponds to the correct orientation.

Note that to calculate the correlation values, we have re-binned the vector peak data from both the orientation plan and experimental peaks into a polar coordinate image with sparse radial bins. It is also possible to perform the correlations of Eqs. \ref{eq:corr} and \ref{eq:corr_mirror} directly on the inputs into Eq. \ref{eq:orientation_plan} and experimental peaks $(q_m,\gamma_m ,I_m)$. However, in our numerical tests, correlations computed from vector inputs were slower than the image correlation approach for all ranges of parameters tested. We attribute this to two factors: first, the polar coordinate images we use have a very small number of radial bins since we only operate on shells which which contain reciprocal lattice vectors. Second, calculating the correlation of all in-plane rotations using Fourier transforms is highly efficient due to the high speed of the fast Fourier transform. This is why we have elected to compute the orientation correlations using radially-sparse polar coordinate images.

Figure 3b shows the correlation values along the $\phi_3$ axis, for the $(\phi_1, \theta_2)$ bins with the highest correlation value in Figure 3a. The symmetry of the correlation values in Fig. Figure 3c reflect the symmetry of the underlying patterns. For the $[0,0,1]$, $[0,1,1]$, and $[1,1,1]$, diffraction patterns, the in-plane angle $\phi_3$ correlation signals have 4-fold, 2-fold, and 6-fold rotational symmetry respectively. By contrast, the asymmetric diffraction patterns with zone axes $[1,1,3]$, and $[1,3,5]$ have only a single best in-plane orientation match. A diffraction pattern with the best fit match overlaid is shown in Figure 3.

The above default values are designed for matching of kinematical diffraction patterns. However, thermal excitation and multiple scattering can lead to non-zero intensities of the "kinematically forbidden peaks," i.e. diffraction signals at reciprocal lattice points where the structure factor is zero. In order to include forbidden peaks, we can include all points where $V=0$ in Eq.\ref{eq:structure_factor_coefs} by setting the structure factor threshold to zero. In a future update of the code, we will use dynamical (i.e. including multiple scattering) structure factor calculations to include peaks which are likely to be excited by multiple scattering. Additionally, we can set} $\omega=0$ in Eqs. \ref{eq:orientation_plan} and \ref{eq:polar_diff_image}, which removes the dependence of the correlation function on the peak intensities entirely, and uses only the peak positions. These steps will calculate the orientation correlation score using only the position of all scattering vectors, including the forbidden peaks.

Figure 3 - Diffraction pattern matching. (a) Projection direction correlation. (b) In-plane rotation correlation for best match shown in (a). (c) Diffraction pattern input spots shown as blue circles, spots corresponding to best-fit orientation match shown as black crosses.

#Matching of Overlapping Diffraction Patterns

In order to match multiple overlapping crystal signals, we have implemented an iterative detection process. First, we use the above algorithm to determine the best fit orientation for a given pattern. Next, the forward diffraction pattern is calculated for this orientation. We then loop through all experimental peaks, and any within a user-specified deletion radius are removed from the pattern. By default, this deletion radius is set to half of the correlation kernel size, i.e.~$0.5 \, \delta$. This value can be modified by the user depending on how close together the diffraction peaks are for a given experiment. Peaks which are outside of the deletion radius, but within the correlation kernel size, have their intensities reduced by a factor defined by the linear distance between the experimental and simulated peaks divided by the distance between the correlation kernel size and the deletion radius. Then, the ACOM correlation matching procedure is repeated until the desired number of matches have been found, or no further orientations are found. Note that while we could update the correlation score after peak deletion, we output the original magnitude of the full pattern correlogram in order to accurately calculate the probability of multiple matches.

#ACOM Integration into py4DSTEM

The ACOM pattern matching described has been implemented into the py4DSTEM python toolkit written by Savitzky, 2021. A typical ACOM workflow starts with using py4DSTEM to import the 4D dataset and one or more images of the vacuum probe. We then use a correlation template matching procedure to find the positions of all diffracted disks at each probe position \citep{pekin2017optimizing}. We use the correlation intensity of each detected peak as an estimate of the peak's intensity. The resulting set of $M$ peaks defined by the values $(q_m, \gamma_m, I_m)$ are stored as a PointList object in py4DSTEM. Because the number of peaks detected at each probe position can vary, we store the full set of all detected peaks in a \emph{PointListArray} object in py4DSTEM, which provides an interface to the ragged structured numpy data.

Most experimental datasets contain some degree of ellipticity, and the absolute pixel size must be calibrated. We perform these corrections on the set of measured diffraction disks using the py4DSTEM calibration routines defined by Savitzky, 2021. We know that the correlation approach is relatively robust against both ellipticity and small errors in the reciprocal space pixel size. However, precise phase mapping may require us to distinguish between crystals with similar lattice parameters; these experiments will require accurate calibration.

We perform ACOM in py4DSTEM by first creating a Crystal object, either by specifying the atomic basis directly, or by using the pymatgen package (Ong, 2013) to import structural data from crystallographic information files (CIF), or the Materials Project database (Jain, 2013). The Crystal object is used to calculate the structure factors, and generate an orientation plan. The final step is to use the orientation plan to determine the best match (or matches) for each probe position, from the list of calibrated diffraction peaks. If the sample contains multiple phases, we perform the orientation plan calculation and correlation matching for each unique crystal structure.

In addition to specifying the orientation plan spanning 3 vectors as in Figure 4a, we define additional methods to describe the space of possible orientations. One such example is fiber texture, where we assume the crystals are all orientated near a single zone axis known as the fiber axis, shown in Figure 4b. We can vary the angular range of zone axes included away from the fiber axis, as well as choose the azimuthal range around this axis in Figure 4b to account for symmetry around the fiber axis. Alternatively, an "automatic" option is provided, which uses pymatgen to determine the symmetry of the structure and automatically choose the span of symmetrically unique zone axes which should be included in the orientation plan, based on the point group (De Graef, 2003)symmetry. These structures and symmetries can be directly retrieved from the Materials Project database entries.

Figure 4 - Defining zone axes included in orientations plan. (a) Zone axis range defined by specifying three directions. (b) Zone axis range defined using fiber direction and angles.

#Simulations of Diffraction Patterns from Thick Samples

One important metric for the performance of an orientation mapping algorithm is how well it performs when the diffraction patterns contain significant amounts of multiple scattering. We have therefore used our ACOM algorithm to measure the orientation of simulated diffraction patterns from samples tilted along many directions, over a wide range of thicknesses. We performed these simulations using the multislice algorithm (Cowley, 1957), and methods defined by Kirkland, 2020 and Ophus, 2017. These methods are implemented in the Prismatic simulation code by Rangel DaCosta, 2021. The diffraction patterns were generated using a acceleration potential of 300 keV, a 0.5 mrad convergence angle, with real space and reciprocal pixel sizes of 0.05 Angstroms and 0.01 1/Angstroms respectively, with 4 frozen phonons. In total we have simulated 3750 diffraction patterns from Cu, Ag, and Au fcc crystals, over 25 zone axes ($[0,0,1]$ to $[3,4,4]$ excluding symmetrically redundant reflections) and thicknesses up to 100 nm with a 2 nm step size. These diffraction patterns were generated using the simulation pipeline and database defined by Munshi, 2021.

#Chemical Synthesis of Twisted AuAgPd Nanowires

The performed synthesis was reproduced with minor modifications from a known method given by (Wang, 2011). All reagents were purchased from Sigma Aldrich. We prepared the following solutions: 500 millimolar of Polyvinylpyrrolidone (PVP, MW 40,000) in Dimethylformamide (DMF), 50 millimolar Gold(III) chloride trihydrate (HAuCl$_{4} \cdot$ 3H$_2$O, $>$49.0% Au Basis) in DMF, \SI{50}{\milli\Molar} Silver nitrate (AgNO$_{3}$) in deionized (DI) water (resistivity $>$18 megaohm/cm), and 400 millimolar L-ascorbic acid ($>$99.0%, crystalline) in DI water. We created the reaction solution in a 4 mL glass vial (washed 3x with DI water and acetone) by mixing 800 microliter DMF, 100 microliter PVP, 20 microliter HAuCl$_{4}$, and 20 microliter AgNO$_{3}$. We mixed the solution for approximately 2 seconds using a Vortex-Genie 2 Mixer set to a value of 10, which spins the reaction solution at a speed of approximately 3200 rpm, then added 100 microliter of L-ascorbic acid solution drop-wise to the mixture while gently swirling by hand. At this point, the color changed from pale yellow to clear. We left the solution at room temperature for 7 days, at which point the solution was light brown/purple. The primary product of this reaction was straight, ultrathin Au-Ag nanowires (2 nm in diameter).

To twist the underlying ultrathin Au-Ag nanowires, we prepared solutions of 1.875 millimolar L-ascorbic acid and 2 millimolar H$_{2}$PdCl$_{2}$ in DI water. In a 4 mL glass vial (3x washed with DI water/acetone), we added 50 microliter of the Au-Ag reacted solution to 640 microliter of the L-ascorbic acid solution. Finally, we added 60 microliter of the H$_{2}$PdCl$_{4}$ solution and allowed the sample to incubate for at least 30 minutes. We purified the reaction solution by centrifuging the product down at 7500 rpm for 4 minutes. We decanted the supernatant, and then rinsed the reaction with DI water 3 times and re-dispersed in DI water. We prepared TEM samples of this material by depositing 10 microliter of purified nanowire solution onto 400 mesh formvar/ultrathin carbon grids.

#4D-STEM Experiments with Patterned Apertures

We collected the experimental data using a double aberration-corrected modified FEI Titan 80-300 microscope (the TEAM I instrument at the National Center for Electron Microscopy within Lawrence Berkeley National Laboratory). This microscope is equipped with a Gatan K3 detector and Continuum spectrometer, and was set to collect diffraction patterns integrated over 0.05 seconds, with 4x binning giving a calibrated pixel size of 0.00424 $Å^{-1}$. We used an accelerating voltage of 300 keV, an energy slit of 20 eV, and a spot size of 6. The beam current was 6 pA. We used a 10 micrometer bullseye aperture (probe size of approximately 1 nm) to form the STEM probe in order to improve detection precision of the Bragg disks (Zeltmann, 2020). We used a convergence semiangle of 2 mrad, with a camera length of 1.05 m. We recorded the experimental dataset using a probe step size of 5 Å, with a total of 286 and 124 steps in the x and y directions.

#Results and Discussion

#ACOM of Kinematical Calculated Diffraction Patterns

For the first test of our correlation method, we applied it to the same patterns calculated to generate an orientation plan for fcc Au. Next, we measured the calculation time and angular error between the measured and ground truth zone axes for each pattern. The results are plotted in Figure 5 for different maximum scattering angles $k_{\rm{max}}$ and angular sampling.

Figure 5 - Zone axis misorientation error. ACOM error of all zone axes as a function maximum scattering angle, and zone-axis and in-plane rotation step size.

The results in Fig.\ref{Fig:ACOM_kinematical} show that the angular error in zone axis orientation is relatively insensitive to the angular sampling. However, the angular error drops by a factor of 10 from approximately $\approx 3^\circ$ to $\approx 0.3^\circ$ when increasing the maximum scattering angle included from $k_{\rm{max}} = 1 \, Å^{-1}$ to $1.5 \, Å^{-1}$, and by another factor of 2-3 when increasing $k_{\rm{max}}$ to $2 \, Å^{-1}$. This is unsurprising, as examining Fig.\ref{Fig:ACOM_corr_match}e shows that there is a large number of visible Bragg spots outside of $k_{\rm{max}} = 1 \, Å^{-1}$, and because Bragg disks at higher scattering angles provide better angular precision relative to low $k$ disks. This result emphasizes the importance of recording as many diffraction orders as possible when performing orientation matching of 4D-STEM data. More spots can be included by collecting data out to higher scattering angles, or by reducing the convergence semiangle to bring weakly diffracting peaks above the noise floor. Setting $k_{\rm{max}}$ beyond the highest angle detected disks will not yield any additional precision but will make the orientation plan larger, so $k_{\rm{max}}$ should be chosen to correspond to the highest scattering angle peaks detected in an experiment.

The inset calculation times reported are for the single-threaded ACOM implementation in py4DSTEM, running in Anaconda on a laptop with an Intel Core i7-10875H processor, running at 2.30 GHz. The calculation times can be increased by an order of magnitude or more when running in parallel, or by using a GPU to perform the matrix multiplication and Fourier transform steps.

#ACOM of Overlapping Diffraction Patterns

A common feature of polycrystalline samples is overlapping grains along the beam direction, leading to overlapping diffraction patterns. To demonstrate the ability of our method to work with overlapping grains, we have generated a combined set of diffraction patterns with three low index zone axes and random in-plane rotations, plotted in Figure 6. Our ACOM code can often correctly return 3 zone axes which match the ground truth values. This example also demonstrates a procedure which could be used to map the location and orientation of multiple phases, even when the diffraction patterns overlap.

Figure 6 - ACOM of overlapping diffraction patterns. Note - there is a bug in the multiple pattern fitting code. Currently being fixed.

#ACOM of Dynamical Simulated Diffraction Patterns

In diffraction experiments using thick specimens, strong dynamical diffraction effects such as multiple scattering can occur. This effect is especially pronounced in diffraction experiments along low index zone axes, where the diffracted peak intensities oscillate as a function of thickness. In order to test the effect of oscillating peak intensities on our ACOM method, we have simulated diffraction patterns for Cu, Ag, and Au fcc crystals, along multiple zone axes. Some example diffraction patterns for the [011] zone axis of Au are plotted in Figure 7a. We see that all diffraction spots have intensities which oscillate multiple times as a function of thickness.

Figure 7 - Dynamical simulated diffraction patterns. (a) Example diffraction patterns for Au oriented to the $[011]$ zone axis for 10–80 nm thick slices. (b) Plots showing the mean zone axis misorientation in degrees as a function of thickness for Cu, Ag, and Au. Each plot shows the errors for correlation prefactors of $q_{{\rm s}} \vert V_{{\rm g}}\vert$ (red) and $q_{{\rm s}}$ (blue).

We performed ACOM by generating orientation plans with an angular sampling of $2^\circ$, a correlation kernel size of 0.08 $\rm{Å}^{-1}$, and maximum scattering angles of $k_{\rm{max}}$ = 1.0, 1.5, and 2.0 $\rm{Å}^{-1}$. We kept the radial prefactor of weighting set to $\gamma$ = 1, and tested peak amplitude prefactors of $\omega$ = 1.0, 0.5 and 0.0. The average zone axis angular misorientation as a function of thickness is plotted in Figure 7b. In total, we performed orientation matching on 3750 diffraction patterns, and a total of 33750 correlation matches on a workstation with an AMD Ryzen Threadripper 3960X CPU (2.2 GHz, baseclock). The typical number of patterns matched per second were of ~80-90, 45-55 and 25-30 patterns/s for $k_{\rm{max}}$ values of 1.0, 1.5, and 2.0 $\rm{Å}^{-1}$ respectively.

As expected, the errors are higher than those achieved under kinematic conditions, and the trend for smaller errors with larger $k_{\rm{max}}$ is also preserved (mean errors of 7.25$^\circ$, 3.09$^\circ$ and 1.39$^\circ$ for $k_{\rm{max}}$ values 1.0, 1.5 and 2.0 $\rm{Å}^{-1}$ respectively, $\gamma$ = 1, $\omega$ = 0.25). We did not observe any dependence of the orientation accuracy on the simulation thickness. Despite the correlation prefactor $|Vg|$ performing well for the examples shown in Figure 3, for the dynamical diffraction simulations along zone axes it was out-performed by prefactors of both $\sqrt{|V_g|}$ ($\omega$ = 0.5) and omitting the peak amplitude prefactor altogether ($\omega$ = 0). We therefore suggest that when mapping samples with a large range of thicknesses, or many crystals aligned to low index zone axes, the position of the diffracted peaks is significantly more important than their amplitudes or intensities. One possible method to increase the average accuracy for a randomly orientated sample while using higher amplitude prefactors is to perform an experiment which recovers more kinematical values for the diffracted peak intensities, for example by precessing the electron beam when recording diffraction patterns (Midgley, 2015 and Jeong, 2021).A precession experiment could however make the diffraction patterns of some grains more dynamical, and thus worsen the orientation accuracy for some probe positions. We note that there is likely no global optimal choice of orientation mapping hyperparameters for all materials and thicknesses, and this may be a worthwhile topic for future investigations.

#4D-STEM ACOM of Twisted AuAgPd Nanowires

We have tested our ACOM algorithm with a 4D-STEM dataset collected for an AuAgPd nanowires. These nanowires are morphologically twisted into double helices via a colloidal growth process as previously reported by \cite{wang2011}. A virtual dark field image nanowires is shown in Figure 8. For each probe position, we have recorded a full diffraction pattern, also shown Figure 8

we have calculated the maximum value across all STEM probe positions to generate a \emph{maximum diffraction pattern}, shown in Fig.~\ref{Fig:exp_structure}b. The beamstop used to block the center beam is visible, as well as various crystalline diffraction rings out to approximately 1.4 ${\rm{Å}}^{-1}$.

Figure 8 - Experimental 4D-STEM datacube.

Download links:

• Raw data (24.7 GB)

• Downsampled dataset (7.82 GB)

• Downsampled vacuum probe (18.1 MB)

• Tutorial dataset (1.95 GB)

• Tutorial vacuum probe (3.6 MB)

\begin{figure*}[htbp] \centering \includegraphics[width=6.4 in]{figure_AuAg_exp_v01.pdf} \caption{{\bf 4D-STEM scan of twisted polycrystalline AuAgPd nanowires.} (a) Diffraction image of probe over vacuum, showing bullseye pattern. (b) Maximum of each pixel in diffraction space over all probe positions. (c) Histogram of all peak locations detected by correlation in \pyFDSTEM{} of (a) with each pattern included in (b). (d) HAADF-STEM image of the sample. (e) 1D histogram of scattering vectors, with fcc AuAg inverse plane spacings overlaid.} \label{Fig:exp_structure} \end{figure*}

After performing the correlation peak finding algorithm in \pyFDSTEM, we have an estimated position and intensity of all detected Bragg peaks. A 2D histogram of these peaks, known as a Bragg vector map, is plotted in Fig.~\ref{Fig:exp_structure}c. Sharp polycrystalline diffraction rings are clearly visible, as well as false positives generated by the beamstop edge. These false positives were manually removed by using a mask generated from an image of the beamstop. A high angle annular dark field (HAADF) image was simultaneously recorded during the 4D-STEM data collection, which is shown in Fig.~\ref{Fig:exp_structure}d.

The final experimental pre-processing steps are to calibrate the diffraction pattern center, the elliptical distortions, and the absolute pixel size. We performed these steps by fitting an ellipse to the $(022)$ diffraction ring, and by assuming a lattice constant of $4.08 \, \rm{Å}$, corresponding to the fcc Au structure \citep{maeland1964lattice}. This process is explained in more detail by \cite{savitzky2021py4dstem}. We assumed that the Ag lattice constant is similar to that of Au. Despite the presence of Pd in the nanowires, there was no significant presence of secondary grains corresponding to the smaller lattice of fcc Pd grains. An intensity histogram of the corrected Bragg peak scattering angles are shown in Fig.~\ref{Fig:exp_structure}e. We have overlaid the 5 smallest scattering angles of Au on Fig.~\ref{Fig:exp_structure}e to show the accuracy of the correction.

\begin{figure*}[htbp] \centering \includegraphics[width=6.4in]{figure_AuAg_orientation_v03.pdf} \caption{{\bf Orientation mapping of polycrystalline AuAgPd nanowires.} (a) Total of measured correlation signal for each probe position. (b) Estimated number of patterns indexed for each probe position. (c) Example of 2 orientations indexed from a single diffraction pattern, collected at the position indicated by the arrow shown in (b), with correlation scores inset. (d) Orientation maps of the 3 highest correlation signals for each probe position. A legend for the crystallographic orientation is shown above, and arrows indicate the direction of the x and y axes, while the zone axis direction is out of the page.} \label{Fig:ACOM_AuAgPd} \end{figure*}

We have performed ACOM on the AuAgPd nanowire sample, with the results shown in Fig.~\ref{Fig:ACOM_AuAgPd} shown for up to 3 matches for each diffraction pattern. For each probe position, the sum of the maximum detected correlation signals for up to three matches are shown in Fig.~\ref{Fig:ACOM_AuAgPd}a. The structure is in good agreement with Fig.~\ref{Fig:exp_structure}, though with additional modulations due to some grains generating more diffraction signal than others. Using a correlation intensity threshold of 0.5, we have plotted the number of matching patterns in Fig.~\ref{Fig:ACOM_AuAgPd}b. The threshold of 0.5 was arbitrary chosen as a lower bound for a potential match, \hl{as the correlation values are scaled by the experimental intensity}. Examples of 2 matches to a single diffraction pattern are plotted in Fig.~\ref{Fig:ACOM_AuAgPd}c. In this figure the correlation score for the first matched pattern was higher than the second. The second match found shows some deformation between the measured and simulated Bragg peak positions, and matches fewer peaks. It therefore produces a lower correlation score, which can be used to threshold the results as in Fig.~\ref{Fig:ACOM_AuAgPd}d. Note that the threshold values for inclusion of any given match into the orientation maps is always user-defined.

Fig.~\ref{Fig:ACOM_AuAgPd}d shows the 3D orientations plotted as inverse pole figures for all probe positions, with the 3 best matches shown. Each image is masked by the total correlation signal, so that low correlation values are colored black. Almost every diffraction pattern with Bragg disks detected was indexed for at least one orientation with high confidence. Additionally, the patterns are very consistent, with a large number of adjacent probe positions recording the same orientation. Some secondary grains are also clearly visible in the second-best match, while very few patterns have been assigned a third match with high confidence.

\begin{figure}[htbp] \centering \includegraphics[width=3.4in]{figure_AuAg_twins_v08.pdf} \caption{{\bf Orientation analysis of grains in AuAgPd nanowires.} (a) Crystal grains, with in-plane (111) planes colored by orientation. (b) (111) planes shared by two overlapping grains.} \label{Fig:ACOM_111_twins} \end{figure}

In order to investigate the grain organization of the AuAgPd nanowires, we have performed clustering analysis on the orientation maps. Grains with similar orientations have been clustered together by looping through each probe position and comparing its orientation to its neighbors. Grains with at least 10 contiguous probe positions are shown in Fig.~\ref{Fig:ACOM_111_twins}a. $(111)$ planes which lie in the image plane are overlaid onto the grain strucure, colored by their orientation. Confirming our observations in Fig.~\ref{Fig:ACOM_AuAgPd}d, only a few grains with substantial overlap were reliably identified. This might be due to the low thickness of the sample (only a single grain along the beam direction), some grains not being oriented close enough to a zone axis to be detected, or multiple scattering deviations in the diffracted signal. There is a noticable bias in the orientation of the $(111)$ planes, which tend to be oriented horizontally near the growth direction of the nanowires.

One hypothesis for the growth mode of these twisted nanowires is that adjacent grains are connected by $(111)$ twin planes, forming local helical structures to give the observed twisted structures. To test this hypothesis, we determined the position of $(111)$ planes from Fig.~\ref{Fig:ACOM_111_twins}a which are shared by two overlapping grains. Fig.~\ref{Fig:ACOM_111_twins}b shows the location of these shared $(111)$ planes (with plane normal differences below $8^\circ$), colored by the normal vector of the plane. Many shared $(111)$ planes were detected, most with normal vectors aligned to the wire growth direction. These observations support the hypothesis that these nanowires are composed of grains connected by $(111)$ twin planes.

These experimental observations demonstrate the efficacy of our ACOM method. In order to improve these results, we will need to collect diffraction data with a wider angular range. This can be achieved by using precession electron diffraction \citep{rouviere2013improved}, multibeam electron diffraction \cite{hong2021multibeam}, or by tilting the sample or beam and recording multiple 4D-STEM datasets \citep{meng2016three}.

Note that this dataset has been binned by a large degree in order to preview in this manuscript. The actual Bragg peak detection and analysis was performed a higher resolution in diffraction space.

#Conclusion

We have introduced an efficient and accurate method to perform automated crystal orientation mapping, using a sparse correlation matching procedure. We have implemented our methods into the open source \pyFDSTEM{} toolkit, and demonstrated the accuracy of our method using simulated diffraction patterns, where we show that lowering or removing the peak-intensity weighting can improve the accuracy for thick samples with substantial dynamical diffraction. We also applied ACOM to an experimental scan of a complex helical polycrystalline nanowire, where we were able to identify shared twin planes between adjacent grains which may be responsible for the twisted helical geometry. All of our methods have been made freely available to the microscopy community as open source codes. We believe that our implementation of ACOM is efficient and accurate enough to be incorporated into automated online TEM software \citep{spurgeon2021towards}. In the future, we will improve our ACOM method using machine learning methods \citep{munshi2021ml}, and we will extend our ACOM methods to include multibeam electron diffraction experiments \citep{hong2021multibeam}.

#Acknowledgements

We thank Karen Bustillo for helpful discussions. CO acknowledges support of a US Department of Energy Early Career Research Award. SEZ was supported by the National Science Foundation under STROBE Grant no. DMR 1548924. AB, BHS, and py4DSTEM development are supported by the Toyota Research Institute. AR is supported by the 4D Data Distillery project, funded by the US Department of Energy. Work at the Molecular Foundry was supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231.