# Score-based denoising for atomic structure identification

Tim Hsu<EMAIL_ADDRESS>Lawrence Livermore National Laboratory, Livermore, CA,
United States  Babak Sadigh  Lawrence Livermore National Laboratory,
Livermore, CA, United States  Nicolas Bertin  Lawrence Livermore National
Laboratory, Livermore, CA, United States  Cheol Woo Park  Lawrence Livermore
National Laboratory, Livermore, CA, United States  James Chapman  Boston
University, Boston, MA, United States  Vasily Bulatov  Lawrence Livermore
National Laboratory, Livermore, CA, United States  Fei Zhou<EMAIL_ADDRESS>Lawrence Livermore National Laboratory, Livermore, CA, United States

###### Abstract

We propose an accurate method for removing thermal vibrations that complicate
the task of analyzing complex dynamics in atomistic simulation of condensed
matter. Our method iteratively subtracts thermal noises or perturbations in
atomic positions using a denoising score function trained on synthetically
noised but otherwise perfect crystal lattices. The resulting denoised
structures clearly reveal underlying crystal order while retaining disorder
associated with crystal defects. Purely geometric, agnostic to interatomic
potentials, and trained without inputs from explicit simulations, our denoiser
can be applied to simulation data generated from vastly different interatomic
interactions. Followed by a simple phase classification tool such as the
Common Neighbor Analysis, the denoiser outperforms other existing methods and
reaches perfect classification accuracy on a recently proposed benchmark
dataset consisting of perturbed crystal structures (DC3). Demonstrated here in
a wide variety of atomistic simulation contexts, the denoiser is general,
robust, and readily extendable to delineate order from disorder in
structurally and chemically complex materials.

##  1 Introduction

In molecular dynamics (MD) of condensed matter, characterization methods for
the simulated atomic configurations aim to unravel meaningful structural
features such as crystalline phases and defects. As the simulations are
typically carried out at finite temperatures, accurate characterization of
structures and defects is complicated by perturbations in atomic positions
induced by thermal vibrations. To this end, increasingly sophisticated methods
have been proposed over the years for identification of local atomic motifs in
simulated configurations [1, 2, 3, 4].

Existing characterization methods usually focus on either ordered crystalline
phases or crystal defects. For example, the common neighbor analysis (CNA)
algorithm [2] identifies simple crystal structures such as the body-centered
cubic (BCC), face-centered cubic (FCC), and hexagonal closed-packed (HCP).
Other commonly used methods for structure identification include bond order
analysis [1, 5], centrosymmetry analysis [6], and adaptive template analysis
[7]. On the other hand, the dislocation extraction algorithm (DXA) [8, 9]
identifies dislocation defects within an a priori known ordered crystalline
environment. All of the mentioned methods rely heavily on domain knowledge,
physical intuition, and heuristics (for review, see for example [3]). As such,
they are often application- and/or structure-specific, and are not always easy
to generalize beyond their original scope of applicability. More recently,
data-driven machine-learning (ML) approaches are being developed for
performing ordered phase classification and sometimes defect detection [10,
11, 12, 13, 14, 15, 16, 17], often employing existing tools such as Steinhardt
order parameters [1] for featurization. While comparatively more
straightforward to develop with modern ML pipelines, these emerging methods
require considerable amounts of carefully curated training data and are often
informed by material-specific physics and domain knowledge which limit
transferability of the trained models.

In this work, we consider the tasks of classifying ordered crystalline phases
and identifying disordered defects to be fundamentally of the same nature,
akin to determining the order-disorder feature signal from thermal noise. We
address this unifying problem with an iterative optimization scheme that aims
to eliminate thermal noise. By minimizing thermal perturbation, our method
significantly simplifies subsequent characterization of both ordered
crystalline phases and crystal defects, thereby placing the two seemingly
distinct tasks on an equal footing.

Figure 1: Our denoising graph network model predicts and subtracts thermal
perturbations in atomic positions. Denoising proceeds iteratively until the
predicted noise becomes approximately zero or after a pre-defined number of
iterations. In this schematic the atoms are shown in 2D Cartesian space for
ease of visualization.

To implement such a denoiser, we trained a graph network model, based on the
equivariant NequIP architecture [18, 19], to denoise heavily perturbed
structures and reveal the underlying order-disorder (Fig. 1). Given a
distribution of pairs of noiseless ideal reference topologies (e.g., BCC, FCC,
and HCP) and their noised (perturbed) counterparts, training our denoiser is
equivalent to learning a score function [20], which in this work is a gradient
field in the atomic coordinate space converging towards points of maximum
likelihood that correspond to the ideal reference topologies (more details in
Section 4.1 and Appendix A). Equipped with this theoretic knowledge, our
denoiser can be considered as an iterative scheme that optimizes perturbed
structures towards ideal topologies. The score function plays a central role
in modern generative models such as the denoising diffusion probabilistic
model (DDPM) [21, 22, 23] for sampling realistic data from a high dimensional
data distribution [24]. Here it is applied for denoising rather than
generative applications. Compared to other existing denoising methods, e.g.
energy minimization (or steepest descent mapping [25]) and vibration-
averaging, our method is unbiased in the sense that the ideal reference
structures have the same probability, and only requires the instantaneous
snapshot as the input. In contrast, energy minimization is intrinsically
biased due to the use of interatomic potential that may favor certain phases
over others, and vibration-averaging requires a tuned averaging window over
multiple snapshots where fast processes may be overly smeared. Further,
utilizing synthetically noised structures as training data, our approach does
not rely on in-depth physical knowledge other than the ideal reference
structures, and is a purely geometric algorithm complementary to existing
physics-based techniques. As such, in contrast to data-hungry approaches, our
denoiser does not need physics simulation data for training.

Prioritizing single-element systems in this work, the denoising capabilities
of our model is demonstrated on several challenging applications, including
identification of transient crystal phases during Cu solidification from melt,
and characterization of dislocations and point defect debris in BCC Ta
undergoing plastic deformation. Importantly, our denoiser does not
overzealously denoise the disordered melt into ordered phases. Further, it is
shown to help reveal and locate point defects, dislocations, and grain
boundaries at high temperatures (approaching melting point) where again the
model is observed not to denoise or rearrange crystal defects into ideal
lattice motifs. Additional demonstrations on two-element SiO2 polymorphs are
also provided. Besides denoising, the underyling neural network architecture
of the denoiser can be extended to classify the denoised atomic environments.
At this stage of development we rely instead on existing methods, mainly CNA
and DXA, to perform the final classification. With appropriate optimization,
we envision that our denoising algorithm would be a robust and highly
efficient filter integrated in the workflows of massive MD simulations for the
purpose of on-the-fly data compression and post-processing analyses.

Iterative applications of score functions allow us to approach perfect
identification of ordered crystalline phases in several case studies with
significant improvement over exiting classifiers with a single application. It
is interesting that similar advantages of iterative use of score functions
were observed in DDPM in an annealed Langevin dynamics setting for generative
purpose [21, 22, 23]. Our findings support the view that iterative models
break down challenging problems into smaller, manageable steps.

##  2 Results

In our approach, denoising a thermally perturbed configuration entails
iteratively subtracting the noises predicted by a machine learned graph
network model (Fig. 1):

${\mathbf{r}}\leftarrow\mathbf{r}-\epsilon_{\theta}(\mathbf{r},\mathbf{z}),$
(1)

where $\mathbf{r}$ is the atomic coordinates, $\epsilon_{\theta}$ is the graph
network model with parameters $\theta$, and $\mathbf{z}$ holds auxiliary
information such as the atom types and unit cell dimensions. Repeated
application of Eq. 1 amounts to finding a converged fixed point corresponding
to a fully denoised structure. Our denoiser is an optimization algorithm that
topologically modifies input (noised) structures towards maximal data
likelihood (further explained in Section 4.1 and Appendix A). By including the
ideal FCC, HCP, and BCC lattices in the training data, our model attempts to
evolve an input noised structure towards one of the three ideal lattices
depending on which lattice type it most resembles topologically. Importantly,
as demonstrated in the results that follow, our denoiser does not excessively
alter the topology of disordered structures, including liquid/melt phase,
point defects, dislocations, and grain boundaries far removed from the ideal
lattice topologies. This property renders our method safe against overzealous
denoising, thus retaining meaningful disordered features in input structures.

Trained with purely synthetic data, our denoiser is applied to a wide variety
of MD-perturbed systems: (1) BCC, FCC, and HCP Cu simulated above the melting
point, as well as liquid/melt Cu perturbed around the melting point; (2) hard-
to-detect transient crystal phases momentarily forming during solidification
of Cu from melt; (3) FCC, HCP, and BCC Cu containing point defects; (4) BCC Ta
containing complex dislocation networks and point defect clusters; (5) BCC Ta
containing grain boundaries; (6) the recently published “DC3” benchmark
dataset for crystal structure identification [15] containing Al (FCC), Ar
(FCC), Fe (BCC), Li (BCC), Mg (HCP), and Ti (HCP) perturbed at the melting
point; and finally (7) SiO2 polymorphs $\beta$-quartz, $\alpha$-cristobalite,
and $\beta$-cristobalite. In (1)–(5), the denoiser is shown to reduce or
eliminate thermal noise, making it trivial to identify the underlying crystal
structures while not destroying meaningful disordered features such as point
defects, dislocations, and grain boundaries. In (6), followed by CNA
classification, the denoiser achieves perfect classification accuracies in all
systems. In (7), the generalizalibity of the denoiser to multi-element complex
materials is validated. The results for each case study are detailed below.

###  2.1 Denoising FCC, HCP, BCC, and liquid/melt Cu

The first demonstration focuses on denoising solid BCC, FCC, and HCP crystals
perturbed by thermal vibrations at above the melting point (3400 K), and
liquid/melt Cu at around the melting point (3000 K), as shown in Fig. 2a.
Before denoising, the popular adaptive CNA (a-CNA) algorithm [3] classifies
most of the solid atoms (82, 77, and 77% in BCC, FCC, and HCP, respectively)
as disordered, i.e. not belonging to any of the three crystal lattice types.
After just one iteration of denoising, the number of misclassified atoms is
significantly reduced (14, 5, and 3% for BCC, FCC, and HCP, respectively). The
following iterations, typically within 5–8 steps, remove the remaining minor
perturbations. The denoised solids resemble perfect FCC, HCP, and BCC
lattices, thus trivializing subsequent phase classification. The Steinhardt
order parameters $\bar{q}_{4}$ and $\bar{q}_{6}$ [5] computed before and after
denoising confirm that virtually all thermal perturbations imparted on the
solids are removed (Fig. 2b). Note that after denoising, the Steinhardt
quantities for BCC/FCC/HCP appear to be a single point, but are in fact about
1,000 points (each point corresponds to an atom) overlapped together.

Figure 2: Iterative application of the denoiser to small cells of BCC, FCC,
HCP, and liquid/melt Cu. (a) Visualization of the structures along the
denoising iterations. (b) Steinhardt features $\bar{q}_{4}$ and $\bar{q}_{6}$
[5] before and after the denoising. (c) Radial distribution function of the
melt phase before and after the denoising. In (a), the structures are shown in
ortho (orthogonal) views along densely packed crystallographic directions,
with additional perspective 3D views for steps 0 and 8. The atoms are colored
according to a-CNA prediction implemented in OVITO [26]. The solid phases and
the melt have been annealed at 3400 K and 3000 K, respectively. ICO stands for
icosahedral coordination.

Interestingly, denoising the Cu melt phase leaves nearly all atoms to remain
disordered, as indicated by a-CNA labeling all such atoms as other or unknown
(Fig. 2a), even though the atomic displacements over the denoising iterations
are roughly the same as that for the solid phases (Supplemental Fig. 2).
Additionally, the very first peak of the radial pair distribution function
(RDF) becomes sharper and splits after denoising, and the peaks at the medium-
range distances also become slightly sharper (Fig. 2c). Further analyses in
Appendix B leave precise effect of the denoiser on melt structure uncertain
and perhaps deserving further scrutiny in future work. The fact that the Cu
melt phase remains disordered even after denoising can be a useful property of
our model not explicitly learned from its training data, which consists of
only perfect and randomly distorted but otherwise ordered crystal lattices. In
simulations involving solid-liquid coexistence, an example of which will be
shown next, we certainly wish our model to denoise only the thermally
distorted crystal lattices while leaving truly disordered phases disordered.

###  2.2 Denoising Cu solidification trajectory

Our denoiser is further tested here on an MD trajectory of Cu solidification
from melt. Previously studied by Sadigh, et al. [27], solid nuclei appearing
in the initial transient stages of Cu solidification are polymorphic,
containing BCC, FCC, HCP, and disordered melt phase simultaneously.
Characterization of such a complex transient behavior is challenging and
presents a useful test case for our method. As shown in Fig. 3a, denoising
four transient configurations of the trajectory results in drastic improvement
in subsequent phase classification by the a-CNA algorithm. This improvement is
manifested in much denser, correctly classified labels (FCC, HCP, and BCC) on
the atoms within the ordered solid nucleus, as well as in considerably sharper
boundaries between the phases. Here again the atoms in the disordered melt
remain disordered as their a-CNA labels remain largely unchanged. Notably, a
few atoms labeled as other (i.e., unknown) are observed within the sharply
defined crystal phases even after denoising. These “unknown” atomic motifs are
likely point defects, e.g. vacancies and interstitials, as will be discussed
in greater detail in the next section. Also, they appear mostly in the BCC
phase likely due to the metastability of the BCC phase in Cu entropically
stabilized under the high pressure (70 GPa) of the simulation. Otherwise the
same BCC phase of Cu is only marginally metastable as manifested in the
appearance of soft modes in its phonon spectrum. At high temperatures close to
melting, such soft modes may well result in some of the atoms within the BCC
phase to significantly deviate from their ideal lattice positions, resulting
in formation of point defects.

Figure 3: Denoising a dynamic trajectory of Cu solidification (314,926 atoms)
significantly improves (a-CNA) phase classification without introducing
unreasonably short interatomic distances. (a) Four consecutive snapshots along
the trajectory, in original, denoised, and energy minmized states, are shown
in ortho views with an additional perspective view for the fourth snapshot.
(b) RDF of the fourth snapshot before and after denoising. In the perspective
view, the atoms classified as other by a-CNA are rendered transparent.

It can be informative to compare denoising to steepest descent energy
minimization (EM), another common method for filtering out thermal vibrations.
In stark contrast to denoising, EM not only reduces the thermal perturbation
but also greatly changes the nucleus structure beyond recognition (third row
of Fig. 3a). Namely, under EM the solid nucleus grows considerably larger, and
the transient BCC phase disappears in favor of the more stable FCC and HCP
phases. The striking difference between EM and denoising can be attributed to
their contrasting assumptions. The EM mapping clearly favors solid phases of
lowest ground state energy and may overzealously nudge atoms towards such
phases, exactly as it happens in the considered example. The denoising
process, from a complementary and purely geometric perspective, relies on an
unbiased, equal prior probability of the reference phases included in
training.

To further investigate whether denoising introduces unwanted or unphysical
artifacts, in Fig. 3b we plot the RDF of the last snapshot from Fig. 3a. The
RDF of the denoised structure generally matches that of the original. Having
been trained to reduce thermal perturbations, the denoiser likely has also
learned not to bring atoms to excessively short distances of each other (that
would be unphysical). The sharp peaks of the denoised RDF are attributed to
the three ordered crystal phases, each of which contributes its own discrete
set of sharply defined interatomic distances.

Note that the denoiser does in fact “denoise” some atoms within the melt phase
into local environments regarded as crystalline (BCC/FCC/HCP) by the a-CNA
classification, as evidenced by the increased number of solid labels in the
melt region after denoising (Fig. 3a). The appearance of such “crystalline”
atoms reflects that even in fully disordered liquid, statistically a small
fraction of the atomic motifs may momentarily resemble a crystal phase. Our
denoiser then acts locally and further enhances such resemblance, thus making
a-CNA (a strictly local classifier itself) to recognize such atoms as
crystalline. Observing that these misclassified atoms are few and far
isolated, we accept this minor artifact as a worthwhile trade-off. A more
elaborate classifier, perhaps accounting for local environments beyond just
the nearest neighbors, should be able to cleanly separate the ordered and the
disordered melt phases in such instances.

###  2.3 Denoising FCC, HCP, and BCC Cu containing point defects

As an example of point defect characterization, the denoiser is applied to
FCC, HCP, and BCC Cu each containing an intentionally inserted extra atom
followed by annealing in MD at 3400 K. These structures were denoised into
ideal lattices with local regions of disorder unknown to a-CNA (Fig. 4). The
Wigner-Seitz defect analysis (OVITO) confirms that these regions of disorder
indeed correspond to point defects. Notably, for the FCC crystal, the thermal
vibrations are sufficiently intense to spontaneously generate two more point
defects, a Frenkel pair of one vacancy and one interstitial.

Figure 4: Denoising facilitates point defect identification in thermally
perturbed FCC, HCP, and BCC Cu crystals. In the original configurations, the
point defects are completely obscured by thermal vibrations but are cleanly
revealed after denoising (circled in black). When applied to the denoised
structures, the Wigner-Seitz defect analysis implemented in OVITO (rightmost
column) correctly assigns mass content to each crystal defect: site
occupancies 0, 1, and 2 correspond to a vacancy, a regular atom, and an
interstitial, respectively. After denoising, the regular atoms in the FCC
crystal (top row) are rendered semi-transparent to more clearly reveal the
point defects.

The example in Fig. 4 demonstrates desirable outcomes of denoising crystal
structures containing point defects. Although the denoiser aims to modify
local atomic motifs towards ideal topology, it cannot do so on regions of
point defects simply due to extra or missing atoms. In such a case, the model
appears to not significantly alter the local topologies around the defects
while denoising the rest of the bulk into an ideal lattice.

###  2.4 Denoising BCC Ta containing dislocations

Our model is similarly effective for denoising lattice dislocations. For a toy
example, a hexagon-shaped dislocation loop inserted into BCC Ta was annealed
at 2500 K ($0.8T_{m}$, where $T_{m}$ is the melting point) and subsequently
denoised (Fig. 5a). Again, satisfyingly, the denoiser does not significantly
alter local atomic configurations near the dislocation loop while cleanly
denoising the surrounding crystal bulk.

Figure 5: Denoising BCC single crystal Ta containing dislocations (and point
defects). (a) A single dislocation loop annealed at 2500 K. (b) A snapshot of
a relatively large MD simulation (1.97 million atoms) of single crystal Ta
undergoing deformation at 2000 K, which results in a network of entangled
dislocation lines and a large number of point defect clusters. The DXA method
for dislocation analysis (OVITO) confirms the presence of dislocations and
delineates the dislocation network topology. The cluster analysis (OVITO) aims
to separate the point defect clusters from the dislocations. The atoms are
normally colored by a-CNA prediction except for the cluster/dislocation
analysis, in which the BCC atoms are rendered transparent.

Note that although the denoiser was trained on ideal and noised Cu lattices,
it is applicable to Ta or any other elemental crystal of FCC, BCC or HCP
lattice structure. This transferability is achieved by simply scaling the
input (noised) structure to match the interatomic distance of the
corresponding Cu phase. The output (denoised) structure would then be re-
scaled back to its original dimensions.

As a more realistic and difficult test, the denoiser was applied to help
reveal a complex dislocation network in a BCC Ta crystal subjected to plastic
deformation at 2000 K (Fig. 5b). Similar to the case of the single dislocation
loop, the dislocation network (as colored by a-CNA) is more sharply defined
after denoising. Subsequent application of the DXA algorithm [9] to the
original and the denoised configurations results in nearly identical
dislocation networks, which testifies to the exceptional robustness of DXA
against thermal perturbation, and confirms that the dislocations are better
captured by a-CNA after denoising.

Despite DXA’s already high performance in dislocation characterization, the
denoiser still benefits or complements DXA by facilitating the
characterization of the point defect clusters that were either left as debris
in the wake of dislocation motion or produced by dragging jogs formed at
dislocation intersections [28]. Focusing on non-BCC atoms clustered within a
cutoff distance of 3.2 Å, a large cluster (corresponding to the dislocation
network) and a high concentration of small clusters are observed. Without
denoising, the small clusters may simply be manifestation of noise based on
visual interpretation. However, after denoising, the small clusters resemble
and likely capture the point defect debris. This is unlikely a case of the
denoiser failing to denoise non-dislocation regions into perfect lattice for
two reasons: (1) the presence of the point defect debris is known a priori,
and (2) the denoiser clearly denoises the non-dislocation region in the toy
case of the single dislocation loop, with virtually no point defects left
(Fig. 5a).

###  2.5 Denoising BCC Ta containing grain boundaries

To test how our method performs on crystals containing grain boundaries, the
denoiser was applied to a Ta bi-crystal containing two tilt boundaries. Prior
to denoising, the bi-crystal was annealed at a high temperature of 2500 K. As
shown in Fig. 6a, denoising does not alter the topology of the defects, and
results in two near-perfect BCC crystals separated by two perfectly planar
grain boundaries, with trace amount of point defects likely emitted from the
boundaries into the bi-crystal interior. A more complex test case is shown in
Fig. 6b, where a polycrystal consisted of 12 grains of BCC Ta had been
similarly annealed at 2500 K. Here again, denoising removes thermal noise
while still revealing a few point defects in the grain interiors.

Figure 6: The denoising process helps reveal grain boundaries obscured by
thermal noise. (a) A bi-crystal of BCC Ta (64,000 atoms) containing two planar
grain boundaries. (b) A polycrystal of BCC Ta (187,921 atoms) containing a
network of grain boundaries. Both crystals were annealed in MD simulations at
2500 K which caused minor coarsening in the polycrystal. After denoising,
trace amount of of point defects becomes visible in both examples. The BCC
atoms are rendered slightly semi-transparent.

###  2.6 Classification accuracy on the DC3 benchmark dataset

Recently Chung et al. [15] proposed an efficient ML approach, named data-
centric crystal classifier (DC3), for crystal structure classification. The
accuracy of this new method was tested on a dataset of configurations spanning
a wide variety of atomistic models, crystal structures, and temperatures
ranging from cryogenic to melting. Compared to other existing methods, DC3 was
shown to achieve the best accuracy on most systems in the dataset. The authors
offered the same dataset as a standard benchmark for crystal structure
classifiers (https://github.com/freitas-rodrigo/DC3). Here we test the
accuracy of our approach, namely denoising followed by a-CNA, on a subset of
the benchmark dataset: Al (FCC), Ar (FCC), Fe (BCC), Li (BCC), Mg (HCP), and
Ti (HCP) at their corresponding melting points (10 snapshots, or about 20,000
atoms for each category). Table 1 shows the classification accuracy comparison
between our approach and DC3 [15]. Within eight denoising iterations, the
performance of denoiser + a-CNA converges to the perfect classification score
of 100% on all but two Fe snapshots (Supplemental Fig. 6a). Confirmed by the
original authors, the two Fe snapshots each contains a Frenkel pair (a vacancy
and an interstitial) that must have formed spontaneously from intense thermal
fluctuations at the melting temperature. Thus, in keeping with its unique
ability to not only denoise thermally perturbed crystal structures but to also
clearly reveal crystal defects, our accuracy is deemed perfect across the
board.

Table 1: Accuracy comparison between our classification approach (denoiser +
CNA) and the DC3 framework [15] on the benchmark dataset therein (only the
BCC, FCC, and HCP phases were considered). 8 iterations were used for
denoising. The accuracy value, computed as the fraction of correctly labeled
atoms, is shown in percentage.

| Al (FCC) | Ar (FCC) | Fe (BCC) | Li (BCC) | Mg (HCP) | Ti (HCP)
---|---|---|---|---|---|---
Denoiser (1 step) + a-CNA | 98.7 | 99.8 | 97.8 | 98.3 | 99.4 | 91.4
Denoiser (8 steps) + a-CNA | 100 | 100 | 99.7* | 100 | 100 | 100
DC3 | 96.9 | 97.5 | 86.8 | 85.8 | 97.4 | 89.4

*Without the point defects in BCC Fe (Supplemental Fig. 6b), the accuracy would be 100%.

###  2.7 Denoising SiO2 polymorphs

Finally, to validate the generalizability of our approach, we trained a
separate denoiser for minimizing thermal perturbation in SiO2 systems. Similar
to the demonstration shown in Fig. 2, the denoiser was applied to MD-perturbed
high-temperature silica polymorphs $\beta$-quartz, $\alpha$-cristobalite, and
$\beta$-cristobalite (Fig. 7). Although we also trained the denoiser with
reference $\alpha$-quartz topology, the inherent similarity between $\alpha$\-
and $\beta$-quartz complicates the task of generating stable and distinctive
MD-perturbed snapshots of $\alpha$-quartz compared to $\beta$-quartz.
Therefore $\alpha$-quartz is not included in this preliminary result.
Regardless, the denoiser removes virtually all thermal perturbation in the
silica polymorphs. Future work may include extending to larger multi-element
systems in the presence of disordered features.

Figure 7: A separately trained denoiser effectively removes thermal
perturbations (simulated at 1000 K) in SiO2 $\beta$-quartz,
$\alpha$-cristobalite, and $\beta$-cristobalite polymorphs. The denoised
structures match the ideal references used during training. 16 iterations were
used for denoising. The silicon atoms are colored in beige, and the oxygens in
red.

##  3 Discussion

In developing our new method we have taken to an unbiased geometric and
statistical perspective on delineating order-disorder features in MD
simulations of solids, a problem particularly difficult at elevated
temperatures approaching the melting point. We regard two seemingly distinct
tasks of classifying ordered phases and locating disordered defects as
fundamentally the same problem, which is addressed by denoising. Based on a
statistical score function, the denoiser presented and tested in this work
effectively reduces and minimizes thermal noise in ordered solids without
impacting isolated disordered defects and the liquid/melt phase. To support
our conclusions, in the preceding sections we applied denoising to reveal the
underlying “anomaly” structures across the entire spectrum of crystal
disorder, namely (0D) point defects, (1D) dislocation lines, (2D) grain
boundary, and (3D) liquid phases, all distinct from the ordered phases that
the denoiser was trained with. For denoising BCC/FCC/HCP, our model was
trained only once and only on the three ordered structures but then shown to
successfully but not excessively denoise the much more complex structures used
for testing. We relate this useful ability of not overly denoising defective
atoms to the inductive bias of our graph network model with its limited number
of message-passing steps. Trained entirely on synthetic data, our model is not
derived from any deep physical insights about topology and geometry of the
reference structures and, as such, does not require careful data curation.
Thus, its extension to other ordered structures should be straightforward.

Our method holds unique advantages (and disadvantages) over two other methods
widely used for reducing thermal perturbation: energy minimization (EM) and
vibration-averaging (VA). With the right parameters and done over a small
number of iterations, EM can lead to similar results as that of the denoiser.
However, over many iterations, EM may grow or shrink certain phases due to the
intrinsic bias associated with the use of interatomic potential. An example of
EM being overzealous and distorting configurations beyond recognition with
respect to the original state is in Section 2.2. The purely geometric
denoiser, on the other hand, does not require a known or developed interatomic
potential. Further, the denoiser is unbiased in the sense that the reference
lattices used in training have equal prior probability. Although the denoising
graph model is more expensive than VA (which entails simple averaging
operations), VA can potentially smear out atomic motion by averaging over a
time interval. Such a distortion can be fairly significant since, to average
out thermal vibrations, VA requires time averaging intervals of hundreds or
even thousands of time steps. Our denoising method, on the other hand, treats
every time snapshot separately and does not coarse-grain over time. As is
often the case, no one method used for denoising is singularly superior to all
other existing methods. We hope that our approach finds its own place among
existing and emerging methods for structure and defect classification, and can
serve as an accurate and efficient pre-processing filter to facilitate
application of more computationally demanding methods of structural analysis
such as DXA.

Based on an equivariant graph network model architecture, our denosing model
is readily extendable to more complex reference structures and materials by
incorporating additional information such as atom types into the graph
embedding. In addition to extending the method to chemically complex systems,
our ongoing and future efforts may also focus on its computational efficiency
and scalability. Finally, beyond atomic structures, the ability of our model
to achieve state-of-the-art classification accuracy through iterative
denoising score functions suggest the idea may be useful in other disciplines
for enhanced accuracy and/or robustness against fluctuations.

##  4 Methods

###  4.1 Theoretical justification of the denoising model

Our approach to denoising builds on the theory of statistical learning of
score functions [20] that establishes equivalence between denoising and score
matching. Consider a probability distribution function $q(\mathbf{x})$ that
exists in principle but is analytically intractable due to the high
dimensionality of the data space $\mathbf{x}\in\mathbb{R}^{d}$ ($d\gg 1$).
Focusing on approximating the gradient of the log-probability density
$\nabla_{\mathbf{x}}\log q(\mathbf{x})$—also known as the score function
[29]—rather than $q(\mathbf{x})$ itself circumvents the often intractable
problem of finding the normalization constant for $q(\mathbf{x})$. Score
matching then amounts to finding an approximating model
$s_{\theta}(\mathbf{x})$ with parameters $\theta$ to match the score function,
with the score matching  loss [29]

$L_{\mathrm{SM}}=\frac{1}{2}\mathbb{E}_{q(\mathbf{x})}\left[\left\|s_{\theta}(\mathbf{x})-\nabla_{\mathbf{x}}\log
q(\mathbf{x})\right\|^{2}\right].$ (2)

Nevertheless, the term $\nabla_{\mathbf{x}}\log q(\mathbf{x})$ is still
unknown. To address this, consider approximating $q(\mathbf{x})$ by adding
isotropic Gaussian noises of variance $\sigma^{2}$ to the (clean) data samples
$\mathbf{x}$, resulting in noised samples
$\mathbf{x}^{\prime}=\mathbf{x}+\sigma\epsilon$, where
$\epsilon\sim\mathcal{N}(0,\mathbf{I})$, and the approximating distribution

$q_{\sigma}(\mathbf{x}^{\prime}|\mathbf{x})=Ze^{-\|\mathbf{x}^{\prime}-\mathbf{x}_{0}\|^{2}/2\sigma^{2}},$
(3)

where $Z$ is a normalization constant. This way, instead of the original loss,
we minimize the denoising score matching loss based on the key insight from
Ref. [20] to train with pairs of clean and corrupted data points:

$L_{\mathrm{DSM}}=\frac{1}{2}\mathbb{E}_{q_{\sigma}(\mathbf{x}^{\prime},\mathbf{x})}\left[\left\|s_{\theta}(\mathbf{x}^{\prime})-\nabla_{\mathbf{x}^{\prime}}\log
q_{\sigma}(\mathbf{x}^{\prime}|\mathbf{x})\right\|^{2}\right],$ (4)

where the new score function $\nabla_{\mathbf{x}^{\prime}}\log
q_{\sigma}(\mathbf{x}^{\prime}|\mathbf{x})$ can be computed via

$\displaystyle\log q_{\sigma}(\mathbf{x}^{\prime}|\mathbf{x})$
$\displaystyle=\log(Z)-\frac{1}{2\sigma^{2}}\|\mathbf{x}^{\prime}-\mathbf{x}\|^{2},$
$\displaystyle\nabla_{\mathbf{x}^{\prime}}\log
q_{\sigma}(\mathbf{x}^{\prime}|\mathbf{x})$
$\displaystyle=-\frac{1}{\sigma^{2}}(\mathbf{x}^{\prime}-\mathbf{x})=\frac{1}{\sigma^{2}}(\mathbf{x}-\mathbf{x}^{\prime})=-\frac{\epsilon}{\sigma},$
(5)

revealing that the score function points from noisy samples
$\mathbf{x}^{\prime}$ to clean ones $\mathbf{x}$. This observation also
implies that learning the score function is equivalent to training a denoising
model. To better see this, note that the denoising score matching loss can now
be simplified into

$L_{\mathrm{DSM}}=\frac{1}{2}\mathbb{E}_{q_{\sigma}(\mathbf{x}^{\prime},\mathbf{x})}\left[\left\|s_{\theta}(\mathbf{x}^{\prime})+\frac{\epsilon}{\sigma}\right\|^{2}\right].$
(6)

After scaling Eq. 6 by a factor of $\sigma$ and incorporating a noise
prediction model $\epsilon_{\theta}(\textbf{x}^{\prime})$ that aims to predict
the applied noise, then the loss function can be written as

$L_{\mathrm{DSM}}=\frac{1}{2}\mathbb{E}_{q_{\sigma}(\mathbf{x}^{\prime},\mathbf{x})}\left[\left\|\epsilon-\epsilon_{\theta}(\mathbf{x}^{\prime})\right\|^{2}\right],$
(7)

finally establishing the connection between the score function model
$s_{\theta}(\textbf{x}^{\prime})$ and the denoising model
$\epsilon_{\theta}(\textbf{x}^{\prime})$ by
$\epsilon_{\theta}(\textbf{x}^{\prime})=-\sigma
s_{\theta}(\textbf{x}^{\prime})$. This clarifies the meaning of
$\epsilon_{\theta}$ in Eq. 1: it is the scaled score function defining the
noise added to clean data. While the noise amplitude $\sigma$ is a hyper-
parameter, it can be estimated by a fitted/trained denoising model itself from
a noisy input. Hereafter, we refer to $\epsilon_{\theta}$ as a score function
for brevity. Along the above steps for connecting between score matching and
denoising, we have omitted certain details for brevity. For a rigorous
formulation, see [20].

The ideal score function is a gradient field in the data space that converges
to clean data points used to train the denoiser. For better intuition, a toy
score function is visualized in Supplemental Fig. 1, which illustrates that
following the score function is the same as denoising. In our context where
ideal FCC, HCP and BCC lattices were used for training, a perturbed (noised)
input configuration may be denoised into one of the three ideal (clean)
structures that it resembles the most. At the same time, a highly perturbed
input configuration bearing no resemblance to any of the three ideal reference
configurations is unlikely to be meaningfully denoised resulting in unknown or
divergent values of predicted noise. The case studies considered in this work
all suggest that our denoising model does not significantly alter such
disordered structures, including melt, point defects, dislocations, and grain
boundaries. This property is instrumental in allowing the denoising model to
reveal underlying crystalline order without impacting meaningful disordered
features in thermally perturbed configurations.

The score function plays a central role in modern likelihood-based generative
models such as the denoising diffusion probabilistic model (DDPM) [21, 22] and
score-based generative model [23], which can be unified under the same
framework [23]. Among its numerous recent achievements [24], DDPM has been
applied to crystal and molecular structure generations [30, 31]. In this work,
however, we apply the score matching method for denoising rather than
generative applications, and focus on a limited number of reference crystal
structures instead of many (thousands or millions) training images/structures.

###  4.2 Model training

Our (clean) data samples are reference crystal structures of interest
represented by the atomic coordinates $\mathbf{r}$ and the auxiliary
information $\mathbf{z}$:
$\mathbf{x}^{(i)}\rightarrow(\mathbf{r}^{(i)},\mathbf{z}^{(i)})$. The noise
prediction model $\epsilon_{\theta}$ was trained with entirely synthetic data,
which is generated by adding Gaussian noises to the atomic coordinates
$\mathbf{r}^{\prime}=\mathbf{r}+\sigma\epsilon$, with
$\sigma\sim\mathcal{U}(0,\sigma_{\max})$ drawn uniformly up to
$\sigma_{\max}\approx 13\%$ of the shortest interatomic distance, adhering to
Lindemann’s law on mean-squared thermal displacement of solids before melting
[32]. Since thermal displacements in classical MD are sums of phonon modes
that follow Boltzmann distributions rather than i.i.d. Gaussians, our working
hypothesis is that the correlated thermal fluctuations have fewer degrees of
freedom than i.i.d. Gaussians, and a model trained with the latter can
adequately handle the former. Our implemented loss function, slightly adjusted
from Eq. 7, is

$L=\mathbb{E}_{\mathbf{r},\mathbf{z},\sigma,\epsilon}\left[\|\sigma\epsilon-\epsilon_{\theta}(\mathbf{r}+\sigma\epsilon,\mathbf{z})\|^{2}\right].$
(8)

The idea of mixing training data with random noises is not new in either
general-purpose or scientific machine learning. Adding a small amount of noise
to the training data, known sometimes as the “noise trick”, is a well-
established data augmentation or regularization technique in general-purpose
machine learning to reduce overfitting and increase model robustness [33]. For
example, Zhou et al. [34, 35] adopted a hybrid training data pipeline of MD
trajectories and Gaussian noise displacements to fit the potential energy
surface of crystalline solids. A similar method was used by Chung et al. [15]
to identify ordered solid phases. The noise trick was also adopted to train
GNN surrogate models for physical simulations [36, 37]. This work, however,
makes denoising noise-corrupted inputs the centerpiece rather than merely a
regularization technique.

The model was trained with randomly drawn $\mathbf{x}\in\\{$FCC, BCC,
HCP$\\}$, using the AdamW optimizer [38] and a learning rate of $2\times
10^{-4}$, over 20,000 weight updates in minibatches of 32 samples. Each
FCC/HCP/BCC cell consists of roughly 1000 atoms. The training was carried out
using PyTorch [39] and PyTorch-Geometric [40]. All other training parameters,
if unspecified in this work, default to values per PyTorch 1.11.0 and PyTorch
Geometric 2.0.4.

###  4.3 Equivariant graph network model (NequIP)

The denoising model output is a vector (on each atom) that should be
equivariant under translation, rotation and mirror operations—the same
requirements for force fields or interatomic potentials. We adopted a
customized version of the E(3)-equivariant NequIP model [18, 19], which
guarantees such equivariance. NequIP is primarily built upon the idea of
equivariant tensor product between two inputs of irreducible representations,
or irreps, resulting in another irreps as the output. Unlike regular tensor
products, the tensor products in NequIP are parametrized by learnable weights
and are therefore termed WeightedTP in this work. Since the exact mathematical
details of the equivariant tensor product can be dense and complex, we refer
to the original work for their precise description [19]. We chose to directly
predict vector outputs (noising displacements) rather than a scalar output.

The main components of our NequIP variant consist of the initial embedding,
the interaction layers, and the final self-interaction layer to produce the
noise output (Fig. 8).

Figure 8: The NequIP model architecture.

In the initial embedding, the structure input is converted to an atomic graph,
with $\mathbf{h}_{i}$ as some attributes for node or atom $i$;
$\tilde{\mathbf{h}}_{i}$ as another set of attributes for the same node $i$;
and $\mathbf{e}_{ij}$ as the vector for the _directed_ edge from node $i$ to
node $j$. Transformed from atom type information by a trainable embedding
matrix, $\mathbf{h}_{i}$ initially only holds scalar information ($l=0$, where
$l$ is the tensor rank or the degree of representation) but is typically
expanded to hold information of higher tensor ranks ($l=1,2,...$) in
subsequent layers. $\tilde{\mathbf{h}}_{i}$, also transformed from atom type
information by a trainable embedding matrix, holds only scalar information and
does not change throughout the model layers. The embedding matrix is not
needed for the single-element systems but necessary for the SiO2 polymorphs.
We used Atomic Simulation Environment [41] and PyTorch-Geometric [40] for the
conversion to graphs.

Each interaction layer consists of several sub-operations: the self-
interaction, the convolution, SkipInit [42], and the gate activation. The
self-interaction updates the attributes of each node $\mathbf{h}_{i}$ via the
WeightedTP operation with $\tilde{\mathbf{h}}_{i}$ and does not aggregate
information from neighbor nodes:

$\mathbf{h}^{\prime}_{i}=\mathrm{WeightedTP}\left(\mathbf{h}_{i},\tilde{\mathbf{h}}_{i}\right)$
(9)

The convolution updates the attributes of each node $\mathbf{h}_{i}$ as the
sum of the WeightedTP operations between the neighboring nodes
$\mathbf{h}_{j}$ and the spherical harmonics of neighboring edges
$Y(\hat{\mathbf{e}}_{ij})$, with the weights learned from the edge distances
$\|\mathbf{e}_{ij}\|$ via a multilayer perceptron (MLP):

$\mathbf{h}^{\prime}_{i}=\frac{1}{Z}\sum_{j\in
N(i)}\mathrm{WeightedTP}_{\|\mathbf{e}_{ij}\|}(\mathbf{h}_{j},Y(\hat{\mathbf{e}}_{ij}))$
(10)

where $\hat{\mathbf{e}}_{ij}$ is the normalized version of $\mathbf{e}_{ij}$
and is therefore a unit vector pointing from node $i$ to node $j$, $N(i)$
denotes the neighbor nodes of node $i$, and $Z$ is a normalization constant.
The MLP contains one hidden layer. The initial layer of the MLP is the basis
function values expanded from the edge distance. For the SkipInit mechanism
[42], the scalar multipliers $\alpha$ are learned from yet another WeightedTP
operation between $\mathbf{h}_{i}$ and $\tilde{\mathbf{h}}_{i}$ (similar to
the self-interaction operation). The gate activation applies equivariant
nonlinearities [43] to the node attributes.

In the end, the final self-interaction layer serves to transform the node
attributes $\mathbf{h}_{i}^{(L-1)}$ from the second to last layer $L-1$, which
may hold scalar, vectorial, and tensorial features at node $i$, into a single
vector describing the noise output:

$\mathbf{h}^{(L)}_{i}=(\epsilon_{\theta})_{i}=\mathrm{WeightedTP}\left(\mathbf{h}^{(L-1)}_{i},\tilde{\mathbf{h}}^{(L-1)}_{i}\right),$
(11)

The complexity of the model is largely determined by the specified irreps
format for the node and edge attributes. For example, an array of 4 scalars
and 8 vectors can be written as 4x0e + 8x1o, with the numbers 4 and 8
describing the multiplicities, the numbers 0 and 1 describing the tensor rank,
and the letters e (even) and o (odd) describing the parity. Higher
multiplicities and tensor ranks can often result in better performance but
also larger memory and computational requirements. We intentionally kept the
model complexity small in favor of scalability to structures of millions of
atoms. The model settings are listed in Table 2.

Table 2: Model settings. Name | Value
---|---
Irreps for initial node attributes $\mathbf{h}$ | 10x0e
Irreps for auxiliary node attributes $\tilde{\mathbf{h}}$ | 10x0e
Irreps for intermediate/hidden node attributes $\mathbf{h}$ | 10x0e + 10x1e
Irreps for edge spherical harmonics $Y(\hat{\mathbf{e}}_{ij})$ | 1x0e + 1x1e + 1x2e
Number of basis functions for expanding edge distance $\|\mathbf{e}_{ij}\|$ | 16
Radius cutoff for edge distance $\|\mathbf{e}_{ij}\|$ | 3.2 Å
Number of hidden neurons for the MLP | 64
Normalization constant $Z$ | 12

###  4.4 Molecular dynamics simulations

This section describes the MD simulations used to demonstrate the capabilities
of our denoising method. These are simulations of (1) BCC, FCC and HCP Cu
structures, both defect-free and with point defects; (2) a solid crystal
nucleus growing inside melted Cu; (3) crystal plasticity in single crystal of
Ta in the BCC phase; (4) Ta grain boundaries; and (5) SiO2 in $\beta$-quartz
and cristobalite polymorphs. All simulations were performed with periodic
boundary conditions using LAMMPS [44].

The MD simulations for Cu using the embedded-atom method (EAM) potential by
Mishin et al. [45] were performed in BCC, FCC and HCP cells containing 1024,
1372, and 1152 atoms respectively. Since the bulk BCC phase is dynamically
unstable in Cu with imaginary phonon frequencies at ambient conditions, the
calculations were performed at a pressure of 60 GPa, where the BCC phase
becomes metastable.

Melting points of the three phases at 60 GPa have been calculated to be 3030 K
for BCC, 3066 K for HCP, and 3073 K for FCC [27]. Although FCC remains the
thermodynamically stable and thus preferred phase below 3073 K, free energies
of three phases are very close under these pressure and temperature
conditions. At slightly higher pressures (71.6 GPa and 85 GPa), the phase
diagram of the model of Cu contains triple points where two of the three solid
phases and the liquid phase coexist [27]. We have taken advantage of these
thermodynamic proximity of three crystal phases and the melt to set up an MD
simulation of a polymorphic critical solid nucleus simultaneously containing
all three solid phases surrounded by melt. The simulation contained 314,928 Cu
atoms and was initiated in an isobaric-isoenthalpic (NPH) ensemble at 70 GPa
from a small near-equilibrium nucleus with coexisting FCC and HCP ordered
regions containing about 200 FCC atoms and 300 HCP atoms, respectively. Upon
switching to an isobaric-isothermal (NVT) ensemble at the same pressure and
temperature 2800 K, the solid nucleus grows and partially transforms to the
BCC phase.

Interatomic interactions in tantalum were modeled using a well-known EAM
potential developed by Li et al. [46]. For simulations of crystal plasticity
in Ta described in Section 2.4, the crystals were created by arranging atoms
in a BCC lattice within in a cubic or an orthorhombic periodic supercell with
repeat vectors aligned along the cube axes of the BCC lattice. Dislocations
were seeded into the crystals in the form of one or several hexagon-shaped
prismatic loops of the vacancy type, following the procedure introduced in
[47]. For the configuration in Fig. 5a, a single dislocation loop was
introduced at the center of a cube-shaped simulation box made of 101,853 atoms
and annealed at temperature 2000 K. The complex network of dislocations shown
In Fig. 5b was generated by initially introducing 12 randomly positioned
dislocation loops into a $\sim 2$ million atoms box, annealing the model at
2500 K and zero pressure and then subjecting the crystal to uniaxial
compression along the [001] crystallographic axis at a “true” strain rate of
$2\times 10^{8}$/s for 2 ns while maintaining pressure near zero in an NPH
ensemble.

For Ta grain boundaries in Section 2.5, the periodic bi-crystal containing two
$\Sigma 5(310)$ symmetric tilt grain boundaries was created by joining two
crystal blocks of different lattice orientations obtained by rotating two
half-crystals in the opposite directions along the common $\langle 100\rangle$
tilt axis. The Ta polycrystal was assembled using atomsk from 12 randomly
seeded grains. Both the bi-crystal and the polycrystal were annealed at 2500 K
and zero pressure.

For silica, the Tersoff potential developed by Munetoh et al. [48] was
adopted. Simulations were performed in the NVT ensemble with unit cell
parameters taken from the Encyclopedia of Crystallographic Prototypes [49, 50,
51], and with temperature ramping from 600 K to 2000 K in steps of 200 K every
10 ps. A 5x5x5, 5x5x5, and 4x4x4 supercells were used for the $\beta$-quartz
(space group 181), $\alpha$-cristobalite (92), and $\beta$-cristobalite (227),
respectively.

## Data Availability

All data required to reproduce this work can be requested from the
corresponding author(s).

## Code Availability

The source code, the trained denoiser models, and a simple demo for this work
are available at http://www.github.com/llnl/graphite.

## Acknowledgements

TH and FZ acknowledge support by the Critical Materials Institute, an Energy
Innovation Hub funded by the U.S. Department of Energy, Office of Energy
Efficiency and Renewable Energy, and Advanced Materials and Manufacturing
Technologies Office. BS, CP, NB and VB are partially supported by the
Laboratory Directed Research and Development (LDRD) program (22-ERD-016) at
Lawrence Livermore National Laboratory. This work was performed under the
auspices of the US Department of Energy by Lawrence Livermore National
Laboratory under contract No. DE-AC52-07NA27344. We thank Dr. R. Freitas for
providing the DC3 benchmark dataset and Dr. A. Stukowski for useful
discussions.

## Author Contributions

TH implemented the denoising graph network model and analyzed the denoised MD
configurations. BS, NB and CP generated the MD data to be denoised.

## Competing interests

The authors state that there is no conflict of interest.

## References

  * [1] Paul J. Steinhardt, David R. Nelson and Marco Ronchetti  “Bond-orientational order in liquids and glasses”  In _Physical Review B_ 28.2, 1983, pp. 784–805  DOI: 10.1103/PhysRevB.28.784
  * [2] J. Honeycutt and Hans C. Andersen  “Molecular dynamics study of melting and freezing of small Lennard-Jones clusters”  In _Journal of Physical Chemistry_ 91.19, 1987, pp. 4950–4963  DOI: 10.1021/j100303a014
  * [3] Alexander Stukowski  “Structure identification methods for atomistic simulations of crystalline materials”  In _Modelling and Simulation in Materials Science and Engineering_ 20.4, 2012, pp. 045021  DOI: 10.1088/0965-0393/20/4/045021
  * [4] Hajime Tanaka, Hua Tong, Rui Shi and John Russo  “Revealing key structural features hidden in liquids and glasses”  In _Nature Reviews Physics_ 1.5 Springer US, 2019, pp. 333–348  DOI: 10.1038/s42254-019-0053-3
  * [5] Wolfgang Lechner and Christoph Dellago  “Accurate determination of crystal structures based on averaged local bond order parameters”  In _Journal of Chemical Physics_ 129.11, 2008  DOI: 10.1063/1.2977970
  * [6] Cynthia L Kelchner, S.. Plimpton and J.. Hamilton  “Dislocation nucleation and defect structure during surface indentation”  In _Physical Review B_ 58.17, 1998, pp. 11085–11088  DOI: 10.1103/PhysRevB.58.11085
  * [7] F.. Sapozhnikov, G.. Ionov and V.. Dremov  “An adaptive template method for analyzing crystal structures and defects in molecular dynamics simulations of high-rate deformations”  In _Russian Journal of Physical Chemistry B_ 2.2, 2008, pp. 238–245  DOI: 10.1134/S1990793108020115
  * [8] Alexander Stukowski and Karsten Albe  “Extracting dislocations and non-dislocation crystal defects from atomistic simulation data”  In _Modelling and Simulation in Materials Science and Engineering_ 18.8, 2010  DOI: 10.1088/0965-0393/18/8/085001
  * [9] Alexander Stukowski, Vasily V. Bulatov and Athanasios Arsenlis  “Automated identification and indexing of dislocations in crystal interfaces”  In _Modelling and Simulation in Materials Science and Engineering_ 20.8, 2012  DOI: 10.1088/0965-0393/20/8/085007
  * [10] QHwan Kim, Joon Hyuk Ko, Sunghoon Kim and Wonho Jhe  “GCIceNet: A graph convolutional network for accurate classification of water phases”  In _Physical Chemistry Chemical Physics_ 22.45 Royal Society of Chemistry, 2020, pp. 26340–26350  DOI: 10.1039/d0cp03456h
  * [11] Kirk Swanson et al.  “Deep learning for automated classification and characterization of amorphous materials”  In _Soft Matter_ 16.2 Royal Society of Chemistry, 2020, pp. 435–446  DOI: 10.1039/c9sm01903k
  * [12] Hideo Doi, Kazuaki Z. Takahashi and Takeshi Aoyagi  “Mining of effective local order parameters for classifying crystal structures: A machine learning study”  In _The Journal of chemical physics_ 152.21 AIP Publishing, LLC, 2020, pp. 214501  DOI: 10.1063/5.0005228
  * [13] Sébastien Becker, Emilie Devijver, Rémi Molinier and Noël Jakse  “Unsupervised topological learning for identification of atomic structures”  In _Physical Review E_ 105.4, 2022, pp. 045304  DOI: 10.1103/PhysRevE.105.045304
  * [14] Andreas Leitherer, Angelo Ziletti and Luca M. Ghiringhelli  “Robust recognition and exploratory analysis of crystal structures via Bayesian deep learning”  In _Nature Communications_ 12.1 Springer US, 2021, pp. 1–13  DOI: 10.1038/s41467-021-26511-5
  * [15] Heejung W. Chung, Rodrigo Freitas, Gowoon Cheon and Evan J. Reed  “Data-centric framework for crystal structure identification in atomistic simulations using machine learning”  In _Physical Review Materials_ 6.4 American Physical Society, 2022, pp. 22–31  DOI: 10.1103/PhysRevMaterials.6.043801
  * [16] V.. Hernandes, M.. Marques and José Rafael Bordin  “Phase classification using neural networks: Application to supercooled, polymorphic core-softened mixtures”  In _Journal of Physics Condensed Matter_ 34.2, 2022  DOI: 10.1088/1361-648X/ac2f0f
  * [17] James Chapman, Nir Goldman and Brandon C. Wood  “Efficient and universal characterization of atomic structures through a topological graph order parameter”  In _npj Computational Materials_ 8.1 Springer US, 2022, pp. 1–12  DOI: 10.1038/s41524-022-00717-7
  * [18] Mario Geiger et al.  “e3nn/e3nn: 2021-05-04”  Zenodo, 2021  DOI: 10.5281/zenodo.4735637
  * [19] Simon Batzner et al.  “E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials”  In _Nature Communications_ 13.1, 2022, pp. 2453  DOI: 10.1038/s41467-022-29939-5
  * [20] Pascal Vincent  “A connection between scorematching and denoising autoencoders”  In _Neural Computation_ 23.7, 2011, pp. 1661–1674  DOI: 10.1162/NECO_a_00142
  * [21] Jascha Sohl-Dickstein, Eric A. Weiss, Niru Maheswaranathan and Surya Ganguli  “Deep Unsupervised Learning using Nonequilibrium Thermodynamics”  In _Icml_ , 2015  arXiv: http://arxiv.org/abs/1503.03585
  * [22] Jonathan Ho, Ajay Jain and Pieter Abbeel  “Denoising diffusion probabilistic models”  In _Advances in Neural Information Processing Systems_ 2020-Decem.NeurIPS 2020, 2020, pp. 1–25  arXiv:2006.11239
  * [23] Yang Song et al.  “Score-Based Generative Modeling through Stochastic Differential Equations”, 2020, pp. 1–36  arXiv: http://arxiv.org/abs/2011.13456
  * [24] Ling Yang et al.  “Diffusion Models: A Comprehensive Survey of Methods and Applications”  Association for Computing Machinery, 2022  arXiv: http://arxiv.org/abs/2209.00796
  * [25] F.H. Stillinger and D.K. Stillinger  “Expanded solid matter: Two-dimensional LJ modeling” Advances in Disordered Materials  In _Mechanics of Materials_ 38.8, 2006, pp. 958–968  DOI: https://doi.org/10.1016/j.mechmat.2005.06.021
  * [26] Alexander Stukowski  “Visualization and analysis of atomistic simulation data with OVITO–the Open Visualization Tool”  In _Modelling and simulation in materials science and engineering_ 18.1 IOP Publishing, 2009, pp. 015012
  * [27] Babak Sadigh, Luis Zepeda-Ruiz and Jonathan L Belof  “Metastable–solid phase diagrams derived from polymorphic solidification kinetics”  In _Proceedings of the National Academy of Sciences_ 118.9, 2021, pp. e2017809118  DOI: 10.1073/pnas.2017809118
  * [28] J.. Stimac, N. Bertin, J.. Mason and V.. Bulatov  “Energy storage under high-rate compression of single crystal tantalum”  In _Acta Materialia_ 239 Elsevier Ltd, 2022, pp. 118253  DOI: 10.1016/j.actamat.2022.118253
  * [29] Aapo Hyvärinen  “Estimation of Non-Normalized Statistical Models by Score Matching”  In _Journal of Machine Learning Research_ 6, 2005, pp. 695–709  DOI: 10.5555/1046920.1088696
  * [30] Tian Xie  “Crystal Diffusion Variational Autoencoder for Periodic Material Generation”  In _ICLR_ , 2022, pp. 18–19
  * [31] Minkai Xu et al.  “GeoDiff: a Geometric Diffusion Model for Molecular Conformation Generation”, 2022, pp. 1–11  arXiv: http://arxiv.org/abs/2203.02923
  * [32] F.. Lindemann  “The calculation of molecular Eigen-frequencies”  In _Phys. Z._ 11, 1910, pp. 609–612
  * [33] Pascal Vincent et al.  “Stacked Denoising Autoencoders: Learning Useful Representations in a Deep Network with a Local Denoising Criterion”  In _Journal of Machine Learning Research_ 11.110, 2010, pp. 3371–3408  URL: http://jmlr.org/papers/v11/vincent10a.html
  * [34] Fei Zhou et al.  “Lattice Anharmonicity and Thermal Conductivity from Compressive Sensingof First-Principles Calculations”  In _Physical Review Letters_ 113.18 APS, 2014, pp. 185501  DOI: 10.1103/PhysRevLett.113.185501
  * [35] Fei Zhou, Weston Nielson, Yi Xia and Vidvuds Ozoliņš  “Compressive sensing lattice dynamics. I. General formalism”  In _Physical Review B_ 100.18 APS, 2019, pp. 184308  DOI: 10.1103/PhysRevB.100.184308
  * [36] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez and Peter W. Battaglia  “Learning mesh-based simulation with graph networks”  In _arXiv_ , 2020, pp. 1–18  arXiv:2010.03409
  * [37] Alvaro Sanchez-Gonzalez et al.  “Learning to Simulate Complex Physics with Graph Networks”  In _arXiv_ , 2020  arXiv:2002.09405
  * [38] Ilya Loshchilov and Frank Hutter  “Decoupled weight decay regularization”  In _arXiv preprint arXiv:1711.05101_ , 2017
  * [39] Adam Paszke et al.  “Pytorch: An imperative style, high-performance deep learning library”  In _Advances in neural information processing systems_ 32, 2019, pp. 8026–8037
  * [40] Matthias Fey and Jan Eric Lenssen  “Fast graph representation learning with PyTorch Geometric” Preprint at https://arxiv.org/abs/1903.02428, 2019
  * [41] Ask Hjorth Larsen et al.  “The atomic simulation environment—a Python library for working with atoms”  In _J Phys Condens Matter_ 29.27 IOP Publishing, 2017, pp. 273002
  * [42] Soham De and Sam Smith  “Batch normalization biases residual blocks towards the identity function in deep networks”  In _Advances in Neural Information Processing Systems_ 33, 2020, pp. 19964–19975
  * [43] Maurice Weiler et al.  “3d steerable cnns: Learning rotationally equivariant features in volumetric data”  In _Advances in Neural Information Processing Systems_ 31, 2018
  * [44] Steve Plimpton  “Fast parallel algorithms for short-range molecular dynamics”  In _Journal of computational physics_ 117.1 Elsevier, 1995, pp. 1–19
  * [45] Y. Mishin et al.  “Structural stability and lattice defects in copper: Ab initio, tight-binding, and embedded-atom calculations”  In _Phys. Rev. B_ 63 American Physical Society, 2001, pp. 224106  DOI: 10.1103/PhysRevB.63.224106
  * [46] Youhong Li, Donald J. Siegel, James B. Adams and Xiang-Yang Liu  “Embedded-atom-method tantalum potential developed by the force-matching method”  In _Phys. Rev. B_ 67, 2003, pp. 125101
  * [47] Luis A Zepeda-Ruiz, Alexander Stukowski, Tomas Oppelstrup and Vasily V Bulatov  “Probing the limits of metal plasticity with molecular dynamics simulations”  In _Nature_ 550.7677 Nature Publishing Group, 2017, pp. 492
  * [48] Shinji Munetoh, Teruaki Motooka, Koji Moriguchi and Akira Shintani  “Interatomic potential for Si–O systems using Tersoff parameterization”  In _Computational Materials Science_ 39.2 Elsevier, 2007, pp. 334–339
  * [49] Michael J Mehl et al.  “The AFLOW library of crystallographic prototypes: part 1”  In _Computational Materials Science_ 136 Elsevier, 2017, pp. S1–S828
  * [50] David Hicks et al.  “The AFLOW library of crystallographic prototypes: part 2”  In _Computational Materials Science_ 161 Elsevier, 2019, pp. S1–S1011
  * [51] David Hicks et al.  “The AFLOW library of crystallographic prototypes: part 3”  In _Computational Materials Science_ 199 Elsevier, 2021, pp. 110450
  * [52] P Ganesh and M Widom  “Signature of nearly icosahedral structures in liquid and supercooled liquid copper”  In _Physical Review B_ 74.13 APS, 2006, pp. 134205

##  Appendix A Visualization of toy score function

Here we briefly provide intuition and visualization for a score function based
on a toy data example. Suppose the data space is only two-dimensional. Let
there be three (noiseless) data points noised with Gaussians of different
standard deviations. Then the score function for the noised data is a gradient
field in the 2D space that converges on the three original data points
(Supplemental Fig. 1). Therefore, the act of denoising is equivalent to
iteratively following the gradient field that eventually points to the clean
data, which in this work is analogous to the ideal reference atomic motifs
such as BCC, FCC, and HCP.

Supplemental Figure 1: Visualization of a toy score function. Given a noisy
dataset based on some clean data samples, the score function for the noisy
data is a gradient field (indicated by the arrows) that converges to the
original noiseless data points. The score function shown here is not exact,
but was learned/estimated with a simple multilayer perceptron using PyTorch.

##  Appendix B The effect of denoising on melt/liquid Cu

Since the denoising model was trained to denoise perturbed solid structures
(namely, BCC, FCC, and HCP), how the model behaves with disordered melt/liquid
structure input is unknown. Here we investigate the effect of the denoising
process on the melt/liquid Cu phase.

Supplemental Figure 2: Histogram of atomic displacement magnitudes from
denoising BCC, FCC, HCP, and melt Cu. This displacement data is based on the
snapshots from Fig. 2, specifically the difference in atomic coordinates
between the original state (step 0) and the denoised state (step 8).

From denoising BCC, FCC, HCP, and melt/liquid Cu shown in Fig. 2, the atomic
displacements of the melt phase is similar to those of the solids, indicating
that all atoms move over roughly the same distance range regardless of the
phase (Supplemental Fig. 2). This is an interesting observation, as the melt
phase remains largely disordered after denoising according to the a-CNA
algorithm and the measured RDF (Fig. 2), even after the same level of
displacement as that of the solid phases. However, the melt Cu structure in
Fig. 2 is relatively small and may not provide enough statistics for a
sufficiently smooth RDF. We therefore applied denoising to a larger cell of
melt Cu phase (256,000 atoms) and measured the RDF, as shown in Fig. 3.
Similar to the result shown in Fig. 2c, the general RDF pattern after
denoising is similar to that of the original, with peak splitting in the short
range and more pronounced ordering in the medium range. We initially
hypothesized that the short range peaks would match those of ideal BCC, FCC,
and HCP structures. However, from Supplemental Fig. 3c, it is unclear whether
these short range peaks closely match those of the solid phases, with some
peaks that seemingly overlap with those of the solid phases and some peaks
that clearly do not.

Supplemental Figure 3: Comparison of the RDF of a denoised melt Cu
configuration with that of the original structure and those of BCC, FCC, and
HCP Cu. (a) Visualization of the original and denoised configurations where
the atoms are colored by the a-CNA prediction. (b) Comparison between the
original and denoised configurations. (c) Comparison between the denoised RDF
of the melt phase and those of the solid phases (BCC, FCC, and HCP).

Therefore, how the denoising model topologically impacts the disordered
liquid/melt phase is not yet fully understood. We posit that the emergence of
the short range peaks and the medium range ordering are attributed to the
atoms being moved into equidistant positions from the nearest neighbors as a
result of the model learning to form ideal atomic motifs. In other words, the
model may be aiming to denoise liquid structures into ideal lattices, but only
to a very limited extent. It is also possible that the melt phase, though
highly disordered, is associated with some inherent structure or order [52],
which in turn influences the behavior of the denoising model. In any case, a
more thorough investigation would be necessary to uncover the precise effect
of denoising on melt/liquid structures. However, such an inquiry is outside of
the scope of this work and is left as another dedicated future effort related
to the inherent order and symmetry in liquid.

##  Appendix C Denoised trajectory of BCC Ta plastic deformation

The plastic deformation trajectory prior to the snapshot shown in Fig. 5b is
also denoised and analyzed here. Using the same analysis for Fig. 5b, the
dislocations and the debris (point defects) are elucidated. This analysis
confirms the movements of the dislocations, as well as the corresponding
increase in the debris concentration (Supplemental Fig. 4), which is known to
increase under plastic deformation. Therefore, our denoising model is shown to
help capture complex networks of dislocations and point defects (or debris) in
perturbed and deformed structures that are otherwise difficult to analyze.
Indeed, analysis of point defects evolution in original (noisy) structures
generally cannot yield meaningful results as a large number of atoms are mis-
classified as defective only due to thermal displacements.

Supplemental Figure 4: Denoising helps reveal that debris (point defects)
concentration increases during plastic deformation as a result of moving
dislocations. (a) The dislocation and the debris are elucidated via the same
analysis applied to Fig. 5b in the main text, namely, with denoising followed
by cluster analysis of non-BCC atoms. (b) Quantitative measurement of the
debris concentration as a function of strain. The debris concentration is
measured as the number of non-BCC clusters (excluding the dislocations)
divided by the unit cell volume. Whether an atom belongs to BCC or other
phases is determined by a-CNA prediction.

##  Appendix D Denoising vs. vibration-averaging

Besides energy minimization, vibration-averaging is another common technique
for reducing thermal perturbations. Here we briefly provide a comparison
analysis between denoising and averaging for the plastic deformation
trajectory mentioned in Fig. 5b, with an emphasis on how denoising/averaging
impacts the frame-by-frame dislocation topological changes. Evolution of
dislocation networks in terms of local topological events can provide key
insights into fundamental dislocation mechanisms at the atomic scale, but
their detailed analysis is made extremely difficult and generally muddled by
incessant, spurious flickers caused by thermal vibrations. Focusing on the
subtle changes in the dislocation network topology, 100 frames of the plastic
deformation trajectory were simulated over small time steps of 0.01 ps
(Supplemental Fig. 5). For a given frame, the vibration-averaging was done by
averaging the atomic coordinates from the past 100 frames (including the
current frame).

Supplemental Figure 5: Comparison between the effects of denoising and
vibration-averaging on the frame-by-frame dislocation network topological
changes based on a short segment of the plastic deformation trajectory studied
in this work. (a) The trajectory made of 100 frames (with 0.01 ps timesteps)
with the atoms colored by a-CNA prediction. (b) The dislocation networks at
the initial frame, extracted for the original, denoised, and averaged
trajectories. (c) The frame-by-frame difference in the number of network links
along the (original, denoised, and averaged) trajectory. Since the topological
changes are hard to visually discern, the difference in the number of network
links is measured and plotted.

Visually, dislocation configurations extracted with DXA from the trajectories
display numerous visible but subtle changes in their topology in going from
one frame to the next. To quantify these changes, we count the frame-by-frame
difference in the number of dislocation links (dislocation lines connecting
two network nodes, i.e. edges of the dislocation network) along the trajectory
(Supplemental Fig. 5c). For the original (noisy) configurations, this number
has a standard deviation of 4.52. However, it is unclear which changes
correspond to “true” topological events and which are flickers induced by
thermal fluctuations. In vibration-averaged configurations, the corresponding
standard deviation is greatly reduced (1.6). In contrast, the difference in
the number of links varies by a standard deviation of 3.0 across for the
denoised configurations. This higher value may suggest that the vibration-
averaging technique is overzealous and eliminate true topological events as a
consequence of slightly smearing out the dynamics. Our denoising approach is
immune to such dynamic distortions. Thus, this example illustrates how
denoising can serve as an accurate pre-processing filter for methods of
structural analysis such as DXA.

##  Appendix E Denoising the DC3 benchmark dataset

The denoiser was applied to the public benchmark dataset introduced by Chung
et al. [15] for crystal structure identification. Structure types not
belonging to {BCC, FCC, HCP} are omitted from the dataset, and only the
configurations simulated at the melting point ($T=T_{m}$) were considered. 8
iterations were used for denoising. The results are shown in Supplemental Fig.
6.

Supplemental Figure 6: Findings from denoising and classifying the DC3
benchmark dataset. (a) CNA classification accuracy along the denoising steps.
(b) Based on the denoised structures (after 8 steps), there is some trace
amount of Frenkel defect pairs in BCC Fe. One such pair (circled) is shown
here. In (b), the atoms labeled as BCC by CNA are rendered semi-transparent.