Realistic Atomistic Structure of Amorphous Silicon from Machine-Learning-Driven Molecular Dynamics

Amorphous silicon (a-Si) is a widely studied noncrystalline material, and yet the subtle details of its atomistic structure are still unclear. Here, we show that accurate structural models of a-Si can be obtained using a machine-learning-based interatomic potential. Our best a-Si network is obtained by simulated cooling from the melt at a rate of 10 K/s (that is, on the 10 ns time scale), contains less than 2% defects, and agrees with experiments regarding excess energies, diffraction data, and Si NMR chemical shifts. We show that this level of quality is impossible to achieve with faster quench simulations. We then generate a 4096-atom system that correctly reproduces the magnitude of the first sharp diffraction peak (FSDP) in the structure factor, achieving the closest agreement with experiments to date. Our study demonstrates the broader impact of machine-learning potentials for elucidating structures and properties of technologically important amorphous materials. A silicon (a-Si) is a fundamental and widely studied noncrystalline material, with applications ranging from photovoltaics and thin-film transistors to electrodes in batteries. Its atomic-scale structure is traditionally approximated in a Zachariasen-like picture with all atoms in locally “crystal-like”, tetrahedral environments, but without long-range order. However, the real material contains a nonzero amount of coordination defects, colloquially referred to as “dangling bonds” (under-coordinated sites) and “floating bonds” (overcoordinated sites). Knowing the properties and abundance of such defects is important, as they can control electronic and other macroscopic properties. We note at the outset that, although defect sites in a-Si may be passivated by hydrogenation (to give “a-Si:H”) in some synthetic conditions, we here focus on the archetypical, hydrogen-free material as made in ion-implantation or sputter-deposition experiments. Even the most advanced experimental approaches do not directly allow the observation of the bulk atomic structure in amorphous materials. Despite significant advances, including in situ NMR techniques and “inverse” approaches such as Reverse Monte Carlo (RMC) modeling of diffraction data, only indirect knowledge can be gained about the local atomic environments, and that only in a statistical sense. For almost three decades, molecular-dynamics (MD) simulations have therefore played a crucial and complementary role, with a-Si being a prominent example. These simulations either use density-functional theory (DFT) or classical force fields. DFTMD describes a system with quantum-mechanical accuracy and can largely correctly capture the structural and bonding subtleties of liquid and amorphous matter. However, it is computationally expensive, and therefore allows only limited system sizes (a few hundred atoms at most) and time scales to be simulated. Indeed, the cooling rates in previous DFT simulations of a-Si (≈ 10 K/s) are orders of magnitude faster than those in experiments. Classical force fields require much less computational effort, giving access to nanometerscale (“device-size”) structural models, both for a-Si and for multicomponent systems derived from it (see ref 25 for but one example). However, they are rarely accurate enough to fully describe the structural variations present in the amorphous state. Capitalizing on today’s “big-data” revolution, machinelearning (ML) algorithms are increasingly used to generate interatomic potentials for atomistic simulations. By “learning from” (or rather, f itting to) quantum-mechanically computed reference data for energies and forces, ML-based interatomic potentials can enable simulations with an accuracy Received: March 24, 2018 Accepted: May 12, 2018 Published: May 12, 2018 Letter pubs.acs.org/JPCL Cite This: J. Phys. Chem. Lett. 2018, 9, 2879−2885 © 2018 American Chemical Society 2879 DOI: 10.1021/acs.jpclett.8b00902 J. Phys. Chem. Lett. 2018, 9, 2879−2885 This is an open access article published under a Creative Commons Attribution (CC-BY) License, which permits unrestricted use, distribution and reproduction in any medium, provided the author and source are cited. D ow nl oa de d vi a 84 .6 8. 93 .4 0 on F eb ru ar y 5, 2 02 0 at 1 2: 37 :0 6 (U T C ). Se e ht tp s: //p ub s. ac s. or g/ sh ar in gg ui de lin es f or o pt io ns o n ho w to le gi tim at el y sh ar e pu bl is he d ar tic le s.

A morphous silicon (a-Si) is a fundamental and widely studied noncrystalline material, with applications ranging from photovoltaics and thin-film transistors to electrodes in batteries. 1−5 Its atomic-scale structure is traditionally approximated in a Zachariasen-like picture 6 with all atoms in locally "crystal-like", tetrahedral environments, but without long-range order. 7−9 However, the real material contains a nonzero amount of coordination defects, colloquially referred to as "dangling bonds" (under-coordinated sites) and "floating bonds" (overcoordinated sites). Knowing the properties and abundance of such defects is important, as they can control electronic and other macroscopic properties. We note at the outset that, although defect sites in a-Si may be passivated by hydrogenation (to give "a-Si:H") in some synthetic conditions, we here focus on the archetypical, hydrogen-free material as made in ion-implantation or sputter-deposition experiments. 10−14 Even the most advanced experimental approaches do not directly allow the observation of the bulk atomic structure in amorphous materials. Despite significant advances, including in situ NMR techniques 15,16 and "inverse" approaches such as Reverse Monte Carlo (RMC) modeling of diffraction data, 17−19 only indirect knowledge can be gained about the local atomic environments, and that only in a statistical sense. For almost three decades, molecular-dynamics (MD) simulations have therefore played a crucial and complementary role, with a-Si being a prominent example. 20−24 These simulations either use density-functional theory (DFT) or classical force fields. DFT-MD describes a system with quantum-mechanical accuracy and can largely correctly capture the structural and bonding subtleties of liquid and amorphous matter. However, it is computationally expensive, and therefore allows only limited system sizes (a few hundred atoms at most) and time scales to be simulated. Indeed, the cooling rates in previous DFT simulations of a-Si (≈ 10 14 K/s) are orders of magnitude faster than those in experiments. 20−22 Classical force fields require much less computational effort, giving access to nanometerscale ("device-size") structural models, 9 both for a-Si and for multicomponent systems derived from it (see ref 25 for but one example). However, they are rarely accurate enough to fully describe the structural variations present in the amorphous state.
Capitalizing on today's "big-data" revolution, machinelearning (ML) algorithms are increasingly used to generate interatomic potentials for atomistic simulations. 26−30 By "learning from" (or rather, f itting to) quantum-mechanically computed reference data for energies and forces, ML-based interatomic potentials can enable simulations with an accuracy that is largely comparable to DFT, but with a computational cost that is orders of magnitude lower, and with linear (order-N) scaling behavior. Comparison to experimental observables is thereby the ultimate benchmark and means of validation for the quality of any ML-based interatomic potential, as we stress that no experimental but only DFT-computed data enter the "learning" process.
We believe that such ML potentials are particularly promising for disordered and amorphous materials, which must be represented by nanometer-scale structural models containing several hundreds or thousands of atoms. A landmark example has been the development of an artificial neuralnetwork potential for the phase-change material GeTe, 31 enabling simulation of the crystallization properties 32 including entire nanowires. 33 We recently introduced a ML potential for amorphous carbon, 34 based on the Gaussian approximation potential (GAP) framework 27 and the Smooth Overlap of Atomic Positions (SOAP) atomic similarity kernel, 35 which captures the intricate structural, mechanical, and surface properties of the material 34 and, more recently, has enabled accurate large-scale simulations of the growth mechanism. 36 Very recent work using neural-network potentials allowed for the atomistic modeling of amorphous Li x Si phases relevant in battery applications. 37,38 Finally, such potentials were used in seminal studies to describe the complex phase transitions between polymorphs of crystalline Si. 26,39 In this Letter, we show how realistic atomistic modeling of a-Si can be enabled by a ML-based interatomic potential, again using SOAP and GAP. We first report on melt−quench simulations with cooling rates much slower (that is, better) than what can be achieved in quantum-mechanical-based simulations, and we show how this leads to a higher-quality and lower-energy structure of a-Si. Our structural models show excellent agreement with experiments probing local structure, including 29 Si NMR shifts and diffraction data for high-quality samples, and open the door for future combined modeling and experimental studies on disordered and amorphous materials.
Simulated quenching from the melt is a widely used technique for generating amorphous model networks. In this, one starts with a liquid and progressively lowers the temperature, "freezing in" an amorphous structure. However, for silicon, this approach is not trivial, due to the change in local environments between the high-coordination metallic liquid and the tetrahedral-like amorphous state. We decided to perform a set of variable-volume and constant-pressure (NPT) quench simulations, in which we varied the quench rate, and thus the run-time, by several orders of magnitude. These were carried out using LAMMPS; 40 details are in the Supporting Information. For the moment, we focus on a system size of 512 atoms in the cell and perform a single simulation at each quench rate. This system size is significantly larger than what has so far been accessible to DFT (64−216 atoms), 20−22,24 but smaller than what is possible for empirical potentials; this will Figure 1. Unlocking slow quenching in molecular-dynamics simulations, using the GAP machine-learning framework, and its application to a-Si. (a) Computing time required for MD simulations with different quench rates. Decreasing the quench rate by an order of magnitude increases the number of required MD steps, and thus the CPU cost, by the same factor. Quench rates of ≈10 12 K/s have so far been the limit for 512-atom DFT-MD simulations, and a system size of 4096 atoms ("4k") has been widely out of reach. Both limits can be overcome using GAP. Our "4k" simulation (magenta) employs an adapted temperature protocol; DFT timing information is extrapolated from a short trajectory; see Supporting Information. (b) Stability of 512-atom a-Si structures, taken at various stages of GAP melt−quench trajectories and subsequently relaxed; energies given relative to crystalline (diamond-type) Si. Experimental data refer to samples freshly deposited ("as-dep.") or annealed at progressively higher temperatures. 12 (c) Angle distribution functions for a-Si GAP structures. Points show original data, sampled from short (5 ps) MD simulations; lines show Gaussian fits; data for different quench rates are vertically offset for clarity. (d) Medium-range order in these a-Si networks, assessed by shortest-path ring statistics. 42 The Journal of Physical Chemistry Letters Letter be addressed directly later on. Our fastest quench rate (10 14 K/ s) corresponds to early, seminal DFT studies, 20−22 whereas our simulations at 10 12 K/s mirror the limit of what is presently possible for DFT-quality MD. By contrast, we here use a recently developed GAP model 30,41 which allows us to increase the simulation time 10-fold beyond that, namely, decreasing the quench rate to 10 11 K/s, while retaining similar accuracy.
While an increase in simulation time by 1 order of magnitude may seem incremental at first sight, the full power of ML potentials becomes apparent when looking at the overall computational effort required (Figure 1a). For demonstration, we performed a brief DFT-MD simulation on a 512-atom a-Si network and use the timing information for a rough extrapolation (Supporting Information). Quenching with a rate of 10 11 K/s would thus require around 16 million core hours, or current nominal costs of $185 000 on the UK national supercomputer. In contrast, the same quench rate in GAP-MD required below 40 000 core hours, equivalent to nominal costs below $500. Using GAP, it would hence be possible to decrease the quench rate even further, but given the results obtained at 10 11 K/s, we subsequently chose to increase the system size instead (see below).
The slow quench rate of 10 11 K/s, "unlocked" here using GAP, is indeed required to generate reliable structural models of a-Si. This is seen in Figure 1b: we took structural snapshots at various increments of the simulations, optimized them into local minima, and plotted their energy (relative to the thermodynamically stable form, diamond-type c-Si) as a function of how far the quench has progressed in time from the liquid to the final a-Si structure. The right-hand side shows the experimental sample stability with increasing annealing and thus ordering, based on calorimetry (as is common, we approximate ΔE ≅ ΔH when comparing theory and experiment). Intermediate quench rates lead to a-Si networks that are as stable as freshly deposited or partially annealed samples (ΔE ≈ 0.17−0.20 eV/atom). By contrast, our slowest quench at 10 11 K/s yields a structure whose stability matches the experimental result for a well-annealed sample from ref 12 (ΔE ≈ 0.14 eV/ atom). The GAP-computed bulk moduli for these a-Si networks range from 62−69 GPa and increase with slower quenching (the material becoming "harder"); the computed Young's moduli increase from 73 to 98 GPa; see Supporting Information.
The benefit of slow quenching is further seen in two of the most common structural indicators used for amorphous solids. In a-Si, the bond angles are distributed around the ideal tetrahedral value (109.5°; Figure 1c). Fitting Gaussian distributions to these data allows us to determine the full width at half-maximum (fwhm), which decreases gradually from 30°to 22°with increasingly slower quenching. The experimental value for the bond-angle deviation of ≈11°(ref 10) is consistent with the half width at half-maximum (HWHM) for our slowest quench. Moreover, the mediumrange structural order is important in covalent amorphous networks, 43 and we quantify it here using shortest-path ring statistics. 42 In diamond-type c-Si, all atoms are in six-membered ring (cyclohexane-like; m = 6) configurations, whereas a-Si also contains a large number of fiveand seven-membered rings, and a lesser amount of smaller and larger ones. All rings with m ≠ 6 depart from the reference crystalline state, and as such are a measure of disorder, but we here distinguish them further as follows. Five-and seven-membered rings are still expected to be energetically viable (supported by their abundance in a-Si), whereas, for example, four-membered rings will be clearly under strain. We therefore label rings with m < 5 as "small-ring defects", and rings with m > 7 as "large-ring defects" (Figure  1d).
In Figures 2a−b, we show computed structure factors, S(Q), which can be compared to diffraction experiments. The third peak (at ≈5−7 Å −1 ) gradually splits into two well-defined subpeaks when moving from the 10 14 K/s (yellow) to the 10 11 K/s data (purple). This is qualitatively consistent with experimental observations: as-deposited samples show a fairly featureless third peak, whereas annealed ones (and also our 10 11 K/s result) exhibit a clear splitting into subpeaks. 14 Even better agreement with the experimental structure factor can be achieved for a larger structural model containing 4096 atoms, which we will show below.
We furthermore computed solid-state 29 Si NMR chemical shifts, δ, for all atoms in the unit cells, thereby characterizing each atomic environment individually. We use established DFT-based algorithms 45−47 and reference all δ values to tetramethylsilane (TMS), analogous to experiments. The results for the different GAP structures are shown in Figure  2c (histograms). Furthermore, due to the broad distribution of δ values in the amorphous state, we fit Gaussian profile functions to these data (lines), as detailed in the Supporting Information. We compare the output of these computations to experiments for pure a-Si prepared by sputter deposition. 44 The latter samples were analyzed via secondary-ion mass spectrometry (SIMS), showing no measurable oxygen contamination

The Journal of Physical Chemistry Letters
Letter and ≈0.2 atom % hydrogen in the samples. 44 This low level of impurity is thought to have little or no impact on the 29 Si NMR results, enabling direct comparison to our simulations. In addition to the numerical values reported in ref 44, we fit a Gaussian profile to the experimental data for the sample annealed at 520°C (before the onset of crystallization at higher temperature). We perform this fit using the same procedure as for our DFT data (Table 1). This yields numerical quality criteria that can be used to assess any given structural model.
Clearly, simulations using the two fastest quench rates (yellow and orange) lead to structures with very large scatter in the computed NMR shifts, as a direct consequence of their distorted atomic environments (and thus large fwhm for the Gaussian fits). The network generated at a slower quench rate, 10 12 K/s (red), agrees more appreciably with experiment with regard to both the broadness and the center of mass for the Gaussian fit (δ DFT = −37 ppm); the latter can be compared to δ exp = −38.3 ppm for as-deposited a-Si, and δ exp = −42.9 ppm for a sample annealed at 580°C. 44 Hence, there is a progressive shift to lower frequency in the experimental data with increasing structural ordering, and this is reproduced by our quenched structure at 10 11 K/s (δ DFT = −51 ppm), both qualitatively and quantitatively (to within a few ppm). The results for the 10 11 K/s quench also compare well with those for a DFT-optimized Wooten−Winer−Weaire (WWW) network of a-Si (δ DFT = −53 ppm; structure taken from ref 48), while those for the faster quenches do not (Table 1).
We now place our melt−quench simulations into a wider context, as there are several different ways of modeling a-Si. First, we survey results of RMC modeling, which is an established means of extracting structural information from diffraction data. 17 Recent work by some of us showed that reasonable restraints can improve the RMC modeling of a-Si. 18 In particular, the SOAP similarity measure, initially developed to encode atomic structure in ML potentials, 35 proved useful for this purpose. 48 SOAP-RMC output, subsequently relaxed using DFT, has thus been shown to provide a high-quality structural model of a-Si. 48 We now take the same structures but anneal them further using GAP: heating to 1100 K, holding, and cooling back to 300 K, for a total simulation time of 50 ps. This relatively short annealing is thought to be appropriate, as a recent DFT-MD study showed that annealing a quenched structure at 10 ps versus 20 ps had no appreciable effect on the outcome. 24 We also performed the same annealing procedure Table 1. 29 Si NMR Parameters of a-Si from Simulation and Experiment

The Journal of Physical Chemistry Letters
Letter for the DFT-optimized WWW model from ref 48; a somewhat similar strategy has been followed before, based on a tightbinding model and a system size of 216 atoms. 23 Finally, we include a state-of-the-art 216-atom structure that was carefully generated in a recent work, by slow quenching using the empirical Tersoff potential and subsequent multistep optimization using DFT (here labeled "Tsf+DFT"). 24 We compare these structures in Figure 3 using three types of quality indicators. First, we report the number of coordination defects (Figure 3a), counting 3-and 5-fold bonded atoms with a bond-length cutoff of 2.85 Å. We then measure the distortion from ideal tetrahedral coordination environments: by fitting Gaussians to the angle distributions and determining their fwhm, and by using a numerical order parameter 49 that was employed earlier for tetrahedral environments in liquid water 50 and chalcogenide glasses 51 (Figure 3b). Beyond nearestneighbor environments, we quantify the medium-range order using shortest-path ring statistics, as above, 42 again considering as "defects" any rings with fewer than five or more than seven members (Figure 3c). In all cases, the GAP-quenched structure with the slowest quench rate (10 11 K/s) exhibits very good figures of merit. Interestingly, the count of five-membered rings (m = 5) decreases continuously in progressively more ordered GAP-quenched structures, but that of seven-membered rings (m = 7) increases instead, as shown on the far right of Figure 3c.
Finally, we prepared a larger a-Si structural model containing 4096 atoms (Figure 4a), using GAP-MD and a variable quench rate between 10 11 and 10 13 K/s, as detailed in the Supporting Information. This system size is in reach for ML-based interatomic potentials, 32 as they scale linearly with system size due to their finite cutoff radius (cf. Figure 1a). Having access to ab initio quality structural models on the 4 nm length scale allows us to study the medium-range order more closely. This fundamental question has been discussed in recent work on nearly hyper-uniform networks, 52,53 in particular, by quantifying the inverse height (H −1 ) of the first sharp diffraction peak in the structure factor at around 2 Å −1 . This quantity is taken as a measure for the degree of structural ordering. 52,53 We compare our structure with the current state of the art, viz., a-Si systems containing 100 000 atoms, 9,54 in Figure 4b. Although the latter were generated with an improved WWW algorithm, not by slow quenching, they allow us to place our work in the context of existing ultralarge structural models. Surprisingly, the latter system size alone does not seem to be needed if the structural modeling itself is sufficiently accurate. Indeed, looking at H −1 , our GAP approach outperforms the previous simulation results in much larger cells, and leads, again, to almost quantitative agreement with experiment (H −1 = 0.58 with GAP, H −1 = 0.57 in experiment; Figure 4b). By comparison, an a-Si structure of the same size (4,096 atoms) but generated using empirical potentials gave a much larger H −1 = 0.81 (ref 54). Moreover, our slowest-quenched GAP-based system, even smaller with 512 atoms/cell, yields H −1 = 0.66, remarkably still outperforming the 100 000-atom structure from ref 14 (H −1 = 0.68). Beyond the first sharp diffraction peak alone, Figure 4b also shows that the agreement in the structure factor between the 4096-atom GAP system and experimental data at larger Q is excellent, and significantly better than for the VBSB 100 000-atom system. 9 In conclusion, we have shown that machine-learning-based interatomic potentials can lead to an unprecedented level of quality in the structural modeling of amorphous materials. We used a Gaussian approximation potential (GAP) to generate high-quality atomistic models of amorphous silicon, quenching from the liquid at a rate of 10 11 K/s, hitherto inaccessible to DFT-quality simulations. These structural models agree convincingly with calorimetry, 29 Si NMR experiments, and Xray structure factors, including the height of the first sharp diffraction peak. We note that ML potentials are critically dependent on the quality of the quantum-mechanical input data, and as of today require significant effort to be developed in the first place; in the present case, our GAP has "seen" diverse liquid and amorphous configurations and interpolates between these. These findings will have implications for future research on disordered and amorphous materials, opening the Reciprocal-space fingerprints in the structure factor, comparing to results for two of the largest structural models to date (containing 100 000 atoms, "100k") and to experimental data from ref 14. The height of the first sharp diffraction peak, H, serves as an indicator for structural ordering; in accord with previous literature, we plot H −1 values in the inset. Data for the structure labeled "VBSB" is taken from Vink et al.; 9 data for the structure labeled "HST" is from a more recent study by Hejna, Steinhardt, and Torquato. 52 The Journal of Physical Chemistry Letters Letter door for quantitatively accurate atomistic modeling with direct links to experiments, for a-Si and beyond.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.8b00902.
Computational details; supplementary data and discussion; coordinate files for structural models (PDF) Additional structural data (concatenated XYZ files) (TXT)