We explore the possibility for reconstruction of the generative physical models describing interactions between atomic units in solids from observational electron microscopy data. Here, scanning transmission electron microscopy (STEM) is used to observe the dynamic motion of Si atoms at the edge of monolayer graphene under continuous electron beam illumination. The resulting time-lapsed STEM images represent the snapshots of observed chemical states of the system. We use two approaches: potential of mean force calculation using a radial distribution function and a direct fitting of the graphene–Si interatomic pairwise potentials with force matching, to reconstruct the force fields in the materials. These studies lay the foundation for quantitative analysis of materials energetics from STEM data through the sampling of the metastable states in the chemical space of the system.

Introduction of aberration correction in Scanning Transmission Electron Microscopy (STEM) has propelled this method to a technique of choice for characterizing a broad range of materials including metals and semiconductors, oxides, 2D materials, and many others.1–4 Beyond multiple qualitative studies visualizing structures of interfaces and localized and extended defects, recent improvements in spatial resolution and the stability of STEM instrumentation have enabled quantitative measurements of atom positions with picometer precision,5 with further advances enabled by segmented detectors and 4D STEM such as sub-picometer precision strain mapping,6 sub-angstrom charge-density mapping,7 electron ptychography for improved spatial resolution,8 and differential phase-contrast imaging.9 This quantitative structural STEM imaging now allows insight into the material’s structures and properties that previously was achievable only on average via diffraction-based methods, with the added advantage of element-, site-, and even charge/valence sensitivity.10 

This rapid growth in quantitative STEM has opened fundamentally new opportunities for exploring physics and chemistry of materials based on local structural measurements. Jia (via TEM)11–13 and Chisholm (via STEM)14 following the earlier work of Pan15 demonstrated direct probing of the symmetry breaking distortions in perovskites, enabling direct mapping of the polarization order parameter fields. This approach was further extended by Jia et al.16 and He et al.17 to probe octahedra tilts in the image plane and, via column shape analysis,18 in the beam directions. These studies have provided insight into ferroelectricity and screening phenomena at the oxide interfaces, emergence of topological defects in ferroelectric superlattices, etc. Mapping of the order parameter fields has now become de facto standard in the field.19–24 

Mapping of the mesoscopic order parameter fields can further be extended to extract the specific aspects of materials physics via matching with the Ginzburg–Landau type models. In one such approach, the analytic solutions for well-defined geometries with unknown boundary and gradient terms can be fitted to the atomically resolved fields to yield the value of the relevant material’s constants.25,26 Alternatively, the gradient coupling terms such as flexoelectric constant can be determined from theory–experiment matching.27 

However, mesoscopic order parameter fields can be defined only for the systems with continuous lattices and a small number of point and extended defects. Correspondingly, the potential of high-resolution STEM data to provide insight into the physics and chemistry of the systems with significant chemical disorder remains largely unexplored. Recently, it was proposed that such data can be mined to define a correlative picture of a material structure, for example, libraries of structural units and defects.28,29 This information can be further used to build the generative physical models. For lattice models, Vlcek et al.30–33 have demonstrated an approach to extract the interaction Hamiltonians from the images of chemically disordered systems using a statistical distance minimization method. However, until now, there have been no attempts to analyze the generative models behind the chemical transformations in STEM associated with the change of chemical bonding patterns.

Here, we demonstrate the analysis of dynamic STEM data to extract interaction potentials from the observations of atomic units. We use dynamic STEM data to acquire multiple snapshots of the dynamic process of chemical evolution in the Si–graphene system and use these to reconstruct the possible interaction potentials.

The graphene samples were fabricated using a wet transfer from the Cu growth substrate to a STEM sample grid and an oxygen baking procedure for cleaning, as detailed elsewhere.34 The dynamic observations were performed in a Nion UltraSTEM US200 operating at 100 kV accelerating voltage. The beam current was set to 10 pA and images were acquired with the high-angle annular dark-field (HAADF) detector to preserve Z-contrast.35–37 In this imaging mode, heavier atoms appear brighter allowing for a straightforward interpretation of image contrast and identification of all atoms. An area of the graphene sample was found where some contaminant material, comprised of mostly amorphous C and Si atoms, was next to an atomically pristine portion of the graphene lattice. An image sequence was then recorded over ∼56 min, capturing the sample evolution under the influence of the 100 kV e-beam in a 200 frame data stack (64 μs pixel dwell time, 31 pm pixel size, 16 nm field of view, and 17 s per frame). Figure 1 shows a selection of the acquired frames with the total accumulated dose indicated in the upper right (calculated as beam current * frame time * # of frames / frame area). We note that the atomic processes involved in the observed evolution occur much faster than can be directly captured in each frame. We observe a significant doping of graphene with the bright Si atoms at the edge of the contaminated region, consistent with previous observations of e-beam driven doping.38 This doping and possible chemical interaction with other species in the contamination (e.g., O) weakens the bonding of the C atoms and promotes structural degradation in this region. We observe a complete separation open between the contaminated region and the graphene lattice. As atoms are continually sputtered from the lattice, we also observe the appearance and movement of defects as well as the retraction of the graphene edge as the total number of C atoms decreases.

FIG. 1.

HAADF-STEM time series of graphene degrading under electron irradiation at 100 kV. (a) Initial state where the mostly pristine graphene lattice can be seen to be surrounded by an area of contamination consisting of primarily amorphous C and Si atoms. (b)–(f) Structural evolution of the graphene lattice during scanning of the electron beam at time steps 1428 s (b), 1819 s (c), 2108 s (d), 2499 s (e), and 3400 s (f). The accumulated total electron dose at each point in time is indicated in the upper right corner of each image.

FIG. 1.

HAADF-STEM time series of graphene degrading under electron irradiation at 100 kV. (a) Initial state where the mostly pristine graphene lattice can be seen to be surrounded by an area of contamination consisting of primarily amorphous C and Si atoms. (b)–(f) Structural evolution of the graphene lattice during scanning of the electron beam at time steps 1428 s (b), 1819 s (c), 2108 s (d), 2499 s (e), and 3400 s (f). The accumulated total electron dose at each point in time is indicated in the upper right corner of each image.

Close modal

The e-beam directed atomic evolution process is non-uniform in time. The graphene lattice is more stable at the beginning (i.e., pristine graphene is much more robust) and transformation from (a) to (b) proceeds over 84 frames (out of 200). The lattice degradation is apparently accelerated at the interface with the contamination, which could be due to chemical reactions with contaminant atoms driven by the e-beam. Finally, the most interesting phenomena are observed in the bulk of the graphene region. During continuous electron beam exposure, graphene does not nucleate multi-atomic voids (holes). Rather, a continuous generation of point defects such as Stone–Wales defects (5–7 rings), etc., is observed. Once formed, the defects apparently migrate to the edge, while preserving the integrity of the bulk. The edge recession observed from (c) to (f) is therefore driven by both defect migration and edge sputtering. This in turn suggests that the majority of the atomic dynamics observed during STEM imaging corresponds to the reversible transitions between possible defect states in the material, and the process is almost ergodic. In other words, the electron beam activates the Si–graphene system allowing it to relax into possible metastable states, and in such a way effectively samples the allowed chemical states of the system which can be energetically reached through the impartation of the beam energy.

We further note that both the transit time of the electron through the graphene sheet (∼attoseconds) and the lifetimes of the electronic excitations in graphene are well below the time scale of individual images. Therefore, the excitations that produce chemical changes in the system are essentially a delta function in time, and observed dynamics corresponds to the chemical relaxation processes. Therefore, the observed atomic configurations correspond to the snapshots of the possible chemical states of the graphene–Si system.

To identify positions of all atoms in every frame of the experimental movie, we utilized a deep fully convolutional neural network (FCNN). The application of FCNN to atom-resolved STEM data is described elsewhere.39 Briefly, the network was trained using simulated data to remove noise from raw experimental images and to separate image pixels associated with different atomic species into different classes. The structures used for model training were obtained from first-principles calculations of the relaxed geometries from the extended library of atomic defects in graphene under different strain.28 The relaxed coordinates were then used as an input into the MultiSlice algorithm40 to produce the simulated STEM images, which were corrupted by different types and levels of noise. The model training was done by utilizing Google's Tensor Processing Units. The FCNN output is a set of well-defined circular features on a uniform background (see the supplementary material). To convert the FCNN output into atomic coordinates, we used a simple center of mass measurement on the extracted features.

These data hence provide the information on possible configurations in the Si–graphene system. Due to the extremely rapid character of electron beam induced processes and associated energy relaxation in graphene and low yield of electron beam introduced changes, these images represent multiple realizations of the metastable atomic configurations in the chemical space in the system. Each of these in turn corresponds to the zero force acting on individual atoms. Here, we explore whether this information can be used to reconstruct the force fields acting between atoms.

As a context to these studies, there have been previous computational studies aimed at studying interactions in graphene. A study by Inui and Iwasaki investigated the interaction energy between a graphene sheet and a silicon substrate.41 The pairwise nonbonded interactions, modeled as the Lennard–Jones (LJ) potential, were σ=3.629Å and ϵ=8.91meV for the C–Si and the C–C parameters are given in Table I. In a LJ potential, σ refers to the distance at which the pairwise potential is 0 and ϵ refers to the depth of the potential well indicating the maximum strength of the attractive interaction. We also include the list of LJ parameters for non-bonded carbon atom interactions used for graphene modeling as reported in a review.42 The Lennard–Jones potential is perhaps the crudest assumption possible, especially since it cannot reproduce the lattice structure of graphene, but it can provide a starting point for comparison.

TABLE I.

Lennard–Jones parameters for pairwise interactions between carbon atoms in graphene sheet. Some parameters taken from Pykal et al.42 

Referenceσ (Å)ε (meV)
Inui and Iwasaki41  3.431 4.55 
Parm 9943  3.399 67 3.7 
OPLS44  3.550 0 3.0 
CHARMM2745  3.550 5 3.0 
Ulbricht et al.46  3.781 08 2.63 
Girifalco et al.47  3.412 14 2.39 
Cheng and Steele48  3.399 67 2.42 
COMPASS49  3.487 87 2.9 
Referenceσ (Å)ε (meV)
Inui and Iwasaki41  3.431 4.55 
Parm 9943  3.399 67 3.7 
OPLS44  3.550 0 3.0 
CHARMM2745  3.550 5 3.0 
Ulbricht et al.46  3.781 08 2.63 
Girifalco et al.47  3.412 14 2.39 
Cheng and Steele48  3.399 67 2.42 
COMPASS49  3.487 87 2.9 

Besides non-bonded interaction energies, different reactive energies in graphene, like bond-breaking, vacancy formation, diffusion, and merging and deformation energies, have been previously reported. Sun et al.50 studied C 1 s binding energy of carbon in different materials using DFT and reported C 1 s binding energy to be 284.80 eV for interior atoms of graphene. Machine learning,51ab initio calculations,51 and quantum based Morse potentials52 have also been used to estimate the aforementioned reactive energies of graphene carbon. Based on previous studies of graphene, dissociation energy was estimated to be 8.34 eV,52 energies for graphene lattice vacancy diffusion and merging were reported to be in the range of 1.1–1.3 eV53 and 1.2–2.1 eV53 respectively. Table 1 of Yang et al. gives a summary of various methods for computing vacancy formation.54 

Here, we used two methods to calculate the pairwise potential directly from the experimental video. One method was based on the structure and the other based on velocities. The first method used pairwise distances from video that were used with kernel density estimation (KDE). Gaussian kernels with a bandwidth of 0.2 were used. KDE gives a probability of pairwise distance, P(r). The radial distribution function (RDF), g(r) was subsequently calculated according to the equation C1P(r)r2, where C1 is a constant determined to ensure that g(r) is unity at large r. g(r) was then used to calculate the two particle PMF, w(2), using the equation

w(2)=kTlng(r)C2.
(1)

This approach of deriving the PMF from the RDF is valid only when we have statistical sampling from an equilibrium ensemble.55 Since the video from the STEM experiment shows ablation and gradual decrease in the particle number, the canonical restriction is not necessarily satisfied here. Nonetheless, we argue that the character of the process suggests significant stationary component, and hence this approach establishes a useful baseline for the two-body potential.

The second approach is to assume Langevin dynamics and fit the potential directly. This requires no assumptions about equilibrium. We assumed a model of pairwise potential using a basis set [B(r)] of 32 Gaussian functions along with a repulsive term, U(r). U(r) is given by

U(r)=u(rr0)(rr0)12,
(2)

where u(r) is a unit step function. The basis set is of the form

B(r)=hiG(si,wi).11+r,
(3)

where hi is the scaled height, si is the mean, and wi is the standard deviation of the ith Gaussian basis set function. We use overdamped Langevin dynamics to calculate the forces directly from the coordinates. The generalized Langevin equation56 is given by

mv˙=Ur+ξ(t)γv,
(4)

where γv is the dampening force and ξ(t) is the random noise term. U is the potential energy function, r is the distance, and v is the velocity. γ is the friction coefficient and ξ(t) is a Gaussian with zero mean. The velocity correlation decay time is given by τ=m/γ. For times, tτ, the inertia term, mv˙, can be neglected. This case is referred to as the overdamped Langevin.57 The corresponding equation is given by

v=γ1Ur+γ1ξ(t).
(5)

We have chosen the overdamped Langevin model since the time lapse between two consecutive samples from the experimental video is 17 s, a significant time period compared to the scale of atomic motion. Figure S1 in the supplementary material gives a quantitative evidence of the velocity autocorrelation function calculated from the experimental video. It shows rapid decay of the autocorrelation between τ=0 and τ=1 supporting the assumption that the correlation time is smaller than the time difference between consecutive frames. This shows that, like the overdamped Langevin model, there is no momentum correlation between frames. Besides the forces arising from pairwise interactions, all the contributions from other dynamical processes are incorporated into the random force term in Eq. (5). Using initial parameter values for B(r), Eq. (5) and a velocity-Verlet integrator with unit time step, we generate the predicted coordinates of the atoms in the next time step. The squared difference between the predicted coordinates (y) and the true coordinates (y) from the next frame of the video is the objective function to be minimized. Since the number of particles in each frame is not constant in the experimental video, we used an index reassignment strategy to trace each atom along the trajectory. To account for the image boundary effects and variable number of particles across the trajectory, we assigned a weight (θi), which is set to either 1 or 0, to each atom. The weight determines if an atom is considered during training. To ensure that all the frames are equally weighted during parameter training the cumulative weight (θ) for each frame was scaled to a uniform value.

The loss function which is minimized during training is given by

χ=(yy)2θ+1nBiNhi2r0,
(6)

where nB is the number of Gaussian functions in the basis set (32). This optimizes the potential function parameters. A similar method of learning pairwise potentials using Gaussian basis set functions has been previously reported.58 This is also like the use of a spline basis-set seen in force matching.59 

We validate our method by applying the discussed pair-potential regression strategy using a toy model with two types of particles, A and B, following overdamped Langevin dynamics of known parameters. The masses of A and B were set to those of C and Si, respectively. Lennard–Jones (LJ) potential was used as the potential in the Langevin equation for the toy-system. To generate the synthetic reference trajectory for the toy system, serving as a proxy for the experimental video, we used the first frame of the experimental video and initialized the positions of the particles for the first frame of the synthetic trajectory. We also tested multi-step integration between each step.60 These results can be found in the supplementary material.

All the calculations for both the methods, PMF calculation using RDF and pair-potential inference using overdamped Langevin dynamics, were based on the experimental video. The distance values had to be scaled to proper units. We used the first peak in the radial distribution function, corresponding to the C–C graphene bond, for scaling the distance to match the standard C–C graphene bond length of 0.142 nm. Figure 2 shows the RDFs for C–C and C–Si in graphene and the corresponding two particle PMFs, wCC(2) and wSiC(2), calculated using this approach.

FIG. 2.

RDFs and the corresponding PMFs for pairwise interactions in graphene. (a) C–C and C–Si RDFs in graphene generated from the processed experimental video. The RDFs were calculated using histograms of C–C and C–Si pairwise distances and fitting Gaussian KDE functions to the histogrammed distances. (b) Two particle PMFs, wCC(2) and wSiC(2), calculated from the RDFs obtained.

FIG. 2.

RDFs and the corresponding PMFs for pairwise interactions in graphene. (a) C–C and C–Si RDFs in graphene generated from the processed experimental video. The RDFs were calculated using histograms of C–C and C–Si pairwise distances and fitting Gaussian KDE functions to the histogrammed distances. (b) Two particle PMFs, wCC(2) and wSiC(2), calculated from the RDFs obtained.

Close modal

Figure 3 compares the two particle PMFs, wCC(2) and wSiC(2), and Si–C and C–C pair-potentials trained using the Langevin assumption. The locations of the first repulsive wells, indicated by the vertical lines in Fig. 3, are comparable. The later minima in the potential could be from lack of three-body terms. The Langevin fit is likely better, because it is based on a regression to the observed particle motions. The PMF, as mentioned above, should only agree if ergodic, equilibrium, NVT sampling is observed. These conditions do not hold, especially since the majority of statistics contributing are from the static graphene lattice. To obtain the well-depth in Fig. 3, we calibrated the energies according to the C–C Lennard–Jones parameters from Inui and Iwasaki (4.55 meV). This results in a C–Si minimum of 1.5 meV, showing that Si fits less favorably into the graphene lattice than carbon which is expected. Regarding the shape of the potentials, the Langevin method yields less extreme potentials but is still showing large peaks, likely because the graphene lattice structure is being projected onto this two-body potential.

FIG. 3.

C–C and Si–C PMFs and pairwise potentials as obtained by the two different methods. The result plotted in orange is the pairwise PMFs obtained using RDF and KDE. The pairwise potentials in blue were learned using Gaussian basis set functions.

FIG. 3.

C–C and Si–C PMFs and pairwise potentials as obtained by the two different methods. The result plotted in orange is the pairwise PMFs obtained using RDF and KDE. The pairwise potentials in blue were learned using Gaussian basis set functions.

Close modal

In this work, we have shown that the videos generated by dynamic STEM experiments can directly be used to reconstruct pairwise potentials. We have explored two techniques for reconstruction, one based only on structures observed and one based on the motion of the atoms. The use of motion (Langevin) requires less assumptions and gives more reasonable approximate two-body pairwise potentials. The structure-based approach (PMF) requires many assumptions and gives unreasonably large forces, mostly because a PMF analysis is inapplicable to the large fraction of static atoms. This work highlights the potential of using state-of-the-art characterization techniques, dynamic STEM, in conjunction with computational modeling methods to learn the underlying physics of different phenomena. This warrants further research along the line of the present work.

See the supplementary material for a detailed analysis of the velocity auto-correlation function and the effects of temperature and random force magnitude on regressed pair-potentials. The supplementary material also includes an example illustrating how the FCNN is used on noisy STEM data to obtain atomic coordinates.

This effort (electron microscopy, feature extraction) is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), Materials Sciences and Engineering Division (O.D., S.J., and S.V.K.) and was performed and partially supported (M.Z.) at the Center for Nanophase Materials Sciences, which is a DOE Office of Science User Facility. Theoretical analysis (M.C. and A.D.W.) is based upon work supported by the National Science Foundation (NSF) under Grant Nos. 1764415 and 1751471. The authors are grateful to Dr. R. Unocic for careful reading and commenting upon the manuscript.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
S. J.
Pennycook
 et al,
Advances in Imaging and Electron Physics
(
Elsevier Academic Press Inc.
,
2008
), Vol. 153, p. 327.
2.
S. J.
Pennycook
and
P. D.
Nellist
,
Scanning Transmission Electron Microscopy: Imaging and Analysis
(
Springer
,
New York
,
2011
).
3.
P. E.
Batson
,
N.
Dellby
, and
O. L.
Krivanek
, “
Sub-angstrom resolution using aberration corrected electron optics
,”
Nature
418
,
617
620
(
2002
).
4.
N.
Dellby
,
O. L.
Krivanek
,
P. D.
Nellist
,
P. E.
Batson
, and
A. R.
Lupini
, “
Progress in aberration-corrected scanning transmission electron microscopy
,”
J. Electron Microsc.
50
,
177
185
(
2001
).
5.
A. B.
Yankovich
 et al, “
Picometre-precision analysis of scanning transmission electron microscopy images of platinum nanocatalysts
,”
Nat. Commun.
5
,
4155
(
2014
).
6.
Y.
Han
 et al, “
Strain mapping of Two-dimensional heterostructures with subpicometer precision
,”
Nano Lett.
18
,
3746
3751
(
2018
).
7.
W.
Gao
 et al, “
Real-space charge-density imaging with sub-ångström resolution by four-dimensional electron microscopy
,”
Nature
575
,
480
484
(
2019
).
8.
Y.
Jiang
 et al, “
Electron ptychography of 2D materials to deep sub-angstrom resolution
,”
Nature
559
,
343
(
2018
).
9.
N.
Shibata
 et al, “
Differential phase-contrast microscopy at atomic resolution
,”
Nat. Phys.
8
,
611
615
(
2012
).
10.
M.
Varela
 et al, “
Spectroscopic imaging of single atoms within a bulk solid
,”
Phys. Rev. Lett.
92
,
095502
(
2004
).
11.
C. L.
Jia
 et al, “
Unit-cell scale mapping of ferroelectricity and tetragonality in epitaxial ultrathin ferroelectric films
,”
Nat. Mater.
6
,
64
69
(
2007
).
12.
C. L.
Jia
,
K. W.
Urban
,
M.
Alexe
,
D.
Hesse
, and
I.
Vrejoiu
, “
Direct observation of continuous electric dipole rotation in flux-closure domains in ferroelectric Pb(Zr,Ti)O(3)
,”
Science
331
,
1420
1423
(
2011
).
13.
C. L.
Jia
 et al, “
Effect of a single dislocation in a heterostructure layer on the local polarization of a ferroelectric layer
,”
Phys. Rev. Lett.
102
,
117601
(
2009
).
14.
M. F.
Chisholm
,
W. D.
Luo
,
M. P.
Oxley
,
S. T.
Pantelides
, and
H. N.
Lee
, “
Atomic-scale compensation phenomena at polar interfaces
,”
Phys. Rev. Lett.
105
,
197602
(
2010
).
15.
X. Q.
Pan
,
W. D.
Kaplan
,
M.
Ruhle
, and
R. E.
Newnham
, “
Quantitative comparison of transmission electron microscopy techniques for the study of localized ordering on a nanoscale
,”
J. Am. Ceram. Soc.
81
,
597
605
(
1998
).
16.
C. L.
Jia
 et al, “
Oxygen octahedron reconstruction in the SrTiO(3)/LaAlO(3) heterointerfaces investigated using aberration-corrected ultrahigh-resolution transmission electron microscopy
,”
Phys. Rev. B
79
,
081405(R)
(
2009
).
17.
A. Y.
Borisevich
 et al, “
Suppression of octahedral tilts and associated changes in electronic properties at epitaxial oxide heterostructure interfaces
,”
Phys. Rev. Lett.
105
,
087204
(
2010
).
18.
Q.
He
 et al, “
Towards 3D mapping of BO6 octahedron rotations at perovskite heterointerfaces, unit cell by unit cell
,”
ACS Nano
9
,
8412
8419
(
2015
).
19.
S.
Das
 et al, “
Observation of room-temperature polar skyrmions
,”
Nature
568
,
368
(
2019
).
20.
C. T.
Nelson
 et al, “
Spontaneous vortex nanodomain arrays at ferroelectric heterointerfaces
,”
Nano Lett.
11
,
828
834
(
2011
).
21.
L. Z.
Li
 et al, “
Observation of strong polarization enhancement in ferroelectric tunnel junctions
,”
Nano Lett.
19
,
6812
6818
(
2019
).
22.
M.
Campanini
,
R.
Erni
,
C. H.
Yang
,
R.
Ramesh
, and
M. D.
Rossell
, “
Periodic giant polarization gradients in doped BiFeO3 thin films
,”
Nano Lett.
18
,
717
724
(
2018
).
23.
G.
Van Tendeloo
,
S.
Bals
,
S.
Van Aert
,
J.
Verbeeck
, and
D.
Van Dyck
, “
Advanced electron microscopy for advanced materials
,”
Adv. Mater.
24
,
5655
5675
(
2012
).
24.
R. J.
Zeches
 et al, “
A strain-driven morphotropic phase boundary in BiFeO(3)
,”
Science
326
,
977
980
(
2009
).
25.
A. Y.
Borisevich
 et al, “
Exploring mesoscopic physics of vacancy-ordered systems through atomic scale observations of topological defects
,”
Phys. Rev. Lett.
109
,
065702
(
2012
).
26.
H. J.
Chang
 et al, “
Atomically resolved mapping of polarization and electric fields across ferroelectric/oxide interfaces by Z-contrast imaging
,”
Adv. Mater.
23
,
2474
(
2011
).
27.
Q.
Li
 et al, “
Quantification of flexoelectricity in PbTiO3/SrTiO3 superlattice polar vortices using machine learning and phase-field modeling
,”
Nat. Commun.
8
,
13936
(
2017
). .
28.
M.
Ziatdinov
 et al, “
Building and exploring libraries of atomic defects in graphene: Scanning transmission electron and scanning tunneling microscopy study
,”
Sci. Adv.
5
,
eaaw8989
(
2019
).
29.
A. G.
Rajan
 et al, “
Addressing the isomer cataloguing problem for nanopores in two-dimensional materials
,”
Nat. Mater.
18
,
129
(
2019
).
30.
L.
Vlcek
 et al, “
Learning from imperfections: Predicting structure and thermodynamics from atomic imaging of fluctuations
,”
ACS Nano
13
,
718
727
(
2019
).
31.
L.
Vlcek
,
A.
Maksov
,
M. H.
Pan
,
R. K.
Vasudevan
, and
S. V.
Kahnin
, “
Knowledge extraction from atomically resolved images
,”
ACS Nano
11
,
10313
10320
(
2017
).
32.
L.
Vlcek
,
W. W.
Sun
, and
P. R. C.
Kent
, “
Combining configurational energies and forces for molecular force field optimization
,”
J. Chem. Phys.
147
,
161713
(
2017
).
33.
L.
Vlcek
,
R. K.
Vasudevan
,
S.
Jesse
, and
S. V.
Kalinin
, “
Consistent integration of experimental and ab initio data into effective physical models
,”
J. Chem. Theory Comput.
13
,
5179
5194
(
2017
).
34.
O.
Dyck
,
S.
Kim
,
S. V.
Kalinin
, and
S.
Jesse
, “
Mitigating e-beam-induced hydrocarbon deposition on graphene for atomic-scale scanning transmission electron microscopy studies
,”
J. Vac. Sci. Technol. B
36
,
011801
(
2017
).
35.
S. J.
Pennycook
and
D. E.
Jenson
, “
High-resolution Z-contrast imaging of crystals
,”
Ultramicroscopy
37
,
14
38
(
1991
).
36.
D. E.
Jesson
,
S. J.
Pennycook
, and
J. M.
Baribeau
, “‘
Column-bycolumn’ compositional mapping at semiconductor interfaces using Z-contrast STEM
,”
High Resolut. Electron Microsc. Defects Mater.
183
,
223
230
(
1990
).
37.
S. J.
Pennycook
, “
Z-contrast stem for materials science
,”
Ultramicroscopy
30
,
58
69
(
1989
).
38.
O.
Dyck
,
S.
Kim
,
S. V.
Kalinin
, and
S.
Jesse
, “
Placing single atoms in graphene with a scanning transmission electron microscope
,”
Appl. Phys. Lett.
111
,
113104
(
2017
).
39.
M.
Ziatdinov
,
C.
Nelson
,
R. K.
Vasudevan
,
D. Y.
Chen
, and
S. V.
Kalinin
, “
Building ferroelectric from the bottom up: The machine learning analysis of the atomic-scale ferroelectric distortions
,”
Appl. Phys. Lett.
115
,
052902
(
2019
).
40.
J.
Barthel
, “
Dr. Probe: A software for high-resolution STEM image simulation
,”
Ultramicroscopy
193
,
1
11
(
2018
).
41.
N.
Inui
and
S.
Iwasaki
, “
Interaction energy between graphene and a silicon substrate using pairwise summation of the Lennard-Jones potential
,”
e-J. Surf. Sci. Nanotechnol.
15
,
40
49
(
2017
).
42.
M.
Pykal
,
P.
Jurečka
,
F.
Karlický
, and
M.
Otyepka
, “
Modelling of graphene functionalization
,”
Phys. Chem. Chem. Phys.
18
,
6351
6372
(
2016
).
43.
J.
Wang
,
P.
Cieplak
, and
P. A.
Kollman
, “
How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules?
,”
J. Comput. Chem.
21
,
1049
1074
(
2000
).
44.
W. L.
Jorgensen
,
D. S.
Maxwell
, and
J.
Tirado-Rives
, “
Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids
,”
J. Am. Ceram. Soc.
118
,
11225
11236
(
1996
).
45.
N.
Foloppe
,
J.
MacKerell
, and
D.
Alexander
, “
All-atom empirical force field for nucleic acids: I. Parameter optimization based on small molecule and condensed phase macromolecular target data
,”
J. Comput. Chem.
21
,
86
104
(
2000
).
46.
H.
Ulbricht
,
G.
Moos
, and
T.
Hertel
, “
Interaction of C60 with carbon nanotubes and graphite
,”
Phys. Rev. Lett.
90
,
095501
(
2003
).
47.
L. A.
Girifalco
,
M.
Hodak
, and
R. S.
Lee
, “
Carbon nanotubes, buckyballs, ropes, and a universal graphitic potential
,”
Phys. Rev. B
62
,
13104
13110
(
2000
).
48.
A.
Cheng
and
W. A.
Steele
, “
Computer simulation of ammonia on graphite. II. Monolayer melting
,”
J. Chem. Phys.
92
,
3867
3873
(
1990
).
49.
H.
Sun
, “
COMPASS:  An ab initio force-field optimized for condensed-phase: Applications overview with details on alkane and benzene compounds
,”
J. Phys. Chem. B
102
,
7338
7364
(
1998
).
50.
C. Q.
Sun
 et al, “
Coordination-resolved C−C bond length and the C 1 s binding energy of carbon allotropes and the effective atomic coordination of the few-layer graphene
,”
J. Phys. Chem. C
113
,
16464
16467
(
2009
).
51.
P.
Rowe
,
G.
Csányi
,
D.
Alfè
, and
A.
Michaelides
, “
Development of a machine learning potential for graphene
,”
Phys. Rev. B
97
,
054303
(
2018
).
52.
B. I.
Costescu
,
I. B.
Baldus
, and
F.
Gräter
, “
Graphene mechanics I. Efficient first principles based Morse potential
,”
Phys. Chem. Chem. Phys.
16
,
12591
12598
(
2014
).
53.
T.
Trevethan
,
C. D.
Latham
,
M. I.
Heggie
,
P. R.
Briddon
, and
M. J.
Rayson
, “
Vacancy diffusion and coalescence in graphene directed by defect strain fields
,”
Nanoscale
6
,
2978
2986
(
2014
).
54.
G.
Yang
,
L.
Li
,
W. B.
Lee
, and
M. C.
Ng
, “
Structure of graphene and its disorders: A review
,”
Sci. Technol. Adv. Mater.
19
,
613
648
(
2018
).
55.
D.
Chandler
, “
Introduction to modern statistical mechanics
,” in
Introduction to Modern Statistical Mechanics
, edited by
D.
Chandler,
Foreword by David Chandler, September 1987 (
Oxford University Press
,
1987
), pp.
288
, ISBN-10: 0195042778, ISBN-13: 9780195042771.
56.
S. A.
Adelman
, in
Advances in Chemical Physics
, edited by
I.
Prigogine
and
S. A.
Rice
(John Wiley & Sons, Ltd.,
2007
), Chap. 2, pp.
143
253
.
57.
G. A.
Pavliotis
,
Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations
(
Springer
,
2014
), Vol. 60.
58.
B.
Rainier
,
C.
Maghesree
,
A.
Dilnoza
,
G.
Heta
, and
W.
Andrew
, “
A GPU accelerated machine learning framework for molecular simulation: HOOMD blue with tensor flow
,”
chemRxiv
8019527.
59.
S.
Izvekov
and
G. A.
Voth
, “
A multiscale coarse-graining method for biomolecular systems
,”
J. Phys. Chem. B
109
,
2469
2473
(
2005
).
60.
A.
Sanchez-Gonzalez
,
V.
Bapst
,
K.
Cranmer
, and
P.
Battaglia
, “
Hamiltonian graph networks with ODE integrators
,” arXiv:1909.12790 (
2019
).

Supplementary Material