Atomic resolution coherent x-ray imaging with physics-based phase retrieval

Abstract

Coherent x-ray imaging and scattering from accelerator based sources such as synchrotrons continue to impact biology, medicine, technology, and materials science. Many synchrotrons around the world are currently undergoing major upgrades to increase their available coherent x-ray flux by approximately two orders of magnitude. The improvement of synchrotrons may enable imaging of materials in operando at the atomic scale which may revolutionize battery and catalysis technologies. Current algorithms used for phase retrieval in coherent x-ray imaging are based on the projection onto sets method. These traditional iterative phase retrieval methods will become more computationally expensive as they push towards atomic resolution and may struggle to converge. Additionally, these methods do not incorporate physical information that may additionally constrain the solution. In this work, we present an algorithm which incorporates molecular dynamics into Bragg coherent diffraction imaging (BCDI). This algorithm, which we call PRAMMol (Phase Retrieval with Atomic Modeling and Molecular Dynamics) combines statistical techniques with molecular dynamics to solve the phase retrieval problem. We present several examples where our algorithm is applied to simulated coherent diffraction from 3D crystals and show convergence to the correct solution at the atomic scale.

Introduction

Coherent diffraction imaging (CDI) is a lens-less imaging technique that uses diffracted x-rays, electrons, or visible light to image a sample^1,2,3. To successfully image a sample, CDI must obtain both the phase and amplitude of the measured photons. During measurement, however, only the signal is recovered and the phase information of the photons is lost^4,5,6. This problem, called the phase problem, is not unique to CDI, and occurs in several other areas of physics, imaging, and signal processing^7,8.

A particular application of CDI that has potential for near-atomic resolution is Bragg CDI (BCDI)^9,10,11. In BCDI, the coherent x-rays scattered from small micron-sized crystals are measured in a three dimensional region around one or more Bragg peaks in reciprocal space. An iterative phase retrieval algorithm can then reconstruct the 3D strain field of the micro-crystal or metallic grain from the measured diffraction intensities^{12,13,14,15,16}. Furthermore, recent advances in BCDI have enabled measurements of multiple Bragg peaks which allows reconstruction of the full strain tensor inside the crystal^14,16,17,18. Much of this work has been enabled by complimentary Laue diffraction microscopy that enables the nanocrystal or sample in question to have its crystalline orientation determined and multiple Bragg peaks easily measured^19,20.

Most state-of-the-art phase retrieval algorithms use an iterative combination of several methods to solve the phase problem^2,21,22. The most fundamental of these methods is an alternating projection method called error reduction (ER)^22,23,24. Because ER is equivalent to gradient descent with a sub-optimal step size²² and the phase problem is naturally non-convex²², ER alone tends to stagnate in local minima. Several global optimization methods have been developed that allow phase retrieval algorithms to escape local minima, such as hybrid input/output (HIO), charge flipping, and shrinkwrap^25,26,27. These phase retrieval methods have been used successfully in recent years to study materials in operando, yielding insights into failure mechanisms of battery materials^{28,29,30,31,32} and imaging strain in platinum nanoparticles during reduction of carbon monoxide for automotive exhaust reduction^33,34,35.

While many successful phase retrieval algorithms exist, few phase retrieval algorithms incorporate materials models. Materials models, such as density functional theory, molecular dynamics (MD), mesoscale modeling, and finite element models, contain a large amount of information about the materials they imitate. If implemented correctly, adding information from materials models to phase retrieval algorithms would constrain the solution to be physically reasonable and significantly restrict the search space of the phase problem. Recent work has indirectly incorporated materials models through physics-aware algorithms based on machine learning principles^36,37. These physics-aware algorithms train machine learning models on data obtained from MD simulations and their corresponding simulated diffraction patterns. However, these algorithms do not explicitly incorporate a physical model during the solving process.

Here, we present a phase retrieval algorithm, called PRAMMol for Phase Retrieval with Atomic Modeling and Molecular Dynamics (see Fig. 1). PRAMMol directly incorporates a physical model during the solving process and was designed with expected diffraction limited storage rings (DLSR) synchrotron upgrades in mind. These synchrotron upgrades will allow for an increase of up to 200 times brightness in coherent x-rays, and may enable sub-Ångstrom resolution^10,11. While atomic resolution has routinely been demonstrated with coherent electron sources^38,39, atomic resolution from CDI has been more challenging to achieve. Because upgraded synchrotrons could enable atomic scale in operando BCDI, PRAMMol was designed to work at the atomic scale and uses MD as its physical model (All LAMMPS^40,41,42 calculations are run with an Al potential⁴². However, the methodology presented here is agnostic to the specific potential used and any monoatomic interatomic potential should achieve similar results.).

Fig. 1: Proposed PRAMMol method.

(left) A coherent x-ray beam in focused on a small nanocrystal or grain of metal, and the intensity around several Bragg peaks (3 in this work) is measured, (right) PRAMMol reconstructs the atom positions using an iterative procedure consisting of three steps, maximum likelihood estimation, cross validation, and molecular dynamics energy minimization. The Bragg peak measurements are incorporated into the algorithm during several different steps, but mainly the Maximum Likelihood Estimation (MLE) optimization and Cross Validation as described in the “PRAMMol Algorithm”.

Results

Time integrated photon flux trials

To demonstrate the power of physics-based phase retrieval and estimate the achieved resolution of PRAMMol across a large range of available photon fluxes, 10 trials were run across 20 values of time integrated photon flux (TIPF)¹⁰ for a total of 200 trials. As described in “Grain Creation”, a 4679 atom face centered cubic (FCC) crystal of ~6 nm diameter was created with a single vacancy at its center (see Fig. 2). For each trial, three diffraction patterns were simulated at the (−1,1,1),(1,−1,1) and (1,1,−1) peaks (see “Diffraction Pattern Simulation” and Fig. 3). While this choice of peaks is arbitrary, three linearly independent peaks are required to successfully reconstruct any three dimensional sample. To decrease the total computational cost, the initial guess of each reconstruction was a low-resolution reconstruction which can be determined through ER/HIO methods. For this work, we accomplished this by creating a coarse grain model of the true object with a resolution of 5 Å. This low-resolution solution was filled with atoms with an FCC lattice with a lattice constant of 4 Å, which on purpose was initiated to not be the exact correct lattice constant. These atoms were placed randomly (i.e., with respect to the true atom positions) but with the correct orientation. The relative lattice angular orientation and lattice spacing can be ascertained from the locations of the Bragg peaks in reciprocal space.

Fig. 2: A simulated FCC crystalline sample of 4679 atoms was created with atomsk with a single vacancy at the center.

a Shows the full 3D grain. b, c, d are slices through the center of the object with the vacancy visible in each of them.

Fig. 3: Three peaks are simulated from the 6 nm sample with a single vacancy at its center.

a Shows the locations of these peaks in reciprocal space. b shows an H-K slice through the $\left(\bar{1},1,1\right)$ peak at the slice with the brightest point in the diffraction pattern. c shows this same slice with Poisson noise added to the diffraction pattern.

After all 200 trials were complete, the resolution was calculated as described in “Resolution Calculation”. Figure 4 shows that PRAMMol is able to achieve atomic reolution at a TIPF of about 10¹⁸ ph μm^−2. Additionally, once atomic resolution is achieved, PRAMMol reconstructs at nearly picometer resolution and continues to improve resolution as TIPF increases.

Fig. 4: To evaluate the PRAMMol method, 10 trials were run across 20 values of TIPF, for a total of 200 trials.

These studies were conducted on the grain with a void from Fig. 2. The atomic position error and resolution are shown as a function of TIPF. PRAMMol achieves atomic resolution at a TIPF of ~10¹⁸ ph μm^−2. This results in resolution and mean average error (MAE) of much less than 1 pm and maximum error (ME) of ~1 pm after atomic resolution is achieved.

Material defect trials

This analysis required the lattice information of the crystal to be known in advance. While this information can be easily determined from the measured diffraction pattern, we also show that PRAMMol is able to successfully reconstruct samples without the lattice information. As described in “Grain Creation”, three FCC crystal samples were created: a 4791 atom sample with no defects, a 4679 atom sample with a single vacancy, and an 8777 atom sample with a screw dislocation. For each sample, a diffraction pattern was simulated with a TIPF of ~10²⁰ ph μm^−2 (see “Diffraction Pattern Simulation” and Fig. 3).

PRAMMol started the reconstruction with a random number of atoms and random atom positions. PRAMMol then reconstructed the atom positions of the samples from their simulated diffraction patterns (modulus squared) and a low-resolution reconstruction of the sample. For the crystal with a single vacancy, PRAMMol converged to the correct atomic positions in ~30 iterations. PRAMMol was also able to find the correct number of atoms in the original object and an unknown scaling constant that encompasses the TIPF and several other experimental factors (Fig. 5). For any single atom in the reconstruction the position error was less than 1 pm (see Fig. 6).

Fig. 5: Reconstruction of sample with a vacancy takes a little more than 30 iterations to converge to the correct scaling constant and number of atoms.

After this, the algorithm iterates to minimize error of specific atom positions. During these last few iterations, only error minimization of the diffraction pattern is occurring rather than resampling atom positions.

Fig. 6: Results of PRAMMol reconstructions of disjoint diffraction patterns with atomic resolution.

PRAMMol effectively reconstructs the grains to extremely high fidelity. Shown here are the X-Y slice (left column), the Y-Z slice (center column), and the X-Z slice (right column) of three different crystals: bulk with no defects (top row), a single vacancy (middle row), and a screw dislocation (bottom row). Also shown is that the atomic positions were retrieved with errors less than 0.01 Å. This shows that the accuracy of PRAMMol relies more on the accuracy of the interatomic potential and noise in the diffraction pattern rather than pixel size or aliasing effects that traditional phase retrieval encounters.

Similarly, for the defect-free sample and the screw dislocation, reconstructions resulted in position errors of less than 1 pm for any single atom. The number of iterations to convergence varied based on how complex the samples were (see Fig. 6 and Supplementary Videos). The reconstruction of the bulk crystal converged to the correct object in around 40 iterations, while the sample with the screw dislocation required 200 iterations to converge.

Temperature trials

In these simulations, x-rays were assumed to scatter off of point particles at a well defined position. For the grains shown in Figs. 2–6, this is a valid assumption because the grains were relaxed and had no temperature dependence (i.e., at 0 K). Experimentally, however, picometer resolution is very unlikely to be achieved for several reasons, one of which is lattice vibrations due to finite temperature. To test the effectiveness of PRAMMol in more realistic experimental conditions, we completed a temperature study with 10 independent reconstruction trials across 10 temperatures for a total of 100 trials. More details can be found in “Temperature-Dependent Diffraction Pattern”. As described in “Grain Creation”, an 8777 atom sample with no defects was created. For each trial, a diffraction pattern with a TIPF of ~10²⁰ ph μm^−2 was created for the sample at the desired temperature as described in “Temperature-Dependent Diffraction Pattern”. Again, the initial guess of each reconstruction was a low-resolution approximation filled with FCC atoms with a lattice constant of 4 Å.

As shown in Fig. 7, PRAMMol was able to reconstruct at sub-Angstrom resolution at or below 400 K. Additionally, the mean average error of the positions of these reconstructions was beneath 0.1 Å. From Fig. 8, we see that PRAMMol is able to successfully reconstruct the interior of the sample, but has some difficulty reconstructing the surface at room temperature. This, however, is somewhat unsurprising because of the unphysical nature of the simulation used to create the diffraction pattern. Because the object is in free space and energy is being continuously added and removed to preserve the correct temperature, there is no guarantee that a nonzero total angular or linear momentum is not introduced. While the object was centered at every iteration, negating the linear momentum, the total angular momentum was allowed to fluctuate, causing the object to rotate more than a sample fixed to a substrate would.

Fig. 7: To assess the PRAMMol in experimental conditions, 10 trials were run accross 10 temperatures, for a total of 100 trials.

For temperatures up to 500 K, PRAMMol effectively reconstructs the object at sub-Angstrom resolution. Above 700 K, PRAMMol loses the ability to reconstruct at near-atomic resolution.

Fig. 8: At 100 and 300 K, PRAMMol effectively reconstructs the entire object at atomic resolution except for a few locations on the surface.

At 700 K, PRAMMol is still able to reconstruct the original object to high fidelity, but occasionally loses atoms. At 1000 K, the object resembles a low-resolution reconstruction that traditional phase retrieval algorithms would output.

Discussion

While atomic resolution BCDI has yet to be demonstrated experimentally, we have shown through simulation that PRAMMol can achieve atomic resolution on small particles at values of TIPF that will be theoretically possible with synchrotron upgrades. Furthermore, we have shown that the only requirements to achieve atomic resolution besides the necessary photon flux is a low-resolution reconstruction found through traditional BCDI methods and an accurate interatomic potential.

Additionally, we have shown that PRAMMol is still able to achieve atomic resolution at ambient temperatures. This temperature study has shown that PRAMMol is proficient well above room temperature. Additionally, rather than lose fidelity completely at temperatures near melting, PRAMMol effectively returns a reconstruction that resembles the output of traditional phase-retrieval with nanometer resolution.

These results are not exhaustive, but they do give compelling evidence to investigate BCDI methods that directly incorporate a physical model. We have shown that PRAMMol, paired with a low-resolution reconstruction found through traditional BCDI methods, has the potential to achieve atomic resolution with less photon flux than traditional BCDI methods alone and that PRAMMol continues to outperform these methods as photon flux increases. Additionally, PRAMMol has several properties that allow it to perform better than projection based methods in certain situations and could extend BCDI to atomic resolution without the need to measure the small angle x-ray scattering signal.

Traditional reconstructions require a regular grid, not allowing for electron density to be defined between them. PRAMMol, on the other hand, defines atomic positions and allows these to be placed anywhere in space, removing any pixel resolution limit that is inherent in traditional BCDI algorithms. Additionally, because atom positions are well-defined point particles, it is straightforward to solve using multiple, disjoint peaks during the reconstruction. While current multi-peak BCDI algorithms are currently used^14,16,17, these methods are not designed for atomic resolution. Also, these algorithms must preprocess the measured data to ensure that the Bragg peaks are centered and that the frame of reference of the data is the same as the frame of reference of the object, introducing errors in the reconstruction.

However, because PRAMMol calculates the discrete Fourier transform explicitly, the computational cost increases to ${{{\mathcal{O}}}}({n}^{2})$. Because of this restriction, the current implementation would have difficulty scaling to larger objects with more than a few hundred thousand atoms. A focus of ongoing work is to combine projection-based iterative phase retrieval methods with PRAMMol. This would allow low-resolution reconstructions to be enhanced at certain voxels with atomic-resolution reconstructions.

Additionally, PRAMMol’s accuracy is highly reliant on the accuracy of the interatomic potential used. Often, embedded atom model (EAM) interatomic potentials have difficulty modeling material properties to high fidelity. Because of this, these potentials may be insufficient to accurately reconstruct samples with PRAMMol. In these cases, machine-learned interatomic potentials (MLIPs) are an effective replacement because they can model material properties to arbitrary precision, albeit with a penalty to computation time.

In summary, this work presents the first practical CDI algorithm that can achieve atomic resolution from a diffraction pattern alone. The previous work by Dietze et al. started with the known phases from their simulated samples and then added noise to determine the effect of limited photon flux¹⁰. However, as we show here, PRAMMol is capable of retrieving the correct atomic positions, number of atoms, and scaling factor from the diffraction intensities alone. Additionally, PRAMMol can obtain the atom positions to extreme accuracy (limited by Poisson noise and numerical precision). Physics-based iterative phase retrieval opens the possibilities of imaging materials at the atomic scale at upgraded synchrotron facilities, substantially advancing our ability to take advantage of their unprecedented photon flux in more efficient and powerful ways.

Methods

PRAMMol algorithm

PRAMMol assumes that the measured intensity, I, of the diffraction pattern, at the position described by the Miller indices hkl in reciprocal space, is a random variable such that

$${I}_{{{{\rm{hkl}}}}} \mathop{\sim}\limits^{ind}Pois({\lambda }_{{{{\rm{hkl}}}}}(C,\overrightarrow{{{{\bf{x}}}}},\overrightarrow{{{{\bf{y}}}}},\overrightarrow{{{{\bf{z}}}}})),$$

(1)

where

$${\lambda }_{{{{\rm{hkl}}}}}(C,\overrightarrow{{{{\bf{x}}}}},\overrightarrow{{{{\bf{y}}}}},\overrightarrow{{{{\bf{z}}}}})=| {F}_{{{{\rm{hkl}}}}}{| }^{2}={\left| \sum\limits_{j = 1}^{n}C{{{{\rm{e}}}}}^{-2\pi {{{\rm{i}}}}(h{x}_{{{{\rm{j}}}}}+k{y}_{{{{\rm{j}}}}}+l{z}_{{{{\rm{j}}}}})}\right| }^{2}$$

(2)

where n is the number of atoms in the sample, F_hkl is the structure factor and is proportional to the electric field at the detector, C is an unknown scaling constant that describes the TIPF and several other experimental factors, and (x_j, y_j, z_j) is the j^th of n atomic positions in the material being imaged. In this formulation, the atomic positions and scaling constant become parameters of a family of probability density functions (Table 1).

Table 1 Common notation used throughout the paper

The simulated diffraction patterns have ignored experimental sampling conditions involving Ewald sphere curvature and the errors that can be introduced in BCDI due to the Fourier slice theorem that have been discussed extensively in the literature^43,44,45. We also assume the kinematic diffraction limit, which we feel is a good approximation for such small crystals as simulated here. Because of their size, we do not expect multiple scattering, absorption, or dynamical diffraction effects to have much of an effect on the simulations. We have also neglected effects such as detector pixel size and orientation. While these experimental challenges will need to be addressed in further application of the work, these issues and methods to model and account for them have been discussed extensively in the literature^17,18,44.

In conjunction with a low-resolution reconstruction, PRAMMol uses maximum likelihood estimation (MLE), cross validation and MD to find estimated atom positions. PRAMMol first finds estimates for atom positions by maximizing the likelihood function (Parameter Estimation: MLE). After maximizing the likelihood, cross validation selects a valid model that results in the highest likelihood on a validation set (Model Selection: Cross Validation). Valid models include any subset of the current atom positions and new atom positions that are generated with MD (Atom Generation). In general, MLE and cross validation alone are not enough to constrain the search space and arrive at the true solution, so this process is performed in conjunction with MD energy minimization (MD Energy Minimization). This three-step process is shown in Fig. 1 and pseudocode is laid out in Algorithm 1. These methods are performed iteratively until no new models are selected during cross validation across 100 iterations.

Algorithm 1

PRAMMol Algorithm

 1: for i in 1:NumIter do

 2: tDiffract ← 0.5% of Diffraction Pattern

 3: vDiffract ← 99.5% of Diffraction Pattern

 4: x,y,z ← BFGS(x,y,z,tDiffract) ⊳ 10 iters. of BFGS maximizing Eqtn. (6)

 5: C ← getScaling ⊳ “Scaling Constant Optimization”

 6: newX,newY,newZ ← GenerateAtoms ⊳ “Atom Generation”

 7: while Object Changes do  ⊳ Loop until the object doesn’t change

 8: ⊳ “Model Selection: Cross Validation”

 9: Sort(x,y,z,vDiffract) ⊳ Sort object atoms by likelihood

10: Sort(newX,newY,newZ,vDiffract) ⊳ Sort new atoms by likelihood

11:  while Likelihood Increases do

12: AddRemoveAtoms(x,y,z,newX,newY,newZ,vDiffract) ⊳ Fig. 9

13: end while

14:  end while

15 RunLammps(x,y,z) ⊳ Fix 50% of object and relax the rest

16: end for

Parameter estimation: maximum likelihood estimation

In statistical theory, point estimation methods are a class of methods that estimate the values of parameters given data that has been sampled. MLE is one such point estimation method that has several useful properties. Maximum likelihood estimates (MLEs) are both consistent and efficient estimators, meaning that they are asymptotically unbiased with a variance that matches the Cramer-Rao lower bound⁴⁶. Because the distribution we wish to estimate is an exponential family and the parameter space contains an open set, the MLE is the minimum variance unbiased estimator. When applied to phase retrieval problems, MLE has been successfully implemented with ptychography^47,48, holography⁴⁹, and speckle imaging⁵⁰.

Assuming the data is sampled independently and the likelihood function is continuous, the definition of the MLE $\hat{\theta }$ is

$$\hat{\theta }=\max\limits_{\theta \in \Theta }L(\theta | \overrightarrow{\iota }),$$

(3)

where Θ is the full parameter space and $\ell (\theta | \overrightarrow{\iota })$ is the likelihood function defined by

$$L\left(\theta | {\vec{\iota}}\right) = \prod\limits_{j=1}^n f_{I_{{{\mathrm{j}}}}}\left(\iota_{{{\mathrm{j}}}} | \theta\right),$$

(4)

where $f_{I_{{{\mathrm{j}}}}}(\iota_{{{\mathrm{j}}}} \vert \theta)$ is the probability density function of the jth data point. Intuitively, the likelihood function quantifies how likely one parameter choice is to fully describe. The supremum of the likelihood is then the most likely value of the parameters given the measured data. In practice, it is often easier to solve the equivalent problem

$$\hat{\theta }=\max\limits_{\theta \in \Theta }\ell (\theta | \overrightarrow{\iota }),$$

(5)

where $\ell (\theta | \overrightarrow{\iota })$ is the log of the likelihood.

From equations (1) and (2), it follows that the log-likelihood for our problem is

$$\ell_{\vec{{\mathbf{I}}}}(C, {\vec{\mathbf{x}}}, {\vec{\mathbf{y}}}, {\vec{\mathbf{z}}} \vert {\vec{\iota}}) = \sum_{{\mathrm{hkl}}} \iota_{{\mathrm{hkl}}} \ln \lambda_{{\mathrm{hkl}}}(C, {\vec{\mathbf{x}}}, {\vec{\mathbf{y}}}, {\vec{\mathbf{z}}}) - \lambda_{{\mathrm{hkl}}}(C, {\vec{\mathbf{x}}}, {\vec{\mathbf{y}}}, {\vec{\mathbf{z}}}).$$

(6)

Note here that ι_hkl is an intensity measurement and λ_hkl is a theoretical value of the intensity at a position described by the Miller indices (h, k, l) in reciprocal space. While it is impossible to find the global maximum of ${\ell }_{\overrightarrow{{{{\bf{I}}}}}}(C,\overrightarrow{{{{\bf{x}}}}},\overrightarrow{{{{\bf{y}}}}},\overrightarrow{{{{\bf{z}}}}}| \overrightarrow{\iota })$ analytically, numerical methods can be used to find maximizers of the log-likelihood. PRAMMol uses BFGS⁵¹, a quasi-newton minimization technique to maximize the likelihood function.

However, the likelihood function has many local maxima and BFGS will in general fail to find the global maximum. Additionally, this MLE algorithm does not allow for a changing number of atoms. To overcome local maxima and find the correct number of atoms, PRAMMol uses cross validation as a model selection technique. Thus, BFGS will only be used on a portion of the diffraction pattern and the rest will be held out to be used as a validation set in model selection. Empirically, we have found that the ideal portion of the diffraction pattern to use in the MLE search is around 0.5%, with 99.5% of the diffraction pattern used in cross-validation.

Model selection: cross-validation

Model selection can be viewed as a complement problem to parameter estimation. In parameter estimation, parameters are selected that are most likely to have generated the data. Model selection attempts to select the most likely model to have generated the data, rather than which parameters. Some common model selection procedures are cross-validation, the likelihood-ratio test, and information criterion, such as the Akaike information criterion (AIC) or the Bayesian information criterion (BIC)⁵². For PRAMMol, cross-validation is the most applicable model selection method because of the complexity and self-regularization of the models under consideration²⁹.

The specific cross-validation fitness metric PRAMMol uses is again the likelihood. For every model under consideration, the one considered the most optimal has the largest likelihood over the validation set. To create multiple possible models for cross-validation to consider, new atom positions are generated with MD as described in the Atom Generation (Atom Generation). These new atom positions are combined with the atom positions found via MLE. The set of models that PRAMMol considers during cross-validation is then any subset of this combined set of atom positions. This is an extremely large set of models that grows exponentially as the number of atoms increases, so cross-validation is only able to compare a subset of all models. To effectively search the space of possible models, PRAMMol uses a grid search algorithm.

The grid search method used in PRAMMol considers two sets, the atom positions in the current object and the atom positions out of the object. The set of atom positions in the object are ordered according to their effect on the likelihood when removed from the object and the set of atom positions out of the object are ordered according to their effect on the likelihood when added to the object. This gives an estimate of which atoms lead to a larger likelihood.

The grid search then attempts to increase the likelihood by removing the worst performing atom in the object, adding the highest performing atom, doing both, or doing nothing. The option with the highest likelihood is chosen and the process is repeated, except that the grid search allows for previous options to be undone. As Fig. 9 shows, this can be thought of as moving across a two dimensional grid and the position is allowed to move one space in any direction. The position on one axis is where to split the set of atoms in the object, with one side of the split staying in the object and the other being removed. The position on the other axis is where to split the atoms out of the object, with one side of the split staying out of the object and the other being added.

Fig. 9: During cross validation, PRAMMol uses two lists of atomic positions, the atoms in the current object and new potential atom positions.

These two are ordered according to their effect on the likelihood, and the cross-validation grid search attempts to find the best point to split the two lists. a Shows one possible split of the two lists, with the indicator showing where to split in both lists. The grid search then has nine choices, three for each list. b Shows a possible move from (a).

When the grid search finds that the best option changes nothing during the search, the grid search terminates and atoms are removed from and added to the object. Then, the new sets of atoms in the object are again ordered by their individual effect on the likelihood, and the grid search is performed again. This process is repeated until the grid search cannot find any improvements after restarting.

MD energy minimization

While the combination of cross-validation and MLE is an effective way to search the parameter space for the correct atomic positions, additional knowledge about the physical system can be added to make the search more effective. To do this, energy minimization using LAMMPS^40,41,42 is applied to a portion of the object keeping all other atoms constant. To choose atoms for energy minimization, PRAMMol first selects all atoms above a user defined energy cutoff. Specific values of an energy cutoff are relative and subject to experimental conditions. For our examples, the cutoff energy used is half the energy per atom of a tetrahedron in free space.

After the high energy atoms are selected, PRAMMol then randomly selects a large fraction of the total number of atoms. A good fraction to choose is around 50% of the object. These atoms then have their energies minimized with MD. This energy minimization process ensures that the current set of atoms remains physically reasonable, but does not enforce this condition too strongly by operating only on a portion of the object. This allows the algorithm to escape local maxima in the likelihood landscape.

Atom generation

The purpose of atom generation in PRAMMol is to generate all possible atom positions that may increase the likelihood of the current atom positions and be physically realistic. To generate atoms for use in cross-validation, PRAMMol uses three methods: Delaunay triangulation, force calculations, and random placement. After these three methods of generation, the atoms are filtered, removing any that are too close to each other or that do not fall in the support of a low-resolution reconstruction. Then, one by one, these atoms are combined with the current object and the energy of these atoms is minimized keeping the current object fixed. After this, the atoms are again filtered for any that are too close to each other or that do not fall in the support of a low-resolution reconstruction.

PRAMMol uses Delaunay triangulation to test atoms positions to potentially add. In a triangulation of three dimensions, each tetrahedron will connect four points. A Delaunay triangulation ensures that a sphere passing through these four points will not encompass any other point of the triangulation. For a given tetrahedron in a Delaunay triangulation, we are then guaranteed that the distance between the sphere’s center and any point in the triangulation is no smaller than the radius of the sphere.

For a metal with a predictably repeating pattern of atom locations, the radius of the circumcenter can be used to locate missing atoms or deviations from the usual arrangement of atoms. For an infinitely large crystal with no defects, the radius of every sphere in a Delaunay triangulation will be exactly the same. If, however, there is a missing atom, the radius will be significantly larger. This gives an excellent way to locate empty spaces in the current object and fill them with prospective new atom positions. Additionally, it will effectively cover most of the object in atom positions without having to search over a large space.

Another way to find missing atoms inside or outside of the object is to look at the atomic forces on each atom. After finding a local maximizer of the likelihood, the atoms will generally not be in a physically realistic location but rather will move to locations where the atoms will be most likely to create the diffraction pattern. This means that there will be atomic forces on the atoms, and in many cases, these forces will be in the direction of the location of the missing atoms. For example, if one atom is missing from the true object in the middle of a perfect crystal, then it is likely that the local maximum of the likelihood will place atoms around the missing one on the lattice positions themselves. This will then create a force on the atoms pointing at the location of the missing one.

For both Delaunay triangulation and force calculation atom generation, the placement schemes are effective but not perfect. To supplement this, atoms are placed randomly close to the object. This allows, over time, for any locations that Delaunay triangulation and force calculations do not reach to be sampled. The combination of these three methods efficiently searches the full range of physically realistic locations for cross-validation.

Scaling constant optimization

Throughout this paper, the positions of the atoms have been thought of as parameters in a probability density function from which the intensity values are sampled. While this is true of the scaling constant as well, an analytical form of the optimal scaling constant can be found when attempting to minimize the squared error of the simulated diffraction pattern rather than maximize the likelihood. Because of this, we do not treat the scaling constant as a parameter to be solved for in MLE and instead use an analytical form that minimizes the mean squared error, E, which is defined as

$$E=\sum\limits_{{{{\rm{hkl}}}}}{\left(\sqrt{{\iota }_{{{{\rm{hkl}}}}}}-| {F}_{{{{\rm{hkl}}}}}| \right)}^{2}$$

(7)

where ι_hkl, F_hkl, and h, k, l are all defined as in equation (2). E can then be minimized with respect to the scaling constant by taking a derivative and setting this equal to zero, as follows.

$$\frac{dE}{dC}=\sum\limits_{{{{\rm{hkl}}}}}-2| {\Omega }_{{{{\rm{hkl}}}}}| (\sqrt{{\iota }_{{{{\rm{hkl}}}}}}-C| {\Omega }_{hkl}| )=0$$

(8)

$$C=\frac{{\sum }_{{{{\rm{hkl}}}}}\sqrt{{\iota }_{{{{\rm{hkl}}}}}}| {\Omega }_{{{{\rm{hkl}}}}}| }{{\sum }_{{{{\rm{hkl}}}}}| {\Omega }_{{{{\rm{hkl}}}}}{| }^{2}}$$

(9)

where

$${\Omega }_{{{{\rm{hkl}}}}}=\sum\limits_{j=1}^{n}{{{{\rm{e}}}}}^{-2\pi {{{\rm{i}}}}(h{x}_{{{{\rm{j}}}}}+k{y}_{{{{\rm{j}}}}}+l{z}_{{{{\rm{j}}}}})}$$

(10)

The scaling constant is calculated once per iteration after the MLE step.

Grain creation

In this work, three samples were generated, a sample with no defects, a sample with a single vacancy, and a sample with a screw dislocation. All of these samples were simulated with atomsk⁵³ using the polycrystal tutorial on the atomsk website and deleting all but one grain. The sample with a single vacancy was then modified by removing one atom at the center and the sample with the screw dislocation was created using the screw dislocation tutorial on the atomsk website. After these grains were created, energy minimization was run in LAMMPS^40,41,42 to bring the samples to an equilibrium state (see Fig. 2).

Diffraction pattern simulation

For all reconstructions performed, the diffraction pattern was simulated over three peaks: the $\left(\bar{1},1,1\right)$, the $\left(1,\bar{1},1\right)$ and the $\left(1,1,\bar{1}\right)$ peaks. For each of these peaks the absolute value of the miller indices sampled ranged from 0.4 to 1.6 in h, k and l with 50 measurement points in each direction. Each hkl value sampled was drawn from a Poisson distribution as defined in Equation (1) (see Fig. 3).

Resolution calculation

To assess the effectiveness of PRAMMol, the sample with one vacancy was reconstructed 10 times across 20 values of TIPF. For one value of TIPF, the reconstructed diffraction patterns of the 10 trials were averaged over constant shells of ∣q∣. Similar to ref. ¹⁰, the resolution is defined as π/q_c, where q_c is the smallest value of q where the mean absolute error of the phase prediction along the shell of constant ∣q∣ is greater than 0.5.

Temperature-dependent diffraction pattern

To simulate a diffraction pattern from a sample at a finite temperature, the sample first underwent a canonical (NVT) ensemble simulation at the temperature of interest for 110,000 iterations. After 10,000 iterations had run, snapshots were collected every 100 iterations for a total of 1000 snapshots. After this, the measured intensity, I, is assumed to be a random variable such that

$${I}_{{{{\rm{hkl}}}}} \mathop{\sim}\limits^{ind}Pois({\bar{\lambda }}_{{{{\rm{hkl}}}}}(C,\overrightarrow{{{{\bf{x}}}}},\overrightarrow{{{{\bf{y}}}}},\overrightarrow{{{{\bf{z}}}}})),$$

(11)

where

$${\bar{\lambda }}_{{{{\rm{hkl}}}}}(C,\overrightarrow{{{{\bf{x}}}}},\overrightarrow{{{{\bf{y}}}}},\overrightarrow{{{{\bf{z}}}}})={\left| \frac{1}{1000}\sum\limits_{s = 1}^{1000}\sum\limits_{j = 1}^{n}C{{{{\rm{e}}}}}^{-2\pi {{{\rm{i}}}}(h{x}_{{{{\rm{sj}}}}}+k{y}_{{{{\rm{sj}}}}}+l{z}_{{{{\rm{sj}}}}})}\right| }^{2},$$

(12)

where s corresponds to one of the 1000 snapshots.