Exploring the configurational space of amorphous graphene with machine-learned atomic energies

Two-dimensionally extended amorphous carbon (“amorphous graphene”) is a prototype system for disorder in 2D, showing a rich and complex configurational space that is yet to be fully understood. Here we explore the nature of amorphous graphene with an atomistic machine-learning (ML) model. We create structural models by introducing defects into ordered graphene through Monte-Carlo bond switching, defining acceptance criteria using the machine-learned local, atomic energies associated with a defect, as well as the nearest-neighbor (NN) environments. We find that physically meaningful structural models arise from ML atomic energies in this way, ranging from continuous random networks to paracrystalline structures. Our results show that ML atomic energies can be used to guide Monte-Carlo structural searches in principle, and that their predictions of local stability can be linked to short- and medium-range order in amorphous graphene. We expect that the former point will be relevant more generally to the study of amorphous materials, and that the latter has wider implications for the interpretation of ML potential models.


Introduction
The study of the amorphous state has long been of fundamental research interest. 1 It is also increasingly important to understand structure-property correlations in amorphous materials, owing to ubiquitous applications in solar cells, 2 transparent electronic devices, 3 or phase-change memories. 4 Whilst bulk amorphous phases are challenging to study structurally, twodimensional (2D) amorphous materials can be directly visualized by atomistic imaging techniques such as high-resolution transmission electron microscopy (HRTEM). [5][6][7][8][9][10][11] And just like graphene is the prototypical ordered 2D material, there is ongoing research interest in its disordered analogue(s). Indeed, Toh et al. recently synthesized a centimeter-scale sample of freestanding monolayer amorphous carbon 12 and characterized the structure based on the interpretation of HRTEM images.
The amorphous forms of carbon have been widely studied using computer simulations. Most commonly, these studies have been carried out in the framework of molecular dynamics (MD), from early work on melt-quenching [13][14][15][16][17][18][19] to direct simulations of thin-lm growth by ion deposition. [20][21][22][23] Very recently, the properties of "amorphous graphite", as a 3D extended arrangement of individual amorphous graphene sheets, were investigated with density-functional theory (DFT) based MD and electronic-structure computations. 24 Earlier simulation studies had emphasized the connection between low-density amorphous carbon and the idealized case of 2D amorphous graphene ("aG" in the following). 25,26 Graphene itself, as a 2D system, is rather well-conned and allows for Monte Carlo (MC) simulations to model topological defects. The Wooten-Winer-Weaire (WWW) algorithm, initially proposed for silicon, 27 remains a simple yet robust approach to generating continuous random network (CRN) models. In this, disorder is gradually introduced through local bond transpositions and structural relaxation. To simulate aG, bond transpositions are introduced as Stone-Wales (SW) defects, 28 with moves being accepted or rejected based on a suitable Metropolis criterion, typically the total-energy difference between the new and old conguration. SW defects are created by a formal in-plane 90°rotation of two atoms around the mid-point of the bond 28 and are the foundational example of (intrinsic) topological disorder in 2D carbon. [29][30][31][32] This approach has long been used to study aG. 33 Various implementations of the WWW algorithm exist: D'Ambrosio et al. showed that the majority of bond transpositions are rejected during annealing, and that an early decision scheme can enhance computation speed by rejecting unfavorable transpositions; 34 Ormrod Morley et al. constructed Metropolis criteria from topological metrics such as ring distributions (that is, energy changes were not considered in that case). 35 Machine learning (ML) based interatomic potential models are increasingly being used to accelerate materials simulations. [36][37][38] ML potentials are typically "trained" with DFT data and can achieve similar accuracy for a small fraction of the cost. One key assumption in many of these methods is that the total energy can be separated into sums of machine-learned atomic energies. 39,40 Whilst being an approximation in the rst place, it was argued that these atomic energies may in fact be amenable to interpretation: a connection between ML atomic energies and local chemical structure was made for partly occupied crystallographic sites in b-rhombohedral boron, 41 and for atomic environments with different coordination numbers in amorphous silicon. 42 ML models for other properties of local atomic environments have recently been investigated, ranging from the local electronic density of states 43,44 to local distortion factors in grain boundaries. 45 The nature of these local ML properties (including atomic energies), and their usefulness in predicting physical properties, remains an interesting research question. (See, e.g., ref. 46 for a discussion of ML atomic energies in a chemically complex system.) In the present work, we explore the congurational space of amorphous graphene based on an ML potential model that gives access to total and local energies. We use an MC bondswitching algorithm where ML atomic energies are used in the acceptance criterion and show that doing so leads to physically sound structural models. Depending on the details of the algorithm, we obtain CRN-like or paracrystalline structures. Our work shows that ML atomic energies can be used in different ways to "drive" MC simulations based on local and nearestneighbor (NN) energy contributions, with implications for research on amorphous graphene and likely on other disordered structures and materials.

Potential-energy models
A common ansatz in developing potential-energy models is that the total energy, E, can be separated into a sum of atomic contributions: where 3 i = 3({r ij }), with i, j being atomic indices, and r ij # r cut . This ansatz can be applied to many systems as short-and medium-range interactions predominantly determine the total energy. However, the question of how to formulate 3({r ij }) is not trivial: in carbon, there is a vast congurational space with a subtle interplay between structure and energetics, such as in dihedral and torsional forces, weak interlayer interactions, and so on. Moreover, locality can depend on the structure: numerical experiments showed a large difference in the locality of atomic forces for diamond versus graphite. 47 (We note in passing that extracting local energies directly from DFT is not trivial, although attempts have previously been made. 48 ) Empirical potentials are typically parameterized for a specic composition and phase of a material. For example, the original Tersoff potential for carbon 13 was parameterized by tting parameters of repulsive and attractive pair potentials to cohesive energies of carbon polymorphs along with the lattice parameter and bulk modulus of diamond. This potential 13 used a bondorder approach, where bond strengths are modied according to the number of neighbors. The original reactive empirical bond-order (REBO) potential was an update of the Tersoff potential incorporating hydrogen, 49 with REBO-II adding further improvements. 50 A long-range Lennard-Jones term was added to REBO-II, creating the adaptive intermolecular reactive bond order (AIR-EBO) potential. 51 The long-range bond order potential (LCBOP) 52 is similar to AIREBO, having a long-range term in conjunction with a bond-order description-albeit here it is built in from the beginning. The environment-dependent interaction potential (EDIP) for carbon 53,54 was initially developed from an earlier silicon EDIP model. 55,56 EDIP has been shown to be successful in describing various properties of amorphous carbon, including the graphitization at low and high densities. [57][58][59] Gaussian approximation potential (GAP) The Gaussian approximation potential (GAP) framework 40 is used to "machine-learn" interatomic potential models from quantum-mechanical data, oen based on DFT. Unlike empirical potentials which are constructed based on physical knowledge, GAP makes a non-parametric t. This means that the model can adjust to complex input data-however, it also means that the selection and quantity of reference ("training") data is critically important. 60 In brief, the local energy for a GAP model is where the sum runs over N t training congurations represented by the local-environment descriptor q t , with a corresponding weighting coefficient a t attributed during tting. K is a covariance kernel which measures the similarity between the input and training congurations, represented by q and q t , respectively. A commonly used approach for the latter task is the Smooth Overlap of Atomic Positions (SOAP) descriptor and kernel. 61 For the present work we use the amorphous carbon potential, GAP-17. 47 This model has been shown to predict energies within tens of meV per atom compared to DFT, as well as providing a good description of structural and elastic properties of amorphous carbon. 47

Initialization
We generate a pristine, 200-atom layer of crystalline graphene using the experimental bond length of 1.42Å (ref. 62 and 63). A spacing of 20Å between layers ensures that there is no interaction between periodic replicas. Where noted, NN-energybased runs start from structural models that have been "thermalized", i.e., disordered, using local energies at b = 2.0 eV −1 (see below). The cell parameters are allowed to relax prior to the rst bond switch, and then kept xed for the duration of the MC simulation. The structures are kept planar, simplifying the present study to the idealized 2D case, and noting that puckering can lower the energy further. 64

Monte-Carlo protocol
A kinetic MC algorithm is used to generate aG structural models. The simulations begin either from cG or from a thermalized structure. At each step, a random atom is chosen along with a random neighbor, determined using a cutoff of 1.85Å. The atoms undergo an SW rotation (90°about the bond center) and the new structure is then relaxed using the conjugate-gradient algorithm. For initial testing of structural relaxations, energy evaluations, and MC runs, we used the Atomic Simulation Environment (ASE) 65 interfaced to quippy 66,67 (https://github.com/libAtoms/QUIP). The production MC runs shown in this work used LAMMPS. 68 The force tolerance for structural relaxations was 1 meVÅ −1 with a maximum of 150 relaxation steps. A topological constraint is imposed where the newly relaxed structure must be 3-coordinate, otherwise the move is rejected. (This constraint is included because there will be moves which create stubborn coordination defects, and we found that these defects hindered the progress of MC annealing where the simulation would get "stuck" in local minima.) If the constraints are met, the new structure is accepted with the following probability: we generate a random number z in [0,1), and if w > z, then the move is accepted, else it is rejected and the previous conguration is kept. In this work, b has no physical relation to temperature since DE is not measured for all atoms about their equilibrium position, as the atoms are being held xed aer the bond transposition and thus thermal uctuations are ignored (unlike with MD). Fixing particle positions at every move results in ergodicity being broken and thus samples are not taken from a Boltzmann distribution. 69 Furthermore, the local-energy and NN-energy framework explicitly do not include all atomic energies which is another source of ergodicity being broken. This is not a concern for the present work as temperature-driven dynamics are not relevant (see ref. 69 for an approach to introducing ergodicity to a similar problem) and results in b becoming a tunable parameter. b was initially chosen to be 2.0 eV −1 to correspond to the study by Toh et al., 12 and then tuned heuristically for the local-and total-energy-based MC runs. We use the atomic energies for the defect pair in the Metropolis criterion, dening the local defect energy as and the NN defect energy as where 3 i is the local atomic energy of the i-th atom, with i = 1 and i = 2 denoting the two atoms in the SW defect pair, and i = 3.6 the topological NNs of the SW defect pair ( Fig. 1c and d).
MC runs based on total energies were also carried out for reference. For the 200-atom systems, 25 independent and parallel simulations were conducted over 10 000 MC steps each for different b values and for the 612-atom system, 20 000 MC steps were taken.

Structural analysis
Ring statistics were determined using a shortest-path algorithm 70 as implemented in Matscipy. 71 Topological metrics typically used in network analysis were applied to the 612-atom structures, namely, Lemaître's law 72 and the assortativity 73 (as discussed below). For SOAP structural analysis, we compared each individual atom in a given aG structural model to an atom in cG. The SOAP parameters were: radial cutoff, 5.5Å; cutoff transition width, 0.5Å; neighbor-density smoothness, s atom = 0.5Å; basis-set convergence parameters, n max = l max = 16; dotproduct kernel raised to a power of z = 4. All structures were visualized using OVITO. 74

Phonon calculations
A geometry optimization was performed on the local, NN and Toh et al. structures keeping the cell xed (F max < 1 meVÅ −1 ). We applied a nite-difference approach with a displacement of 0.01Å. The forces in the displaced structures were computed using GAP-17 and second-order force constants were computed using phonopy. 75 A 50 × 50 × 1 k-point mesh was used to determine the vibrational density of states (VDOS). We observed minor imaginary modes in the aG models, presumably due to the constraint of planarity; these are omitted from VDOS plots and will be investigated in future work.

Results and discussion
Local energies for a single defect Fig. 2 characterizes the structure and energetics of a single SW defect in graphene. GAP-17 predicts that the energies of the SW defect pair of atoms are considerably lower in energy than those of the NN atoms. This behavior is in marked contrast with that for the empirical potentials ( Fig. 2b and Table 1): the latter suggest that the SW-pair atoms are relatively higher in energy (with the exception of LCBOP), with atoms in the 5-membered rings being higher in energy overall and those in 7-membered rings being lower. LCBOP and GAP-17 predictions are similar when averaged over the pentagon-pair environments, only differing slightly in structural parameters. REBO and AIREBO show slightly larger SW-pair energies and more strained structures ( Table 1).
The total defect formation energy DE total for GAP-17 is similar overall to that predicted by most empirical potentialsand to DFT values, which we obtained using CASTEP. 79 In an earlier study, hybrid-DFT and quantum Monte Carlo calculations predicted values of 5.69 eV and 5.92 AE 0.03 eV, respectively. 30 EDIP, LCBOP, and REBO were parameterized for sp 2 environments, and their DE total is in close agreement with that for GAP-17. AIREBO was also parameterized for sp 2 environments, however the added long-range term may not be optimized for 2D planar carbon. In contrast, D3 local and D3 NN are not in close agreement, where D3 local differs the most across all the potentials (except for REBO and AIREBO which are parameterized identically for short-range terms). Overall, GAP-17 predicts the lowest energy for the NN defect formation energy among the potential models investigated.
Structural parameters as predicted by the respective potentials are also shown in Table 1. The defect-pair bond length (d 1 ), defect-NN bond length (d 2 ), and defect bond angle, q (Fig. 1b) are in close agreement across all potentials-with the exception of EDIP, which notably increases the cell parameters upon relaxation. This is shown by the difference in cG bond length (d 0 ) between EDIP and the other potentials when the cell parameters are allowed to relax, with EDIP predicting an elongated d 0 . GAP-17 and LCBOP are in close agreement with regard to q, predicting a slightly wider bond angle for the 5j7 pair, whereas EDIP predicts it slightly lower, notwithstanding the fact that the EDIP defect formation energy is in close agreement with that for GAP-17. Furthermore, GAP-17 is in excellent agreement with the bond lengths and angles from our DFT calculations.  Structure and formation energy of a single SW defect in graphene from ML and empirical potentials as well as DFT. The table shows the computed bond length in pristine graphene (d 0 ), the defectpair (d 1 ) and defect-NN bond lengths (d 2 ), the defect bond angle (q; cf. Fig. 1b), and three energies associated with defect formation: the local energy change in the rotated atoms only (D3 local ), the NN energy change (D3 NN ), and the total defect formation energy (DE total ) Using local energies to drive Monte-Carlo annealing Fig. 3 shows the energy proles of independent parallel MC simulations via successive SW transformations, which we use to create structural models of aG. Panel (a) shows the evolution of the average energy of the ensemble when using (only) the atomic energies of the SW defect pair in the Metropolis criterion. Clearly, using a b value of 2.0 eV −1 leads to a highly disordered 3-coordinate structure (as shown in the inset), with energy convergence occurring aer about 1000 MC steps. Such structures contain high-energy 3-and 4-membered rings with signicant strain. We call the resulting structures "thermalized" and use them as starting points for the runs characterized in blue in Fig. 3b, as indicated by an arrow. Nearest-neighbor (NN) atoms are dened topologically in the present work, and so we include the sum over the SW pair and its bonded neighbors, including six atoms in total (eqn (5)). The energy difference before and aer the SW transformation is used in the Metropolis criterion. Results from two protocols, starting either from cG or from the previously thermalized structures, converge to similar energy values (green and blue curves in Fig. 3b). It seems that runs from both starting points tend toward a mutual limiting distribution. This is reected in Fig. 3b, as both curves have means and standard deviations in very close agreement aer approximately 4000 MC steps. However, it is apparent that this protocol is not ergodic, as discussed in the Methods section. The principle of ergodicity states that all possible congurations of the system should be attainable 80 and thus it is clear that given the topological constraints imposed, it is impossible for certain proposed moves to be accepted even though they may be energetically favorable. Fig. 3c illustrates the effect of varying b on the resulting energy prole when using the local-energy criterion. It is clear It is seen that a smaller b (higher simulation "temperature") leads to a higher-energy structure and faster convergence. We emphasize that despite the high standard deviations, the lowest occurring energy value is DE = 0 eV at. −1 corresponding to ideal crystalline graphene. (d) Heuristically tuning b to find a mutual limiting distribution between MC approaches (within statistical fluctuations). that using a lower b value corresponds to a higher-energy nal structure as well as rapid convergence: for b = 2.0 eV −1 , convergence occurs in fewer than 1000 MC steps (black line). Increasing b results in lower-energy structures and slower convergence, as shown by the b = 50.0 eV −1 and b = 60.0 eV −1 results. For b = 60.0 eV −1 , convergence has not occurred aer 10 000 steps, however the mean appears to be tending towards that for b = 50.0 eV −1 , implying that there might be a minimum level of disorder that is "stable" for the local-energy criterion. Additionally, as b was increased beyond 60.0 eV −1 , proposed moves were always rejected, further suggesting a minimum convergence energy using this framework. Fig. 3d shows the results from three different runs with heuristically tuned b values for local-, NN-, and total-energy criteria, respectively, converging towards a mutual distribution (within statistical uctuations). This result allows the direct comparison of structures generated by the different frameworks side-by-side, as the effect of b has been removed. Such a comparison is particularly instructive for lower-energy structures, as we will show below: they show a richer congurational space, including medium-range order (as compared to the highly disordered high-energy structures), and this space may be traversed differently by different MC runs.

The congurational space of amorphous graphene
Having explored different protocols for generating aG structural models, we next created larger-scale structures: the system size was increased to 612 atoms and structures from local-and NNenergy-based searches were studied. Additionally, the 610-atom structural model published by Toh et al. (ref. 12) was included for comparison. Fig. 4 shows these three congurations. As we keep the cells xed throughout the Monte-Carlo searches (Methods), we computed the tensile stress for the nal structures, using GAP-17, and found them to correlate with the degree of disorder. The structure generated using local energies had an in-plane tensile stress of 3.53 GPa; the NN structure had 2.85 GPa, and for the Toh et al. structure we obtained 0.87 GPa. For the local-energy framework, the absence of crystal-like pockets is clear as seen in panels (a) and (d). The structure resembles a CRN, with chains of 5-and 7-membered rings running across the cell, and there are more large rings compared with the other structures. The pronounced disorder is likely a result of the low ML energy for the SW defect itself (Fig. 2a), and thus of the low energy cost for these transpositions (given that medium-range order is not captured in the local energies alone).
With NN energies used in the Metropolis criterion, panels (b) and (e) suggest that using an NN criterion encourages small pockets of crystal-like regions forming, indicating that direct contributions from NNs maintain medium-range order, whilst retaining an amorphous structure at b = 2.0 eV −1 . Visually, the locally averaged energies (up to NNs) in panels (a-c) suggest that using just the local defect-pair energy gives more regions of higher energy (shown in yellow in Fig. 4a) compared with the NN framework (Fig. 4b). This appears consistent with the nature of CRNs versus paracrystalline structures, and color-coding by NN energies shows a clear difference between the two structures.
In the NN-energy-based structure (Fig. 4b), there are pockets of crystallinity, indicated by regions of ordered 6-membered rings. Interestingly, there is an aggregation of more disordered regions, with 5-membered rings surrounding larger 7-and 8membered rings. For the structure shown in Fig. 4c, the authors started from a randomized, hard-sphere constrained structure and worked towards a paracrystalline sample using the AIREBO potential. In the resulting structure, we note the presence of coordination defects, since the randomized initial structure was not topologically constrained. There is also a 4-membered ring (pink in Fig. 4f). It is evident that this sample is paracrystalline with regions of locally crystal-like order separated by 5j7 grain boundaries, and that it is more ordered than the NN-based structure (Fig. 4e). Larger defects appear to congregate as seen, for example, in the top of the gure. These defective environments are high in energy, as shown by the color-coding. Fig. 4g-i shows the distribution of NN energies for the respective structures. The average energy for the atoms in the local-energy-based structure is 0.33 eV at. −1 with a standard deviation of s = 0.12 eV at. −1 , and the energy of the NN-based structure is 0.28 AE 0.12 eV at. −1 above pristine graphene. Both values are close to the mean energy resulting from the independent runs for the 200-atom system (Fig. 3b and d): hence, an individual 200-atom run will not likely be sufficient to describe aG, but an ensemble of multiple independent runs will be-just like many small-scale structural models of 3D amorphous carbon have been used in the tting of GAP-17. 47 The distribution in the locally averaged energies is similar overall between Fig. 4g and h. For Fig. 4i, we nd a bimodal distribution in the energy, indicating that this structure shows higher paracrystalline order compared to that generated by the NN framework. There exist pockets of more ordered, more energetically favorable regions in between more disordered ones, as seen in the structural model shown in panels (c) and (f). Very recently, a paracrystalline sample of diamond has been synthesized. 82 This discovery, along with the synthesis and characterization of paracrystalline graphene 12,83 suggests that the landscape of disorder in carbon may have a link between "fully" amorphous and crystalline. This is further reected in the NN energy distributions in Fig. 4, where panel (g) clearly shows a single distribution, panel (i) a bimodal distribution, and panel (h) a small hint of one. These ML energy distributions may provide a quantitative distinction between CRN and paracrystalline graphene, where a prominent bimodal distribution indicates the latter and the lack thereof indicates the former.

Quantifying disorder
We show further, quantitative structural indicators in Fig. 5. The bond-angle distribution (panel a) characterizes mediumrange order. The local-energy-based structures (gray) display a wide range of bond angles, reecting the relatively high level of disorder. The absence of a clear peak at 120°also emphasizes the reduction of medium-range order as is characteristic of CRNs. The NN curves (green) are narrower and centered around 120°, consistent with locally crystal-like environments. A shoulder peak at z109°shows the strong presence of 5membered rings. The Toh et al. structure has a pronounced peak at 120°with a smaller shoulder at z109°, indicative of the larger degree of locally crystal-like order. The bond angles in the 200-atom models (solid lines) agree well with those in the 612atom ones (dashed).
The count of shortest-path rings is another metric for medium-range order. Ring counts for the local-energy MC runs reect the large disorder, showing more 5-and 7-membered rings than in the other structures. The structures from the NN runs have ring counts centered around 6, as for the Toh et al. structure, with a larger count suggesting increased crystallinity in the latter case. As with the bond-angle distributions (Fig. 5a), the ring statistics for the 200-atom versus 612-atom structures are in close agreement within each other, within the standard deviation of the values for the former (Fig. 5a).
SOAP is a structural similarity metric for atomic environments. 61 2D plots can reveal correlations between SOAP similarity on the one hand, and locally averaged energies on the other hand. 42 Fig. 5d shows the distribution of locally averaged energies (up to NN) and SOAP similarity to cG. The wider the spread in both axes, the greater the structural disorder. Comparing local-energy-with NN-energy-based data, both KDE curves are in close agreement, being skewed slightly toward the le (SOAP) and to higher energies for the runs using local energies only. This is expected, since the ML model predicts relatively low energies for the SW pair (Fig. 2), allowing structures to become more disordered. When averaging over NN environments, the distribution shis and narrows as we observe fewer highly disordered environments. Data for the structure from Toh et al. 12 are provided to show how ordered this structure is. Defective environments are clearly identied with a series of data points at higher energies. As with Fig. 4, there is a clear bimodal distribution in the KDE curve. As seen for the structural indicators characterized in Fig. 5a and b, the dashed lines representing the 612-atom structures agree well with the set of separate 200-atom structures.
In addition to the well-established ring statistics and the SOAP similarity, we analyzed the structural models with two topological metrics typically used in network theory, viz. Lemaître's law 72 and the assortativity metric. 73 Lemaître's law connects the second moment of the ring distribution (m 2 = <k 2 > − <k> 2 , where k is the ring size) with the fraction of sixmembered rings, p 6 . Since we impose the constraint that all atoms must be 3-fold connected, it follows that the distribution can be explained well by a single maximum-entropy distribution leading to a characteristic curve. 72 As lim p 6 /0.6 from p 6 = 1, m 2 increases linearly as 1 − p 6 in the region of . This line is extended (red line) to show that most of the points from both local (b = 50.0 eV −1 ) and NN frameworks are located on this curve. This curve corresponds to the maximum-entropy solution if only 5-, 6-, and 7-membered rings are present, which is the case in both these frameworks. m 2 then increases exponentially beyond this where structures from the local framework at b = 2.0 eV −1 are located. These structures show greater disorder and a wider spread of m 2 values, as expected due to the presence of 3-and 4-membered rings, where the maximumentropy solution follows an exponential prole. For the 612atom structures, the local framework (b = 50.0 eV −1 ) based structure is more disordered than the NN-based one and the Toh et al. structure (orange). The latter does not yield a datapoint on the curve as it has a (small) number of 2-coordinate sites and hence would be located on a different Lemaître curve. 35 The assortativity, r, measures how likely a large ring is to be next to a smaller ring (disassortative, r < 0) or to other large rings (assortative, r > 0). In many physical systems, one nds a preference for disassortative congurations. This is, here, reected in the range of r values (Fig. 5e). For the 200-atom systems, we nd no discernible correlation between r and p 6 . For the 612atom ones, both the local and NN frameworks yield structures with similar r values. The Toh et al. structure has a slightly less negative r, implying a slight preference for a more random arrangement. This may be due to coordination defects in the structure skewing r toward zero. We nally considered properties that can be compared to previous or future experiments. Vibrational properties of graphene are readily characterized by Raman spectroscopy 84,85 and computations, 26,86 and we show predicted VDOS in Fig. 6a. For pristine graphene (as reference), we obtain a characteristic peak at 1510 cm −1 and a smaller one at 1410 cm −1 , in close agreement with previous computations 26, 86 and experimental spectra. 84,85 As disorder increases, the 1510 cm −1 peak shis to higher frequencies, as shown in Fig. 6a: to about 1520, 1570, and nally 1610 cm −1 .
The peaks at 1510-1610 cm −1 seen in Fig. 6a for the disordered structures are close to Raman data by Toh et al., with the corresponding peak at 1558 cm −1 . 12 The blue shi of the 1510 cm −1 peak with increasing disorder is a well-known effect and has been observed in VDOS and Raman spectra. 26,86 Furthermore, we observe an "amorphous" peak in the acoustic (low-frequency) region, at 814-836 cm −1 , red-shiing with increasing disorder.
Another way of validating against experiment is to evaluate the radial distribution function (RDF) and X-ray diffraction (XRD) patterns: in disordered structures, characteristic peaks will become wider and less intense. 12,87 However, when comparing amorphous structures, differences in RDF and XRD are subtle: 18,87 in Fig. 6b, the second RDF peak at 2.41Å broadens with increasing disorder, but this effect is subdued for the remaining peaks at 2.82 and 3.62Å for the local and NN structures. The Toh et al. structure has longer-range order, re-ected in more prominent RDF peaks at greater distances (4.23 and 4.93Å). The local and NN structures have peak locations comparable to that of the pair correlation function recreated from TEM in ref. 12.

Conclusions
Atomic energies predicted by the GAP ML framework can be used to drive Monte-Carlo structural exploration in principle. We have shown this by creating structural models of amorphous graphene, one of the prototypical disordered systems in physics and chemistry. We found that using (only) ML atomic energies leads to structures resembling continuous random networks, whereas including nearest-neighbor energies tends to drive the simulations toward some degree of paracrystalline order. We suggest that histograms of local energies, as shown in Fig. 4, can give insight into the degree of "randomness" in different amorphous networks: deviations from random order are visible as peaks at either low (paracrystalline) or high (coordination defects) energy.
The demonstrated ability to use ML local energies in MC annealing indicates potential for future research on amorphous materials. Describing local environments using ML methods can provide insight into the relation between atomic structure and energetics, and therefore structural stability, in amorphous materials-including amorphous carbon, which has emerging applications in biosensing 88 or batteries. 89 The fact that NNaveraged energies yield reasonable, partly paracrystalline structural models may be attributed to the fact that they provide "smoothing" over the variance in local atomic energies. This nding is consistent with earlier ndings for the electronic DOS 43,44 and might have wider consequences for ML predictions of local properties, which are yet to be fully explored.

Author contributions
Z. E.-M. developed the Monte-Carlo protocols and carried out the computational work. V. L. D. initiated and coordinated the study. All authors made substantial contributions to data analysis and discussions. Z. E.-M. and V. L. D. wrote the paper with input from M. W., and all authors approved the nal version.

Conflicts of interest
There are no conicts to declare.