Towards an atomistic understanding of disordered carbon electrode materials

Machine-learning and DFT modelling, linked to experimental knowledge, yield new insight into the structures and reactivity of carbonaceous energy materials.


Supplementary discussion (I): PDF data
Figure 2c of our Communication shows PDF data for a GAP-generated porous carbon structure and for two representative experimental samples. It is stated in the main text that our simulation "reproduces all general features". Here, we discuss in more detail the remaining deviations between computed and experimental results. Fig. 2c in the main text). From top to bottom: Pair distribution functions for a porous carbon structure (scaled down arbitrarily by dividing by a scaling factor of 2), for a carbide-derived carbon from ref S1(shown in Fig. 1d and labelled "TiC-CDC-600" there), and for a hard carbon from ref S2 (shown in Fig. 1a and labelled "pristine" there). Vertical offsets of +6 (-3) Å -2 for the GAP (hard carbon) data are chosen to ease visualization, respectively.

Fig. S1 (as supplement to
There are four main differences between simulated and experimental PDFs. These differences are highlighted by shading in Fig. S1: • The first peak (A) is less broad in the GAP data than in both experiments. Our calculated value refers to a single structure (from an MD simulation at 300 K), but has otherwise not been thermally broadened. There are two possible routes by which one could do so. We calculated an averaged PDF for ten different structures from a short (1,000 steps ≡ 1 ps) MD trajectory, but this led to only insignificant changes. We also investigated the effect of artificially adding broadening, by convoluting the simulated PDF data with Gaussians of different widths (0.11 Å, as suggested in Ref. S3, as well as smaller values). The latter procedure, however, led to loss of features: for example, S3 the previously well-defined peak at 2.8 Å was smoothed out to appear like a shoulder of the peak at 2.5 Å only. We therefore refrain from applying Gaussian broadening to our calculated data.
Some previously presented computational results, obtained with classical interatomic potentials, do show broader distributions in the first PDF peak. We note, however, that these simulations do not follow a quenching trajectory down to 300 K but are quenched only to higher temperatures. This precludes a direct comparison to our data.
The quench approach we are using is not a mimetic one, as it here serves only to generate an initial amorphous precursor (in the spirit of commonly used modelling approaches such as in Ref. S3). Other structure-generation strategies might be interesting to explore in the future, and could likewise be coupled to the GAP framework.
• The overall magnitude of the first peak (A) is visibly larger in the GAP data than in experiments (recall that scaling has been applied for plotting): numerical integration yields ≈ 1.98 for the GAP structure, compared to ≈ 1.15 for the experimental "porous" carbon (TiC-CDC). To rule out computational artifacts, we integrated over the first peak to obtain the coordination number (CN), which gave 3.25, close to the value of 3.0 in ideal graphene and graphite. For the TiC-CDC, the same procedure yields a CN close to 2 when assuming a density of 1 g cm -3 . Note, however, that densities for these experimental samples are difficult to determine, and might be underestimated here (leading to an underestimation of the CN). To account for higher sample densities, we also analysed the PDF for a GAP structure at ≈ 1.5 g cm -3 , which led to similar conclusions as the ≈ 1.0 g cm -3 one.
• A shoulder on the low-r side of the second peak is seen in the GAP data (B) but not in experiments. This stems from the presence of 5-membered rings, as clearly identified by our analysis in Fig. 2d. We re-iterate that the count of 5-/7-membered rings in our structures is high, exceeding that in more graphite-like structures, but in turn allowing us to sample diverse local environments. Future work will deal with likewise more ordered GAP-derived structural models.

S4
• The magnitude of the peak at 3.8 Å is somewhat higher in the GAP data (C) than in experiments. We believe that a similar argument might hold as made for the higher first peak (A). This signal has been linked to one of the cross-ring distances in two adjacent 6-membered rings, S1 so an overestimation as in case B is unlikely to occur here.
• Additional contributions are seen in the computed data around 4.7 Å (D), whereas no such signals appear for the experimental samples. Although it is difficult to make a final statement due to the disordered nature of our structures and the (after all) limited sample size, we did inspect interatomic distances (Fig. S2), and found several occurrences corresponding to the distance region highlighted as (D) in which odd-membered rings are involved. We therefore believe that this region of the PDF can serve as another "fingerprint" for disorder in carbon networks, and that these contributions there will diminish upon further annealing. This will be investigated in future work.  S2. Structural fragment from a GAP-generated carbon structure, emphasizing interatomic distances that correspond to the additional PDF signals highlighted as (D) in Fig. S1. Pairs of atoms are colour-coded individually, and it is seen that several of them involve structural defects or 5-/7-membered rings (one pair of such odd-membered rings is highlighted by an arrow).

Supplementary discussion (II): Pore sizes
To gain further insight into the nature of the pores in our low-density structural models, we used the POREBLAZER 3.0.2 software S4 (with default parameters) to calculate geometric pore size distributions. In brief, this is realised by using a two-steps Monte Carlo procedure. In the first step, a random point a is chosen in a fine 3D grid of points which divides the simulation cell into bins. The program checks for overlaps with the neighbouring atoms considering a van der Waals diameter of 3.431 A for the carbon. If point a is not overlapping with the solid, a second Monte Carlo step is done to find the largest sphere that encompasses point a and does not overlap with the carbon structure. This process is repeated many times and the cumulative pore volume function Vp(r), i.e. the void space that can be covered by spheres of diameter r or smaller, is built. The geometric pore size distribution is then obtained by differentiating Vp(r). Details of the approach and references to further literature are in Ref. S4.
Results for the small structures (≈ 200 atoms) and a density of ≈ 1 g cm -3 are shown in Fig. S3 as the most representative case. Most of these structures have a single pore occupying a large part of the cell (the precise volume being dependent both on the carbon network and on the optimised mass density); a few structures show a bimodal distribution. The observed pore sizes range from 0.8 to 1.3 nm which is close to what is observed experimentally for TiC-CDC samples of similar densities. S1,5,6 The generated structures thus provide a range of models compatible with the experimental results, both in terms of local structure ( Fig. 2 and S1) and porosity. Nevertheless, it is important to point out that a single small structure (at the size of what is accessible to DFT simulations), having a limited number of pores or even one, is insufficient to represent correctly the wider distribution of pores observed experimentally for bulk materials.  S4. Electronic structure of disordered carbon systems, generated in GAP-MD simulations, at densities of 1.0 g cm -3 (a), 1.5 g cm -3 (b), 2.0 g cm -3 (c) and of idealised graphene for comparison (d). The panels on the left-hand side give total electronic densities of states (with the Fermi level, EF, set as the energy zero), as well as projections onto the s (blue) and p (teal) valence orbitals. The panels on the right-hand side show a close-up of the energy window around EF, and additionally the inverse participation ratio (IPR) for each Kohn-Sham eigenvalue is given (green): larger values indicate localised states, whereas low values indicate higher delocalisation. The IPR around EF is notable for the disordered systems (panels a-c) but vanishingly small for the idealised graphene system (panel d).
Having access to optimised carbonaceous structures allows us to feed these into DFT computations. In the main text, we focus on an exemplary study of Na intercalation (Fig. 3-4), but the pure frameworks and their electronic structure are amenable to (single-point) DFT analyses as well. Other than in the main text, we are here interested in the band gap of the pure structures and therefore employ a more computationally expensive, higher-level DFT method for this particular purpose.
As another example, and to illustrate the usefulness of the approach, we here investigate the electronic structures of systems containing ≈ 1,000 atoms at densities around 1.0 g cm -3 , 1.5 g cm -3 and 2.0 g cm -3 , respectively; for comparison, we study a similarly sized idealised graphene system as well (Fig. S4). These computations were performed using CP2K and the hybrid HSE06 functional, S7,8 which is known to provide accurate electronic band gaps in solids and has been applied to smaller amorphous carbon structural models in a recent study. S9 Overall, we find that, as the density increases, and the atomic structure increases in likeness to graphene (cf. Fig. 1 in the main text), the electronic structure tends toward graphene too. All porous carbon structures exhibit a computed band gap of zero eV, and the states near the Fermi level, EF, correspond to p-orbital contributions (Fig. S4a-c), in line with the simplified notion of π systems that do not hybridise with the valence 2s orbitals. The total DOS of graphene (Fig.   S4d) has a sharp drop-off near EF, and the remaining finite DOS is merely a consequence of the Γ-point sampling and Gaussian broadening applied to the electronic levels. As the density of the amorphous structures increases from 1.0 to 2.0 g cm -3 , the absolute DOS at EF decreases by approximately two thirds, again emphasizing the closer correspondence to idealised graphene.
The states near EF are further analysed using the Inverse Participation Ratio (IPR), which quantifies the localization of the Kohn-Sham eigenstates and is inversely proportional to the number of atoms on which a given eigenstate is localised. The number of states with an IPR value > 0.2 in the 5 eV range below EF decreases from 15 to 10 to 2 as the density increases from 1.0 to 1.5 to 2.0 g cm -3 . Visual inspection of the wave functions revealed that localised states near EF are mainly localised on regions of disorder, such as odd-membered rings or undercoordinated sites.
Amorphous porous carbons are used as electrodes in supercapacitors, as ions can enter the pores. When applying a voltage to the electrodes, a high electronic conductivity is desirable; however, theoretical and experimental work has shown that the conductivity of amorphous carbons is much lower than that of graphene. S10-12 The localised states near EF seen in Fig. S4 indicate that electron hopping between localised states is required. As the localisation near EF decreases with increasing carbon density, the electronic conductivity is likely to increase.

Supplementary discussion (IV): Na intercalation from DFT-MD
As stated in the main text, we performed DFT-MD simulations (annealing at 1,000 K for 10 ps and subsequent cooling to 10 K) for the same carbon structure but with different numbers of intercalated Na atoms. Figure S5 provides an overview of the resulting structures, with an emphasis on the behaviour of Na atoms. In the case of low filling (a), the atoms are isolated, consistent with the notion of positively charged "Na + " ions repelling each other. On increasing the concentration of Na atoms slightly (b), we observe the formation of a five-membered cluster at the endpoint of the MD trajectory, whereby we define a "cluster" through connectivity, drawing bonds up to a maximum Na-Na distance of 3.8 Å. A further increase of the Na concentration leads to the formation of a larger, nine-atom cluster in the large pore of the same host structure (and only there). Figure 4 in the main text shows how the occurrence of such Na-Na contacts and, ultimately, clustering, is concomitant with lowered charges on the atoms involved: that is, with a gradual back-transfer of electrons and a transition from "Na + " to "Na 0 " (see discussion in the main text).

Fig. S5.
Structural snapshots after DFT-MD annealing and cooling, in parallel runs, for a carbon framework with different numbers, NNa, of intercalated atoms. The upper panels show the entire structure (carbon atoms as wireframe, Na atoms as yellow spheres; viewed down the a-axis of the simulation cell). The lower panels show close-ups of the resulting clusters, with Na-Na distances given in Å. We draw bonds between Na atoms with a spacing of 3.8 Å or closer.
It is stated in the main text that the Na filling and high-temperature MD simulations leave the carbon structure "largely unaffected" (p. 3). To support this statement, we provide here structural drawings of the same carbon framework in different scenarios, omitting Na atoms for clarity (Fig. S6). From left to right, we show i.
an exemplary, DFT-relaxed structure resulting from the ab initio random structure searching (AIRSS; Ref. 43 in the main text) like procedure used in Fig. 3 (insertion of a single Na); ii. a structural snapshot after 5,000 timesteps (5 ps) of high-temperature (1,000 K) DFT-MD annealing of an Na14C206 structure (i.e., a case with much higher filling); iii. the same after 10,000 timesteps (10 ps), indicating the stability over time of the carbon framework; iv. the same structure after cooling in MD, which allows the Na atoms to settle into their local minima.

Computational details
Partial ring analysis. The contributions from different ring types (as shown in Fig. 2d in the main text) were analysed for internal distances only (i.e., distances inside the rings). The R.I.N.G.S software S13 was used to identify all 5-, 6-and 7-membered rings in the GAP-generated structures and to provide lists of the carbon atoms belonging to these different rings. This output was then used to calculate carbon-carbon distances for each of the lists, thus separating the PDF contributions from 5-, 6-and 7-membered rings. We could check easily that the sum of all these contributions is equal to the total PDF calculated in the usual way.
GAP modelling. The QUIP / quippy code, which is freely available for non-commercial research at http://www.libatoms.org, was used for all GAP-driven MD simulations. Amorphous precursors were generated by rapid quenching from the melt as described in our previous publication S14 and in line with established procedures in the community. S15 The following protocol was then used to generate porous and partly graphitised structures: • Heat to 3,000 K over 10 ps, with a stepwise increase in temperature every 1 ps; • Anneal at 3,000 K for 100 ps; • "Prune" the resulting structure of unfavourable fragments: o identify all atoms with a CN of < 2 ("dangling-bond defects") and o identify all atoms with CN = 2 whose both neighbours similarly have CN = 2 (atoms in close-to-linear "sp-like" chains); remove all atoms so identified; o repeat the above steps until no unfavourable fragments remain in the structure.
• Anneal at 3,000 K for a further 100 ps; prune as above; anneal for a final 100 ps; • Cool to 300 K over 10 ps, with a stepwise lowering of the temperature every 1 ps; • Optimise the resulting structure, using the following protocol: o create five copies of the cell with the lattice parameters scaled by -2%, -1%, 0%, +1%, +2%, respectively; for each of these copies o perform an MD simulation with exponential quench from 300 K to 10 K over 5 ps; o then optimise atomic positions using the LBFGS algorithm; S16 o fit the resulting E(V) data to the Birch-Murnaghan equation of state; S17 o perform a final MD quench and optimization as above at the so-determined optimum unit-cell volume.
The timestep was 1 fs in all simulations.

S11
The above protocol was carried out for initial densities of 1.0, 1.25, 1.5, 1.75, and 2.0 g cm -3 , respectively; the final densities varied slightly from sample to sample, due to the removal of atoms in some cases and the optimization of the unit cell; results are given below. We studied 10 systems in parallel at an initial system size of 216 atoms/cell, giving rise to 50 candidate structures, which were subjected to DFT computations as below. We furthermore created 3 systems with 1,000 atoms in the cell and initial density settings of 1.0, 1.5, and 2.0 g cm -3 , using the same GAP annealing and optimization procedure but not subjecting these to further DFT optimization, due to the large cost of performing structural relaxations at this system size.
Among the 5 × 10 = 50 "small" structural models (≈ 200 atoms/cell), four showed the formation of unfavourable chain motifs after optimization, and these were therefore excluded. One of the DFT relaxation. DFT computations were carried out using the projector augmented-wave (PAW) method S20 as implemented in VASP. S21-24 The plane-wave energy cut-off was 500 eV, the SCF halting criterion was ∆E < 10 -6 eV, and reciprocal space was sampled at Γ. The PBE functional S25 together with pairwise Tkatchenko-Scheffler dispersion corrections S26,27 was used for structural optimisations, such as to minimise residual forces on atoms below 0.01 eV Å -1 .

Single-point energies.
It is known that many-body dispersion (MBD) effects, beyond pairwise additive terms, play an important role in carbon nanostructures. S28 We therefore computed single-point energies for the optimised pure carbon structures using the MBD approach by Tkatchenko and co-workers, S29-31 as implemented in VASP by Bučko and co-workers. S27,32 In contrast, intercalation energies for single Na atoms (translated into voltages), were computed at the PBE level, as these pertain to relative energy differences between two largely similar carbon networks and to the Na-C interaction instead.

Structural details
Due to the large file sizes, all structures are provided as a separate ZIP file (available as ESI), containing POSCAR (VASP input) as well as XYZ format.
Relevant details of the generated 216-atom structures are given in Table S1 below, namely: a) the initial setting for sample density (at which the amorphous precursor was generated); b) a running index (as used in the attached ZIP file to identify individual structures); c) the "optimised" density, that is, after removal of defects during annealing and subsequent structural optimization (including optimization of the cell volume); d) the total VASP-computed energy at the PBE+MBD level (see Computational details section above); e) the relative energy, given as the difference from the most stable structure found, and per atom; f) a qualitative identification of the structure type (based on visual inspection), labelling each entry as "porous" (P) structures containing voids, mainly found at lower density, or "graphite-like" (G) showing extended defective 2D sheets, mainly found at higher density (see Fig. 1 in the main text for an exemplary illustration). Table S1. Overview of ≈200-atom structures generated in this work (see text for details).