An Atomistic Carbide-Derived Carbon Model Generated Using ReaxFF-Based Quenched Molecular Dynamics

We report a novel atomistic model of carbide-derived carbons (CDCs), which are nanoporous carbons with high specific surface areas, synthesis-dependent degrees of graphitization, and well-ordered, tunable porosities. These properties make CDCs viable substrates in several energy-relevant applications, such as gas storage media, electrochemical capacitors, and catalytic supports. These materials are heterogenous, non-ideal structures and include several important parameters that govern their performance. Therefore, a realistic model of the CDC structure is needed in order to study these systems and their nanoscale and macroscale properties with molecular simulation. We report the use of the ReaxFF reactive force field in a quenched molecular dynamics routine to generate atomistic CDC models. The pair distribution function, pore size distribution, and adsorptive properties of this model are reported and corroborated with experimental data. Simulations demonstrate that compressing the system after quenching changes the pore size distribution to better match the experimental target. Ring size distributions of this model demonstrate the prevalence of non-hexagonal carbon rings in CDCs. These effects may contrast the properties of CDCs against those of activated carbons with similar pore size distributions and explain higher energy densities of CDC-based supercapacitors.


Introduction
Carbide-derived carbons (CDCs) are a class of porous carbons with well-ordered porosities and heterogeneous, short-range graphitic ordering [1]. The most common synthesis approach uses systems on the order of a few hundred or thousand atoms are sufficient for capturing short-range features such as bonding, ring formations, and a small number of structural features and pores. Larger system sizes on the order of tens or hundreds of thousands of atoms are needed to better sample a range of structural heterogeneities and pore size distributions.
As an alternative approach, quenched molecular dynamics (QMD) has viably reproduced atomistic models of amorphous and/or porous carbons [53][54][55], silica [56][57][58], and metals [59,60]. QMD generates metastable structures by lowering the temperature of a system over the course of a simulation. Although this process lacks a direct physical analog (i.e., in practice, CDCs are typically formed via chemical etching), QMD is, nonetheless, capable of generating realistic atomistic models of porous materials. This technique allows quench rate and initial and final temperatures to be modulated in order to control some of the structural properties of the CDC, such as porosities and RDFs. The inherent parallel nature of molecular dynamics simulations allows them to explore large system sizes. Previous QMD results have used force fields such as the Adaptive Intermolecular Reactive Empirical Bond Order (AIREBO) [61], the Environment-Dependent Interatomic Potential (EDIP) [62], the Reactive Summation State (RSS) [53], and the Tersoff potential [63] to generate nanoporous carbon structures for use in molecular simulation. While many impactful studies [49,52,54,55,[64][65][66] have been done with these forcefields, the generated structures have often failed to accurately reproduce experimental pore size distributions and always depend, at least in part, on the accuracy of the force field itself. Carbon-based materials, and CDCs in particular, present a unique challenge for these forcefields because of their nanoscale heterogeneity. Individual forcefields can sufficiently describe specific structural moeities, such as short-range graphitization at low density [67], but are often not sufficient to capture the full range of non-ideal features observed in CDCs [68].
One possible way to increase the accuracy of computationally derived CDC structures is to employ sophisticated reactive forcefields, such as ReaxFF [69,70]. ReaxFF uniquely includes three-and four-body, non-bonded van der Waals, and electrostatic interactions in conjunction with a long-range bond order concept, where parameter sets are trained to reproduce a range of first-principles data describing reaction energies as well as reaction barriers. This allows ReaxFF to largely approach the accuracy of first-principles calculations, but at significantly (several orders of magnitude) reduced computational cost, although, we note more costly than other reactive and non-reactive force fields [52]. The use of ReaxFF for the formation of CDC structures is limited, with only a single recent study [71] examining the formation of amorphous carbon structures using QMD. This study examined the carbon hybridization, ring formation, and pore size distribution as a function of density and quench rate, and demonstrated the potential that ReaxFF could be effective in generating model CDCs; however, the study was not explicitly focused on the formation of experimentally relevant CDCs. EDIP, which has been shown to exhibit similar behavior to ReaxFF at low density [67], has been used to generate CDCs via removal of metal atoms from a carbide lattice and subsequent annealing at the target density, showing close agreement with experimental pore size distributions.
Here, we further explore the use of the ReaxFF force field to develop novel, structurally relevant atomistic models of CDCs. Our combined computational and experimental approach directly compares experimentally-derived structural information to properties of models derived from simulation. Here, we examine the effects of quench rate on CDC structure, including a slower rate than was used in prior ReaxFF studies [71]. We also investigate the role of a novel post-quenching compression step as a means to adjust pore size distribution. CDC structure is examined, including pore size distributions, pair distribution functions, and ring sizes. These analyses make comparisons to X-ray scattering, nitrogen sorption, and electron microscopy results. We demonstrate close agreement between experiment and simulation, including pore size distributions, for our approach.

Quenched Molecular Dynamics Simulations
Reactive molecular dynamics simulations were performed using the ReaxFF [69,72,73] force field implemented in the LAMMPS [74] simulation package (Version 27-Jul-2013, Sandia National Laboratory, Albuquerque, NM, USA). ReaxFF is a reactive potential that models the breaking and formation of chemical bonds through bond order, which is calculated from interatomic distances and used in calculating a bond energy. Three-and four-body, non-bonded van der Waals, and electrostatic interactions are also included. Parameter sets are developed by fitting to structural and thermodynamic properties to quantum density functional theory (DFT) data. The specific ReaxFF parameter set used in this study was optimized for modeling carbon materials [70].
Initial configurations were generated by assigning 20,000 carbon atoms to random positions in a periodic simulation cell of length 7.488 nm. This corresponded to a bulk density near 0.95 g/cm 3 , consistent with previous studies [49,55]. In an attempt to mitigate biases arising from an ordered initial state, random positions give the system a fluid-like initial state. Quenched molecular dynamics simulations were performed from an initial temperature of 3500 K to a final temperature of 3000 K. The use of temperature ranges higher than used in experiment is a common approach in ReaxFF studies to bring reaction kinetics into a time scale accessible in molecular simulation [70,75,76]. This temperature range is near the phase transition from liquid carbon to graphite at low pressures reported by experiment [77,78]. While the melting point was not rigorously examined with this force field, the transition appears to occur within the range expected from experiment. In simulations of the lengths considered here, major structural changes do not occur after the system has condensed from a fluid to solid state, so additional quenching below 3000 K was not considered. In order to explore the effects of quench rate, simulations varied from 5 ps to 500 ps in length, giving quench rates between 1 and 100 k/ps. These rates are studied for their effects on the structures and are not to be interpreted physically or compared directly to synthesis parameters. Temperature was controlled using a Nosé-Hoover thermostat with a damping constant of 10 fs. All quenching simulations were performed in the canonical (NVT) ensemble and a 0.5 fs time step was used throughout. This is a larger timestep than what is normally used in ReaxFF studies (0.25 fs); however, these systems include only carbon atoms and therefore lack any O-H or C-H vibrations, which are much faster than C-C vibrations.
For select samples, an additional compression step was applied after the quenching step using the NPT ensemble. The NPT integrator follows the work of Shinoda et al. [79], combining the work of Martyna et al. [80] and Parrinello and Rahman [81], as implemented in LAMMPS. The ReaxFF force field was again applied with a time step of 0.5 fs. Isothermal temperature control and isobaric pressure control at 3000 K and 20,000 atm were used with damping constants of 10 fs and 100 fs, respectively. Selection of this pressure is discussed in the Supplementary Material, including the relationship between system density and compression pressure in Figure S1. Like the quenching step, there is no direct physical analog to this step, but it was used to drive structural properties toward the experimental target. Effects of the post-quenching compression step are discussed in context below. For both compressed and non-compressed samples, no additional structural relaxation was performed.

Grand Canonical Monte Carlo Simulations
In order to study the properties of QMD-generated CDC models as adsorbent materials, grand canonical Monte Carlo (GCMC) simulations of nitrogen isotherms were performed with LAMMPS. A single-site nitrogen model was utilized as the fluid and both fluid-fluid and fluid-solid interactions were described with 12-6 Lennard-Jones potentials. A cut-off distance of 15 Å was used with no long-range corrections. Carbon and nitrogen atoms were described using C = 28.0 K, N = 95.2 K, σ C = 3.4 Å, and σ N = 3.75 Å. Lorentz-Berthelot mixing rules were used, giving C/N = 51.6 K and σ C/N = 3.575 Å. The carbon substrate was fixed in position; only moves involving nitrogen sites were considered. Construction of simulated isotherms for each carbon structure involved many short simulations of 2 × 10 6 translation attempts and 2 × 10 6 insertion/deletion attempts. At a fixed temperature of 77 K, a range of 40 pressures ranging from 1 × 10 −6 atm to 1 atm for the bulk fluid were considered. Beginning with a bare carbon structure at the minimum pressure of 1 × 10 −6 atm, the aforementioned short simulations were repeated at the same pressure until changes in the number of nitrogen atoms in the system were sufficiently close to zero to be considered at equilibrium. Subsequent pressure values incremented by a small amount and this process repeated until all pressures in the prescribed range were considered.

CDC Synthesis and Annealing
Porous CDC particles were synthesized using a previously described procedure [5]. Silicon carbide (SiC) particles (1.0 µm-5.0 µm diameter, 99.1+ % purity, Alfa Aesar, Haverhill, MA, USA) were placed into a quartz boat and loaded into a tube furnace. The furnace ramped up to 900 • C under flowing Ar gas. Once the system reached this temperature, pure Cl 2 gas flowed over the furnace at a rate of 0.5 L min −1 for 4 h and etched away the metal carbide using the following reaction [82]: Subsequently, the system cooled down to 600 • C and H 2 gas annealed the leftover carbon particles for 2 h [26]. Finally, the material cooled down to room temperature under flowing Ar.
The material was placed into a graphite crucible and loaded into a vacuum furnace (Solar Atmospheres, Souderton, PA, USA). The furnace outgassed for 24 h until a high vacuum (1 × 10 −6 torr) was reached [7]. Under continuously pulled vacuum, the furnace ramped up at a 10 • C min −1 rate up to 700 • C, held at that isothermal condition for 11 h, and, subsequently, cooled to room temperature [83]. This step removed all non-carbon functional groups, impurities, and moieties (adsorbed from air exposure or byproducts of CDC synthesis) from the CDC surface.

Gas Sorption and Porosimetry
Gas sorption measurements were carried out using Quadrasorb sorption analyzer (Quantachrome Instruments, Boynton Beach, FL, USA). Samples (0.050 g-0.075 g) were loaded into glass cells and outgassed at 120 • C for 24 h prior to analysis. The instrument collected gas adsorption and desorption isotherms (adsorbed volume vs. relative pressure) using N 2 adsorbate gas in the 0.00075 P P 0 -0.9995 P P 0 range (P 0 = 760 torr). Measurements were collected at a 77 K isothermal condition using a liquid N 2 coolant bath. The Brunauer-Emmett-Teller equation (BET) was used to calculate the BET specific surface area (SSA) of the CDCs [84]. The equation was calculated using a linear regression for 0.05-0.30 P P 0 values and normalized by the sample mass to obtain m 2 g −1 values. Quenched Solid Density Functional Theory (QSDFT) data reduction analysis provided information on a secondary (DFT-derived) SSA value, cumulative pore volume cm 3 g −1 values), and pore size distributions (PSDs, pore diameter vs. dV dr plots) [30]. The data reduction used only the adsorption isotherms and assumed slit pore configurations [85]. Quantachrome's Quadrawin software (Version 5, Quantachrome Instruments, Boynton Beach, FL, USA) provided all BET calculations and DFT data modeling kernels and computational analysis.

Ex Situ X-Ray Total Scattering
Synchrotron X-ray total scattering data were collected at the 11-ID-B beamline at the Advanced Photon Source (APS) at Argonne National Laboratory (Lemont, IL, USA). Ground powder samples were loaded into polyimide capillaries (1 mm inner diameter) and measured in transmission mode at ambient conditions (photon energy 58.65 keV, 0.2114 Å) using an amorphous silicon image plate detector (Perkin Elmer, Fremont, CA, USA). Scattering data of a CDCs annealed at 700 • C were collected for a total of 15 min and parasitic scattering from an empty polyimide container (duplicating the assembly without a sample) and an Ni powder standard (99.99%, Sigma Aldrich, Saint Louis, MO, USA) were collected using the same setup. The program Fit2D [86] was used to calibrate the sample to detector distance and detector alignment with data from a CeO 2 powder standard. Raw scattering data was integrated into Q-space spectra, applying a mask and polarization correction during integration. The normalized total scattering patterns, S(Q) were produced in PDFgetX2 [87] (Version 1.0, Michigan State University, East Lansing, MI, USA) by subtracting polyimide container scattering, utilizing the appropriate sample composition, and applying standard corrections for the area detector setup [88]. Pair distribution function patterns, PDF G(r), were calculated via Fourier transformation of the total scattering data utilizing a Q range of 1.0 Å −1 (Q min ) to 26.0 Å −1 (Q max ) for the Ni dataset, and was fit in PDFgui [89] in the range 1 Å-50 Å to calibrate the instrument Q damp = 0.038(1) and Q broad = 0.021(1) effects on real-space data. The same Q range was used for sample data. Because the diffractometer is only capable of capturing a small portion of such pore signals with the limitation of Q min near 0.5 Å −1 , a Fourier transform with Q min of 0.5 Å −1 will thus result in a weak oscillating signal with a wavelength of approximately 13 Å in the real-space domain from the 2π Q inverse relationship. For the purpose of comparing PDFs to QMD-generated models without complications arising from incomplete measurements of pore structures, a Q min value of 1.0 Å was selected.

Abberation Corrected Scanning Transmission Electron Microscopy
The structure of the CDC was characterized using an abberation corrected Nion UltraSTEM (Nion Company, Kirkland, WA, USA) operating at 60 kV. Scanning transmission electron microscopy (STEM) imaging provides [90,91] atomic resolution characterization of the overall CDC morphology following the chlorine heat treatment and subsequent vacuum annealing.

Results and Discussion
A quenched molecular dynamics routine was used to generate atomistic models for CDCs. Images provide qualitative analyses of sample morphology and structural properties, while quantitative characterization includes pair distribution functions, pore size distributions, ring size distributions, and ring neighbor correlations. Comparisons of QMD-generated CDC models to experimental results are made in appropriate corresponding discussions.
A summary of the models is listed in Table 1. Throughout, models generated through simulation named 'QMD-x' where 'x' is the quench rate used to produce it (in units of 1 K ps −1 -100 K ps −1 ). A 'c' is appended to the name if a post-quench compression step was included. The experimental sample is referred to as 'Expt.'. Table 1. Summary of sample structures. Models are named 'QMD-x' where 'x' is the quench rate used to produce it (in units of 1 K ps −1 -100 K ps −1 ). A 'c' is appended to the name if a post-quench compression step was included. Q: Quench rate, K ps −1 ;p: mean pore size, nm; L: Box length, nm; ρ: Particle density, g cm −3 .

Imaging and Spectroscopy of CDC Structures
Before evaluating QMD-generated CDC structures, we first consider imaging of these structures from experiments. Figure 1 shows representative annular dark-field (ADF) STEM images of the CDC structure. They reveal that the CDC structure is composed primarily of high curvature graphene sheets.
The STEM images also show regions where the graphene sheets come together along their edges to form slit-like geometrical configuration. Moreover, due to the atomic-scale chemical sensitivity of ADF STEM imaging, single atoms with an atomic number greater than carbon are clearly visible in the CDC. Shown in Figure S2, Electron Energy Loss Spectroscopy (EELS) measurements from these brighter atoms confirm that these are silicon atoms, likely the result of incomplete chlorine treatment. Raman spectra were also collected on the experimental CDC sample and presented in Figure S3. Prominent D and G bands, characteristic of disordered carbon, are observed.

QMD-Generated CDC Structures
Snapshots of final-state structures produced with QMD are shown in Figure 2 and can be compared to representative STEM images of experimental samples shown in Figure 1. Simulation boxes are cubic with characteristic lengths near 7 nm. Therefore, while mesoscale features cannot be analyzed, properties on the order of few Å can be considered. In particular, we focus on the aggregation of hexagonal and non-hexagonal rings into sheet-like networks and the size distribution of pores. These properties are quantified in  Clear trends with respect to quench rate are observed: faster quenches produce more amorphous structures and slower quenches produce more ordered structures. This is in agreement with prior QMD studies of amorphous carbons [53,55,71]. Shown in Figure 2a, a quench rate of 100 K ps −1 generates a highly amorphous structure with many stringy domains and few rings. Only a small number of rings are found, as quantified in Figure 4, and no specific ring size dominates the distributions. This low quantity prevents the formation of large sheets with graphitic features, which had been observed in CDCs synthesized at high temperature and that had exhibited prominent graphitic domains. The pore size distribution (Figure 3) shows that all porosity is below 1 nm and most pores are approximately 0.5 nm. Despite this, the observed nature of the carbon bonding indicates that these nanopores are not necessarily directly representative of those found in synthesized CDCs.
Slowing the quench rate gives more time for rings to form, as seen in Figure 2b, and some networks begin to reach the size of several rings and resemble small graphene-like sheets. Rings are mostly hexagonal, but a number of Stone-Wales and similar defects are observed. Stringy domains are not eliminated but appear at notably smaller frequency. The overall structure remains amorphous, but the greater quality of ring networks generates pores of greater definition. These pores are also larger, as seen in Figure 3. Many pores are sub-nanometer, but there is also a broad shoulder between 1 and 2 nm.
As shown in Figure 2c, further reducing the quench rate by a factor of 10 to 1 K ps −1 generates a more ordered structure. Most carbon atoms assemble into hexagonal rings, which aggregate into large, curved sheets that enclose nanopores. Some sheets collapse onto each other and form domains reminiscent of [002] stacking in graphite, but most sheets are single-layer and comprise pore walls. These pores exhibit a wider size distribution. While a fraction of pores is in the subnanometer diameter range, Figure 3 shows similar quantities of pores up to 2.5 nm in size. This emergence of larger pores with decreasing quench rate has been observed in previous QMD studies [55,71,92]. It bears repeating that quench rate in simulation lacks a direct physical analog; therefore, this trend cannot directly be attributed to any features of synthesis via halogen etching despite how physically realistic the generated structures may be. In summary, slower quenches produce structures with physically realistic carbon bonding but larger pores. Because the small pore size of CDCs gives rise to many important physical properties, it is worth attempting to modify the structures generated from slow quenches to have smaller pores. To this end, we applied a post-quench compression step to the QMD-10 and QMD-1 structures via NPT simulations at 20,000 atm and 3000 K, generating structures QMD-10c ( Figure 2d) and QMD-1c (Figure 2e).
Comparisons between Figure 2b,d and Figure 2c,e shows that the compression step qualitatively decreases pore sizes, particularly of larger pores. Some rearrangement is observed in Figure 2d, but for the most part, and particularly in Figure 2e, pores are modified without collapsing existing pores and/or forming new ones. Qualitative visual observations do not show any alterations of ring bonding. Thus, the combination of ReaxFF with a post-compression step results in pore size distributions that more closely match experiments as well as a recent study that generated CDCs using the EDIP force field and, rather than QMD, annealing of a system after removal of metal atoms from a carbide lattice [65]. Specifically, this compression step combines properties of QMD-generated models at different quench rates: the local bonding structure is dominated by six-membered rings, but the pore size distribution does not largely extend in the mesoporous regime.

Pore Size Distribution
Pore size distributions of QMD samples are shown in Figure 3. Consistent with previous QMD studies [55,92], the pore size and dispersity increase with slower quench rates. The compression steps (samples QMD-10c and QMD-1c) are shown to quantitatively increase the population of sub-nanometer pores and, to a lesser extent, decrease the dispersity. Although the post-quenching compression step fails to fully reach the monodisperse sub-nanometer target, agreement with the experimental pore size distribution is improved, most notably for QMD-10c.
We note that, while experimental characterization results demonstrate a strong predominance of sub-nanometer pores, a small shoulder at 1-1.5 nm and a small peak between 2.0 and 2.5 nm are observed. As such, the larger pore size distributions observed in some of our or systems may indeed be realistic. However, they may form in higher frequency than in experimentally produced CDCs.

Ring Size Distributions
Although a robust quantification of the number of non-hexagonal rings is not currently accessible via empirical characterization of experimental systems, it can be quantified in simulation. The relative quantities of rings of sizes 4 through 8 in QMD-generated structures are shown in Figure 4. Ring counts are normalized to the number of hexagonal rings in a fully graphitic sample of the same size (10,000 rings in these 20,000 atom systems). By comparing compressed and non-compressed samples, the compression step is shown to have a small impact on the number of rings but little impact on the distribution. On the left are ring size distributions from samples produced by QMD alone and on the right are ring size distributions from samples that also underwent a post-quench compression step.
The quench rate is shown to have a significant impact on the quantity of ring formation and type of rings formed. At the fastest quench rate, few rings are formed with no particular size dominating the distribution. At the slowest quench rate, the total number of rings approaches the graphitic limit and hexagonal rings are the most common. The ring size distribution for the samples of medium quench rate mixes properties of the limits: the total number of rings is large and hexagonal rings are the most common, but not by a large margin. For all QMD-generated structures, the number of four-membered rings is negligible and the number of eight-membered rings is non-zero but not significant. For all quench rates, the number of five-and seven-membered rings is similar, indicating that their existence may be structurally coupled.
To further investigate the relationships between non-hexagonal rings, neighbor-neighbor correlations of non-hexagonal rings were analyzed and shown in Figure 5. For example, the number of 7-membered rings neighboring a 5-membered rings was counted and shown as a '5-7' neighbor type. With slower quench rate, 5-7 and 7-7 neighbor types increased in quantity while 5-5 neighbor types decreased in quantity. This trend is consistent with common graphene defects such as the Stone-Wales 55-77 defect [93] and the 555-777 divacancy [41]. Other common non-hexagonal defects found in graphene follow a consistent theme of 7-membered rings neighboring each other and 5-membered rings, but rarely 5-membered rings neighboring each other [94].

Sorptive Loading
Nitrogen isotherms generated from GCMC simulation are shown in Figure 6 and compared to a representative experimental isotherm. Characteristic Type I isotherms are observed: most loading occurs at low relative pressures and little change in loading occurs at moderate relative pressures. A partial condensation phase transition is only observed in the QMD-1 sample, which is expected to be the result of mesopores. The loading spikes in experimental samples near P = P vap are the result of larger non-porous features formed from the packing of the CDC nanoparticle and not the result of mesopores found in bulk CDC phases. The lower loadings in compressed samples compared to non-compressed samples is a result of an artificially high density, although the shapes of the isotherms are not noticeably affected.

X-Ray Total Scattering
Further structural analysis relies on a synchrotron X-ray source for total scattering/pair distribution function analysis. The S(Q) of an annealed CDC sample is shown in Figure 7 and compared to simulated data of an ideal graphene sheet of 10 × 10 unit cells. We note here that the large intensity uptake observed in the low Q region arises from the surface scattering and pore structures present in the CDCs. Comparisons to QMD-generated structures in real space are discussed below. Most of the peaks observed in the experimental sample correspond well to the diffraction pattern of the graphene structure, the exceptions being peaks at 1.4(1) Å −1 and 1.83(2) Å −1 (∼3.43 Å and ∼4.5 Å in real space). The 1.83(2) Å −1 peak corresponds to graphite (002) reflection typically observed at 26 deg 2θ with Cu-Kα radiation. This peak, however, is not completely resolved and sits alongside a broad diffuse peak centered at 1.4(1) Å −1 . Such features suggest the presence of poorly coordinated graphene sheets, to which this peak is attributed. In addition, as mentioned in the experimental section, the large intensity uptake at Q < 1 Å −1 by surface and/or pore structure scattering is omitted for subsequent PDF analysis. The PDF G(r) is related to the RDF g(r) as therefore, a PDF, unlike an RDF, may fall below zero for some values of r. The length scale over which PDFs provide strong resolution (a few Å), corresponds to carbon atoms separated by a small number of bonds. Pair distribution functions from QMD-generated structures are shown in Figure 8. The software package DISCUS [95] (Version 4.2, Michigan State University, East Lansing, MI, USA) used in these calculations in order to include the effects of some biases in the experimental apparatus, including Q max , Q damp , and Q broad . For the local carbon bonding comparisons, only a small isotropic thermal displacement factor is introduced (B iso = 0.2 Å 2 , approximately 0.05 Å mean-square displacement) because the QMD-generated structures are large enough to contain the desired static disorder. The overall intensity scale factor is also adjusted to compare with observed PDF. To maintain consistency, all simulated PDFs used the same parameters. The peaks of these PDFs correspond well to the distances between carbons in an idealized graphene-like structure and follows previous interpretation in which no strong correlation can be observed between different sheets and/or pore walls [96,97]. The first three peaks at 1.43(2), 2.45(2) and 2.83(2) Å correspond to neighbors in a hexagonal ring and peaks at 3.75(5), 4.28(5) and 4.97(8) Å correspond to the location of carbon atoms in neighboring rings. The well-ordered nature of these peaks indicates, as expected, that hexagonal rings are favored. However, significant peak asymmetry and broadening observed at the third and further peaks do indicate the presence of non-hexagonal rings and the resulting sheet curvature effects. Together, these define the extent of structural coherence (approximately 15 Å) in the experimental PDF data. Compared to simulated PDFs, slower quenches provide better agreement with experiment, though QMD-1 is not able to totally capture graphitic ordering at distances greater than ∼5 Å. Consistent in both experimental and simulated PDFs is the general feature in peak broadening and the lack of graphite [002] correlation near 3.43 Å, the latter indicating a lack of sheet-sheet stacking despite other graphitic properties.

Effects of Initial Density
We consider the effects of initial density on the final structure. The QMD-10c sample, which has an initial density of 0.95 g cm −3 , was considered as a reference. Simulations were performed with identical inputs with the exception of the initial density. Values of 0.8 and 0.9 g cm −3 were considered. After compression, these two systems reached densities of 1.56 and 1.55 g cm −3 , respectively, as compared to 1.467 g cm −3 for the system originally at 0.95 g cm −3 .
The pore size distributions of each of sample, before and after compression at 20,000 atm, are shown in Figure 9. Pore size distributions of samples compressed at other pressures are shown in Figure S4. For all un-compressed samples, the distribution is broad, covering the range 0.3 nm-2.0 nm. There is also a small increase in the mean pore size as the density decreases. The compressed samples, however, exhibit similar pore size distributions. Each has a peak near 0.5 nm and a shoulder that disappears near 1.0 nm. The uniformity of these distributions suggested that the initial density does not strongly impact the porosity of the structure after compression, for these relatively high density systems, as quench rate is likely to be the dominant factor for structural evolution.
We also consider the pair distribution functions of these samples, shown in Figure 10. For samples with a range of initial densities, compression has a small impact. Peak heights are not noticeably changed, but compression shifts each peak slightly to the left. This is attributed to the high pressure needed to increase the density of each sample and not differences in the nature of the carbon bonds. Comparisons of samples before and after compression show that density has no noticeable impact on the pair distribution functions. In summary, compression is shown to have negligibly small impacts on short-range carbon interactions while driving most of the pore size distribution below 1 nm.
We should note, however, that they are still quite small compared to experiments.  As shown in Figure 1, experimental samples include a variety of highly amorphous domains and features such as graphitic slit pores and single, isolated graphene sheets, in addition to other heterogeneous nanoscale domains. As a result, experimental samples will show localized variation in structure and density. Furthermore, experimental CDC samples could be thought of as a collection of individual structures comprising larger particles. The models generated here represent a separate collection of structures that may overlap with those found in experiments. Thus, prudent studies of the behavior of CDCs should consider the collective behavior of a range of nanoscale structures to better understand the resulting physical properties observed in experiments.

Conclusions
The use of the ReaxFF reactive force field in generating an atomistic model for carbide-derived carbons was explored using a quenched molecular dynamics (QMD) scheme. The effects of quench rate and system compression were explored. Comparisons were made to an experimental carbide-derived carbon (CDC) material with similar structural properties. These models that use ReaxFF and employ a post-quenching compression step provide clear improvement upon prior CDC models generated using QMD in terms of local bonding and pore size distribution and also closely match those recently published using a process involving removal of metal atoms from a carbide and subsequent annealing [65].
Slower quench rates and a compression step yield good agreement with experimental pore size distributions (PSDs), pair distribution functions (PDFs), and qualitative morphological features in abberation-corrected scanning transmission electron microscopy (STEM) images. The sorptive behavior closely matches CDCs and typical nanoporous adsorbents. The prevalence of non-hexagonal rings and their neighbor correlations are also explored. System compression, after the quenching step, is shown to decrease average pore size toward a sub-nanometer distribution without disrupting short-range features such as ring size. Overall, QMD-10c and QMD-1c are presented as viable atomistic models for CDCs.
Future approaches may further improve the quality of this model. Specifically, an effort that combines a QMD routine with a reverse Monte Carlo stage could improve fits to PDFs and, if a suitable error norm can be constructed, also generate more realistic PSDs. The development of highly parallelized Monte Carlo software would greatly aid in this, as the computational cost of ReaxFF or other complex potentials necessitates parallelization. If future experimental probes more accurately quantify the distribution of ring sizes, they will provide another set of target data to which a simulation routine could be tuned.
This work provides the framework for use of these structures as model adsorbents and electrodes in molecular simulation of adsorption and electrochemical capacitors, respectively. Specifically, use of these models introduces non-idealities such as pore size dispersity, surface curvature, and non-hexagonal rings.
Supplementary Materials: The following are available online at http://www.mdpi.com/2311-5629/3/4/32/s1. Figure S1: System densities over the length of compression simulations with different reference pressures; Figure  S2: ADF STEM image and EEL spectra collected on a SiC-CDC sample after vacuum annealing at 700 • C; Figure  S3: Raman spectrum of experimental sample; Figure S4: Pore size distributions resulting from the use of barostats with high reference pressure.