Computer Simulation of Carbonization and Graphitization of Coal

This study describes computer simulations of carbonization and graphite formation, including the effects of hydrogen, nitrogen, oxygen, and sulfur. We introduce a novel technique to simulate carbonization,"Simulation of Thermal Emission of Atoms and Molecules (STEAM),"designed to elucidate the removal of volatiles and density variations in carbonization residue. The investigation extensively analyzes the functional groups that endure through high-temperature carbonization and examines the graphitization processes in carbon-rich materials containing non-carbon"impurity elements". The physical, vibrational, and electronic attributes of impure amorphous graphite are analyzed, and the impact of nitrogen on electronic conduction is investigated, revealing its substitutional integration into the sp$^2$ layered network.


Introduction
The global graphite shortage has steered research towards alternatives to natural graphite for current and emerging technologies [1].Among these materials, coal has garnered natural attention.Although the detrimental effects of coal combustion for energy generation on health, climate, and the environment have been extensively documented in previous studies [2,3,4], contemporary apprehensions regarding the scarcity of natural and synthetic graphite (currently on the list of critical materials [5]) have prompted a fresh appraisal of coal's potential as a carbon resource [6,7,8].In the United States, research on coal utilization has extended beyond energy production to include using coal-based cryptocrystalline graphite as electrode materials in lithium-ion [9,10,11] and aluminum-ion [12] batteries.Additionally, there is a growing interest in using coal as filler material in composites to enhance the mechanical properties of polymers [13], as well as to improving the electrical properties of metals [14].The later exemplified by successful cases using polycrystalline graphene composites with copper [15] and aluminum [16].
Synthetic graphitization of coal begins with carbonizationthe controlled heating of coal within the temperature range of 800 K to 1500 K, in the absence of oxygen.This elevated temperature triggers the release of volatile compounds including gaseous hydrocarbons, oxides of carbon, organo-nitrides, and organo-sulfides, yielding a solid residue known as coke [17].
The consensus is that subjecting the coke to higher temperatures (2500 K ∼ 3200 K) induces the cleavage of aliphatic chains and the formation of polyaromatic compounds through free radical crosslinking [18].These polyaromatic compounds exhibit interconnected carbon atom rings, resulting in structures having sp 2 hybridization, akin to graphite [19].However, research suggests that non-carbon elements in graphitic precursors (e.g., coal) such as hydrogen, nitrogen, oxygen, and sulfur, remain present in the final graphite material and interfere with its formation [20,21,22].For instance, Franklin demonstrated that the presence of oxygen can render the graphite precursor nongraphitizing [23].Bourrat and co-workers showed that sulfur can endure in the coke after carbonization of some graphite precursors for temperatures above 2000 K, forming crosslinks between the graphitic layers [24].
Established experimental techniques such as Raman Spectroscopy [25], Fourier transform infrared spectroscopy [26], electron paramagnetic resonance [27], and others, can inform on the structure of the materials obtained at different stages from coal carbonization to its graphitization, but they fall short in clarifying how low concentrations of non-carbon atoms influence properties like the electronic conduction in coal-derived nanostructures.
Molecular simulation is an ideal method to gain insight into the atomistic attributes of coal and its derived nanostructures.Recent studies led by Thapa and colleagues [28,29] have provided a new perspective on the graphitization mechanism of carbonaceous materials, painting an atomistic portrait of amorphous graphite formation at elevated temperatures in elemental carbon.This work challenged the notion that layering in graphite emerged solely from Van der Waals forces, and indicated that layering is due significantly to wavefunction mixing, involving interactions among π electrons in the galleries.This insight, also observed in other layered carbon allotropes such as multi-shell fullerenes [30] and multi-walled nanotubes [31], emphasizes the importance of molecular simulations to understand the chemistry and bonding in complex systems.
Our study examines carbonization and the emergence of layered nanostructures derived from coal-like computer models.We include carbon, hydrogen, oxygen, nitrogen, and sulfur.A new method for simulating carbonization, known as the "Simulation of Thermal Emission of Atoms and Molecules (STEAM)," is presented and explored in detail.Analysis of gas evolution from bituminous coal models during carbonization, and the functional groups present in the coke are conducted.Additionally, ab initio density functional theory (DFT) [32] is employed to probe graphitization in models with "coal-like" elemental compositions containing 5% and 10% non-carbon constituents (by weight).The analysis encompasses an exploration of the vibrational modes, electronic structure, and transport in impure amorphous graphite.Furthermore, the study investigates the impact of carbon-substituted nitrogen on the electronic transport in the layered amorphous nanostructure, utilizing the space-projected conductivity technique [33,34].

Computational Methods
Molecular dynamics (MD) simulations were conducted using two methods: (1) The reactive force field (REAXFF) interatomic potential [35], including hydrogen, carbon, nitrogen, oxygen, and sulfur, and (2) DFT with plane-wave pseudopotentials for all the aforementioned elements except hydrogen.The REAXFF calculations were executed within the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [36], while the Vienna Ab initio Simulation Package (VASP) [37] was employed for DFT computations.A single k-point (Γ) and periodic boundary conditions were used for all DFT calculations.When relevant, simulations followed the canonical or the isobaric-isothermal ensemble at specified temperature and/or pressure, controlled by a Nosé-Hoover thermostat and/or barostat [38,39,40].MD simulations using REAXFF and DFT employed timesteps of 0.25 fs and 1.5 fs, respectively.The subsequent sections elaborate on the specific simulation protocols for carbonization and graphitization.

Carbonization
An initial model of coal was made.
Starting with ChemDraw ® , 3D computer models of Pittsburgh No. 8 coal (hereafter referred to as P8 coal) were constructed from Solomon's proposed P8 coal, derived from structural and thermal decomposition data [41].P8 coal is a metallurgical, highvolatile (A) bituminous ranked coal (hvAb) [42].Figure 1a showcases the 2D representation of the P8 Solomon model, characterized by its chemical formula: C 161 H 140 O 15 N 2 S 2 .
The aliphatic and aromatic structures were initially delineated in 2D using ChemDraw ® .Subsequently, these structures were subjected to crude optimization through the MM3 force field [43] to achieve a stable steric conformation (refer to Figure 1b).The PACKMOL software package [44] was used to construct supercell models, each comprising 5 units (1600 atoms) of the macromolecular unit structure.In the figures, hydrogen, carbon, nitrogen, oxygen, and sulfur atoms are visually distinguished by white, gray, blue, red, and yellow colors, respectively.The supercell models were extensively annealed using REAXFF to optimize their structure.New configurations were achieved through conjugate gradient relaxation, leading to local energy minima.To obtain realistic 3D coal models that reflect the density of P8 coal, the supercells were further annealed under isothermal-isobaric conditions (NPT), allowing density optimization.By doing so, the resulting models satisfy the periodic boundary conditions and therefore, can be used in accurate DFT calculations.
Vitrinite reflectance (R o ) is a well-established gauge of the interplay between temperature and pressure during coal maturation [45].Utilizing R o = 0.78 for P8 coal [46] translates to a corresponding burial temperature of around 390 K [47,48].Consequently, the NPT simulation maintained a constant temperature of 400 K. Low pressure values (0.2, 0.4, and 2.0 GPa) were selected in alignment with experimental data from Reference [45], where a similar pressure range was used to establish the relationship between R o and pressure.
Figure 2a shows the density time evolution for the simulations at 0.2, 0.4, and 2.0 GPa.Other relevant thermodynamic quantities were also continuously monitored to ensure simulation stability.The energy and pressure time-series plot at 0.4 GPa is depicted in Figure 2b.Models obtained at 2.0 GPa and 0.4 GPa were selected for this study due to their optimal densities (ρ) which align closely with experimental (1.46 g/cm 3 ) and estimated particle (1.22 g/cm 3 ) densities for Pittsburgh No. 8 coal, as reported by White and co-workers [42].The optimal densities obtained at 0.4 and 2.0 GPa are shown in the second column of Table 1.
The models were equilibrated at 500 K for 50 ps under constant temperature (NVT) conditions.Subsequently, the tem- perature was ramped to 1273 K over 100 ps.We introduce a new simulation protocol designed to capture the emission of volatiles and variations in density during coke formation in the carbonization process.This technique, referred to as "Simulation of Thermal Emission of Atoms and Molecules (STEAM)," is an iterative approach that incorporates successive NVT and NPT cycles.During the NVT phase, volatile molecules within a specified molecular mass range are systematically removed become non-bonded to the rest of the network.Meanwhile, the NPT phase ensures density convergence and maintains periodic boundary conditions in the system.
In our implementation of STEAM within this study, both the NPT and NVT phases extended over 125 ps each.Throughout the NVT phase, molecules with a molecular mass below 50 g/mol -formed due to bond cleavage -were removed at a rate of up to 5 molecules every 1.25 ps -modeling the outgassing.Completion of the STEAM process was determined by three criteria: (1) Absence of atom/molecule deletion in an NVT step, (2) predominance of high-molecular-mass molecules (preferably a substantial coke-forming molecule), and (3) maintenance of consistent pressure during NPT steps after the NVT phase in (1).Our simulation protocol is summarized as follows: 1. Formulate coal model and generate supercell of macromolecular coal units.Figure 2c shows coal structure formed at 0.4 GPa (ρ = 1.21 g/cm 3 ) after item 3 was completed, and the coke structure obtained after item 10.A corresponding figure for the coal and coke formed at 2.0 GPa can be found in Figure S1 in the supplementary material.
In the initial phases of carbonization, rapid emission of light gases was observed, succeeded by expulsion of heavier gases before the eventual formation of coke.This pattern is illustrated by the density-time plots for the NVT phase in the overview of the STEAM process at 2.0 GPa and 0.4 GPa, shown in Figure 3 and Figure S2 respectively.Within the NVT step, the consistent volume scaling demonstrates that the release of low-mass molecules corresponds to minor step heights, while highermass molecules result in more prominent step heights.Noteworthy among the gases emitted are hydrogen (H 2 , stemming  Subsequent bond cleavage analysis revealed that the primary source of oxygen-related volatiles was hydroxyl (OH) present in carboxyl (R-COOH) or alcohol (R-OH) functional groups.In specific cases, ether (R-O-R ′ ) bond cleavage form formaldehyde (HCHO), which subsequently combined with •OH radicals to produce CO and CO 2 .As carbonization progressed, prior to coke formation, hydrogen sulfide (H 2 S) emerge from mercaptans (R-SH).The delayed release of H 2 S was also observed by Whittaker and Grindstaff [22]: The evolution of the sulfurous gas induces an internal pressure that results in the irreversible expansion of graphitic precursors in the late stages of carbonization, or the early stages of graphitization [22].In rare instances, CS (with potential for evolution into CS 2 ) was liberated from heterocyclic thiophenes.Furthermore, bond fragmentation within aromatic rings resulted in aromatic hydrogen generation and triggered the liberation of hydrogen cyanide from ring nitrogen.
The compositions of the coke, as detailed in Table 1, revealed that over 65% of the non-carbon elements were emitted as gases, contrasting with the retention of about 80% of carbon atoms.As depicted in Figure 2c and Figure S1, some carbon ring structures (5, 6, and 7-membered rings) were retained in the coke (indicative of the early stages of graphitization).The coke matrix retained pyrrolic and pyridinic nitrogen due to hydrogen atom removal from 5-and 6-membered heterocyclic aromatic rings, respectively.Aliphatic ethers (R-O-R ′ ) also persisted, serving as bridges between aromatic and aliphatic carbon structures.Cyclic ethers also formed as oxirane and oxolane structures.The Carbonyl components in P8 coal gave rise to emerging ketones (R 2 C=O), and the Organosulfur compounds manifested as thioethers (R-S-R ′ ) and thioketones (R 2 C=S).

Graphitization
To study graphitization, we found that utilizing DFT precision forces is required [28,29]; thus, in this section we employed VASP.Our approach involves randomly distributing 200 atoms of carbon, nitrogen, oxygen, and sulfur in a cubic box to achieve a density of 2.4 g/cm 3 .Henceforth the non-carbon elements will be referred to as "impurity elements".Six models were created, incorporating 5% and 10% impurity elements (three models for each concentration).The oxygen to nitrogen and oxygen to sulfur ratio was 3:1 for both cases.This ratio was deliberately selected to reflect the elemental concentration in the coke (as detailed in Table 1), while maintaining a carbon-rich environment.The models were annealed for 360 ps at 3000 K -the common graphitization temperature [49,50,51].Subsequently, the models were relaxed to an energy minimum configuration using the conjugate gradient algorithm.Note that hydrogen is excluded.This is because hydrogen is typically "burned off" in the carbonization phase, and frankly, the timestep required is too short for practical VASP simulations.
The interactions between the electrons and ions were described using the Perdew-Burke-Ernzerhof projectedaugmented-wave pseudopotential [52,53].We set a cutoff energy of 400 eV for the plane-wave-basis used to expand electronic wave functions during the molecular dynamics simulation, and a higher cutoff of 520 eV was used for the electronic structure calculations.For where it applies, the layered nanostructures containing non-carbon elements are referred to as impure amorphous graphite.
Figure 4a -d is a chronological depiction of the layering process observed in layered nanostructure with 5% impurity concentration.In contrast, the models with 10% impurity concentration lack distinct layers; instead, they exhibit inter-layer atomic bonding involving carbon atoms near oxygen atoms (see Figure 4e).Wang and Hu highlighted the connection between layering defects in graphite and the presence of oxygen, revealing that even at low concentrations, graphite oxidation significantly disrupts its layered structure [54].
Preceding any discernible indications of layering (indicated by the brown arrows in Figure 4b), a process involving rearrangement of the non-carbon atoms occurred.Nitrogen atoms readily adopt graphitic nitrogen (N-3) structures (forming bonds with three neighboring carbons in sp 2 configuration).Once formed, this configuration was maintained throughout the simulation.In specific cases, nitrogen atoms disrupt ring connectivity in the layers by bonding with only 2 carbon atoms after substituting a sp 2 carbon atom, thereby forming pyrrolic (N-5) or pyridinic (N-6) structure by carbon substitution in a 5-or 6-membered ring, respectively.Oxygen atom form ether-bridges (C-O-C) or ketone (C=O) structures that terminate ring connectivity.Similarly, sulfur atoms display a preference for thioether (C-S-C) and thioketone (C=S) structures.Notably, sulfur conformation in impure amorphous graphite mirrors sulfur's inclimation to form C-S-C bond at the edges (grain boundaries) in crystalline graphite [55].
In early graphitization stages, sulfur atoms initially establish crosslinks between layers (see Figure 4c).However, upon energy optimization, sulfur atoms demonstrate a preference for bonding within the layers (refer to Figure 4d).From experimental observations, Kipling and co-workers suggested that crosslinking of sulfur between layers exist in graphitized materials from sulfur-containing precursors [56] This was later challenged by Adjizian and co-workers in a DFT study that found such sulfuric crosslinking structures were energetically unfavorable in graphite [55].Our results offer a nuanced perspective: crosslinking does occur during intermediate graphitization stages, but the energy-optimized sulfur conformation is not as a crosslink between the layers.Additionally, our observation from the structures with 10% non-carbon elements also contradicts Kipling's notion that higher sulfur concentrations could promote inter-layer connections; instead, distinct layers simply failed to form, atleast on the timescale of our simulation.
The functional groups in impure amorphous graphite were also identified in coke post-carbonization (see Section 3.1), indicative of their energy stability within carbon layers during coal graphitization [20].This notion is reinforced by NMR analysis of graphite oxide, revealing stable ether (C-O-C) structures and unstable hydroxyl (C-OH) groups, with the transient C-OH condensing into consolidated C-O-C linkages [57].Furthermore, a predominant sulfur-doped graphite conformation reported by Li et al. is the theophinic (C-S-C) structure [58].
We analyzed the local structure in impure amorphous graphite using the radial distribution function (Figure 4f and g).Key details about the primary peak, which represents the nearest distance between carbon and non-carbon pairs, are provided in Table 2. Our analysis is limited to this range due to the low impurity concentration in the models, which inherently lacks extensive structural information beyond immediate neighboring atoms.Notably, the nearest neighbor values derived from our models closely correlate with experimental data for the observed functional groups, as discussed earlier.Another feature is the emergence of a distinct peak around 1.24 Å in the radial distribution function plot of the 5% impurity model (Figure 4f), corresponding to the C=O bond length in carbonyl groups [59].This peak becomes more prominent in the 10% impurity model's plot (highlighted by gray dashed lines in Figure 4g), suggesting a higher oxygen concentration encourages the formation of carbonyl functional groups in impure amorphous graphite.The vibrational signatures in impure amorphous graphite were extracted and compared with that of pristine graphite [28] by computing the vibrational density of states (VDoS) using the harmonic approximation.This entailed calculating the Hessian matrix by evaluating forces from atomic displacements of 0.015 Å along six directions (±x, ±y, ±z).The VDoS (g(ω)) is computed as follows: where δ(ω − ω i ) is approximated as a Gaussian with a standard deviation of 1.5% of the maximum frequency observed.N and ω i represent the number of atoms and the eigen-frequencies of normal modes, respectively.The VDoS for the pristine and impure amorphous graphite have similar structures but with shifts in peak location (refer to Figure 5a).A shoulder appearing at the low-frequency end (around 5 THz) matched in both models.Moreover, the peak with a shoulder around 21 THz in the pristine graphite model was resolved as a single peak in impure amorphous graphite.While there are noticeable similarities in the vibrational structures between pristine graphite and impure amorphous graphite, it's important to recognize that the classification of vibration modes applied to pristine graphite, including acoustic and optical modes, cannot be directly extended to describe amorphous systems.To address this disparity, we rely on the concept of the phase quotient (Q p ) [64], which acts as a metric to discern whether vibrations can be characterized as acoustic (+Q p ) or optical (-Q p ) modes.The phase quotient is derived as [65]: where N b is the number of valance bonds, and u i p and u j p are the normalized displacement vectors for the p th normal mode.The index, i, ranges across all carbon atoms, and j enumerates neighboring atoms (C, O, N, S) linked to the i th carbon atom.A purely acoustic (optical) vibration manifests as a phase quotient of +1 (-1).
In Figure 5b, Q p profiles are illustrated as crimson and black curves, representing carbon/carbon (C-C) and carbon/noncarbon (C-X) vibrations.Particularly noteworthy is the stronger in-phase (acoustic mode) coherence of C-X vibrations at very low frequencies, in contrast to C-C vibrations.This coherence in C-X vibrations transitions rapidly to an out-of-phase pattern -akin to optical modes -even surpassing the out-ofphase behavior of C-C vibrations at higher frequencies.
In graphite, the acoustic mode is commonly linked with the low-frequency region < 26.8 THz [67].The Q p for C-C vibration in impure amorphous graphite shifted to the optical mode (-Q p ) at ≈ 21 THz.However, a local minimum emerged at Q p = -0.15,corresponding to a frequency of 24.4 THz (indicated by the first black line in Figure 5b).The vibration mode briefly reverted to the acoustic mode, peaking at 27 THz (the vibration mode-switching frequency in graphite), before Q p transitioned back towards the optical mode.Q p became negative (optical mode) at ≈ 34 THz, indicated by the second black line in Figure 5b.Two potential explanations underlie this behavior.Firstly, the turning point of the phase quotient for C-X, occuring at a lower frequency of 22.6 THz, may have influenced the C-C vibrational mode, inducing a switch from optical to acoustic vibrational behavior.Secondly, the region experiencing a shift in slope within the phase quotient plot corresponds to the transition from bending to stretching vibrational characteristics.This transition is evaluated using the bending/stretching quotient (S p ) proposed by Marinov and Zotov [68]: Here, u i p and u i p are the same as in Equation 2, and ri j is the unit vector parallel to the m th bond.S p → 0 indicates bond-bending character, and S p → 1 represent predominantly bond-stretching vibrations.The region between 24.4 -34 THz, marked by two black vertical lines in Figure 5c, signifies a transition from bending to stretching vibrational character.Significantly, this transition aligns with the region of C-C vibration mode switching, as indicated by Q p (see Figure 5b), implying a direct influence of vibrational character transition on vibrational modes in impure amorphous graphite.
In Figure 6a-c, the electronic density of states (EDoS) is provided for three configurations: (a) amorphous graphite, and impure amorphous graphite with (b) 5% and (c) 10% impurity concentrations.The Fermi level (E f ) in the plots have been shifted to zero (indicated by the green dashed line), and only regions, E f ±2 eV, are displayed here.The complete electronic structures, plotted in Figure S3a, indicate minimal alterations to the overall electronic structure of amorphous graphite, when compared to the impure amorphous graphite in Figure S3b and c for 5% and 10% impurity concentrations.We use the inverse participation ration (IPR) to study localization of the electronic states.The IPR (ζ n ) is defined as: where γ i n represents the contribution of the i th Kohn-Sham state to a given eigenvector (ϵ n ).A high IPR value (ζ n → 1) indicates that the wave function is localized on very few atoms.Conversely, a low value ( ζ n → 0) indicates that the wave function is delocalized (distributed over many atoms).
The gray dashed lines in Figure 6a -c represent the IPR for the electronic bands.In amorphous graphite, the low IPR suggests the absence of localized states near E f .However, the introduction of impurity elements leads to localized states near E f (see Figure 6b and c).
We identified the elemental species responsible for the highly localized states in impure amorphous graphite, focusing on IPR  E f ) is shifted to zero and indicated by green dashed lines, and only the region between E f ± 2 eV is shown.The complete EDoS is provided in Figure S3.SPC calculation [66], presented as a grayscale heatmap on selected layers containing nitrogen, is shown in (e); the darker regions represent the electronic conduction paths in the layered nanostructures.The average electronic conductivity (σ, in S/m) for all 3 models is plotted in (f).values ≳ 0.5.In the layered nanostructure with 5% impurity concentration, two such states emerge, separated by 0.07 eV, and they are exclusively associated with a single ketone oxygen atom.When the impurity concentration is increased to 10%, additional localized states surface.In Figure 6c, states 1, 2, and 3 (marked by brown arrows) localize on distinct oxygen atoms forming ketones, while state 4 localizes on a thioketone.State 5, with IPR ≈ 0.48, is also centered on the same thioketone as state 4, albeit separated by 0.2 eV We note that the nitrogen atoms in impure amorphous graphite do not contribute to the localized states around the Fermi energy (E f ).Extensive research has focused on intentionally doping nitrogen into various carbon allotropes [69,70].For example, nitrogen has been introduced into graphite (or graphene) to create nitrogen-doped quantum dots [71], and into diamond to form nitrogen-vacancy centers for potential supercomputing applications [72].With this context, we computed the electronic structure for impure amorphous graphite containing 5 weight-percent of only nitrogen atoms, as shown in Figure 6d and Figure S3d.The models were constructed by repeating the simulation protocol detailed in Section 3.2, and a representation of the nitrogen-containing impure amorphous graphite model is provided in Figure S3e.
As shown in Figure S3d, the overall electronic structure of the layered structure with 5% nitrogen, exhibits minimal differences when compared to that of amorphous graphite (see Figure S3a).Notably, two states emerge in proximity to the Fermi level (refer to Figure 6d).However, these states (with IPR values of 0.32 and 0.51) do not exclusively localize on specific atoms; instead, they are distributed across a small number of atoms, implying an absence of true localization.This observation is further illustrated in Figure S3f for the state with the higher IPR value (IPR = 0.51).
Interestingly, while nitrogen atoms substitute for carbon atoms within the sp 2 configuration of the layers, they do not contribute to localizing electronic states near E f .This leads to the question: Does nitrogen merely emulate carbon behavior within the layers, and how does it impact electronic conductivity in this layered nanostructure?To address this, we compute the electronic conduction and transport in the nitrogencontaining impure amorphous graphite, utilizing the spaceprojected conductivity (SPC) methodology [66].
The SPC framework utilizes the Kubo-Greenwood formula for electronic conductivity [73,74] and the Kohn-Sham singleparticle states (ψ i,k ) to project the electronic conductivity from a given atomic arrangement onto a spatial grid (see Reference [66] for details on the SPC method).Figure 6e depicts the electronic transport path using a grayscale heatmap within selected nitrogen-containing regions.Darker shades indicate high electronic conduction, while lighter areas indicate limited or absent electronic transport.Notably, the presence of nitrogen atoms along a conduction pathway acts to impede or completely halt electronic transport.To examine the potential impact of structural attributes on electronic conductivity originating from nitrogen within the layered nanostructure, we computed the radial distribution function (RDF) and bond angle distribution (BAD) for nitrogen-containing impure amorphous graphite (Figure S4a and b, respectively).The most prominent RDF peak, representing the average nearest neighbor distance calculated for carboncarbon (C-C) and carbon-nitrogen (C-N) bonds, were closely aligned -separated by only 0.015 Å.The bond angle distribution (BAD) for central carbon (C-C-C) and central nitrogen (C-N-C), as depicted in Figure S4b, displayed similar pattern, with both having a maximum peak around 118• (close to the 120• angle in graphite).These suggest that nitrogen substitutions within the layered nanostructure maintain the precise sp 2 layered arrangement of the carbon atoms.Consequently, the impedance of electronic conduction pathways by nitrogen atoms does not stem from their structural attributes within impure amorphous graphite.
Our prior research highlighted that electronic conductivity in amorphous graphite favors path with interconnected 6membered (hexagonal) rings [28].Building on this understanding, we now observe that the introduction of nitrogen atoms, even within a 6-membered ring arrangement, disrupts electronic conduction within the planes.This insight suggests that purposefully incorporating nitrogen into graphite derived from other carbonaceous materials could provide unique avenues for electronic conduction.The average electrical conductivity (σ) was calculated for both amorphous graphite and impure amorphous graphite with varying impurity weight percentages (see Figure 6f).The electronic conductivity of pure carbon amorphous structure measured significantly lower by a factor of 100 compared to crystalline graphite [28].The electronic conductivity of nitrogen-containing impure amorphous graphite was ≈ 20% lower than that of amorphous graphite.The incorporation of oxygen and sulfur impurities further decreased the conductivity values.In comparison to amorphous graphite, the electronic conductivity dropped by 58% and 64% in impure amorphous graphite with 5% and 10% impurity concentrations, respectively.These estimates are qualitative but indicative of the transport consequences of the impurities.

Conclusions
The thrust of this work is to offer a first-of-its-kind, pragmatic perspective into the structure, vibrational behavior, and electronic properties of impure amorphous graphite, derived from coal-like atomistic models through high-temperature transformation.Employing a novel simulation protocol with the REAXFF potential, designed to replicate the initial carbonization stages leading to single-molecule coke formation, this research unveils enduring functional groups containing noncarbon elements during the pyrolysis of bituminous coal at 1200 K. Leveraging ab initio DFT, the graphitization of models featuring coal-like elemental compositions, including carbon atoms with 5 and 10 weight-percent of nitrogen, oxygen, and sulfur atoms, is explored.The findings emphasize the overall vibrational and electronic structure of the impure amorphous graphite bear significant similarity to that of amorphous graphite.However, the presence of nitrogen, particularly in graphitic form, impedes electronic conductivity within layers.This research advances both methodologies and insights for exploring coal's alternative applications beyond energy production, offering valuable contributions to its integration into cutting-edge electronic technologies.
We include additional figures to supplement our main discussion.Figure S1 displays the initial P8 coal configuration post-density optimization at 0.4 GPa and 2.0 GPa, along with resulting coke models after STEAM implementation.These models retain 5-and 6-membered rings with residual non-carbon elements.Figure S2 outlines the STEAM process for 0.4 GPa coke, while a similar plot for 2.0 GPa coke is in the manuscript.
In Figure S3a-d, we present the electronic structure of (a) amorphous graphite and impure versions with (b) 5% and (c) 10% non-carbon elements.Figure S3d illustrates the electronic structure for a layered nanostructure with 5% nitrogen impurity.A representative model of nitrogencontaining amorphous graphite is shown in Figure S3e, and atoms localized at the state with IPR = 0.51 are depicted in Figure S3f.Radial distribution function (RDF) and bond angle distribution (BAD) for nitrogen-containing amorphous graphite are presented in Figure S4a and b, respectively.
Figure S1: The starting and final configurations from both the 0.4 GPa and 2.0 GPa simulations are displayed.Following the process of carbonization, the ultimate configuration reveals the presence of five-and six-membered rings alongside sp carbon chains.A small number of hydrogen, nitrogen, oxygen, and sulfur atoms are retained in coke.The hydrogen, carbon, nitrogen, oxygen, and sulfur atoms are colored white, gray, blue, red, and yellow, respectively.

S2
Figure S2: An overview of the simulation cycle illustrating the formation of coke under a pressure of 0.4 GPa is shown.The main manuscript elaborates on a detailed, systematic methodology that was employed to achieve these results.Notably, in the eighth cycle, we took the precautionary step of reiterating the final NVT phase.This measure was taken to ensure that no new bond cleavage events could occur during this stage of the simulation.The coke formed from this particular simulation is also presented in Figure S1.

Figure 1 :
Figure 1: (a) The 2D coal model proposed by Solomon [41], and (b) the same model converted into a 3D stable structure using ChemDraw ® .Polar hydrogen atoms are colored in white, carbon in gray, nitrogen in blue, oxygen in red, and sulfur in yellow.

Figure 2 :
Figure 2: (a) The density time evolution of various pressures.The time development for the energy and pressure convergence at 0.4 GPa is shown in (b).The P8 coal model obtained from the NPT simulation shown (a) and (b), as well as the coke obtained after the NVT/NPT iterative simulation are also shown in (c).

4 .
Conduct NVT simulation, ensuring adequate atom removal intervals to allow a reasonable probability for molecule fragment recombination (e.g., start with 1000 timesteps). 5. Execute the NPT step while closely monitoring density convergence.6. Iterate between NVT and NPT phases until no bond cleavage occurs in NVT step.7. Re-run the last NVT step from item 6 to confirm absence of new bond cleavages.8. Verify that mainly high-molecular-mass molecules remain, ideally a single molecule.9. Perform an additional NPT step to attain a fully converged coke density.10.Employ conjugate gradient relaxation to acquire an energy-minimized configuration.

Figure 3 :
Figure 3: Iterative molecule removal process at 2.0 GPa.The NVT step in cycle 6 was repeated to ensure consistency.A similar plot obtained at 0.4 GPa is shown in the supplementary material.

Figure 4 :
Figure 4: (a -d) Chronology of formation of the layered nanostructure with 5% impurity concentration included (using plane-wave DFT).The brown arrows in (b) pinpoint the putative planes.The forming structure from the model with 10% impurity elements is displayed in (e).The corresponding radial distribution function obtained for impure amorphous graphite with (f) 5% and (g) 10% impurity concentrations.

Figure 5 :
Figure 5: (a) Vibrational density of states (VDoS) for impure amorphous graphite with 5% non-carbon elements (IaG-5%) compared with that of pristine graphite (pG).The phase quotient (Q p ) describing the vibrational modes between carbon atoms (C) and between carbon and non-carbon atoms (X) is shown in (b).The bond stretching and bending (S p ) in amorphous layered structure is shown in (c).

Figure 6 :
Figure 6: Electronic density of states (EDoS) calculated for (a) amorphous graphite (only carbon atoms present), impure amorphous graphite with (b) 5% and (c) 10 % impurity concentration, and (d) 5% nitrogen impurities only.In (a) -(d), the Fermi level (E f ) is shifted to zero and indicated by green dashed lines, and only the region between E f ± 2 eV is shown.The complete EDoS is provided in FigureS3.SPC calculation[66], presented as a grayscale heatmap on selected layers containing nitrogen, is shown in (e); the darker regions represent the electronic conduction paths in the layered nanostructures.The average electronic conductivity (σ, in S/m) for all 3 models is plotted in (f).

Figure S3 :
Figure S3: The electronic density of states (EDoS) for pure amorphous graphite (containing only C atoms) is shown in (a).The EDoS for impure system is displayed for (b) 5% impurity concentration, (c) 10% impurity concentration, and (d) only 5% nitrogen impurities.The green dashed lines indicate the shifted Fermi level, set at zero for all cases (ad).Despite the introduction of impurities, the fundamental electronic structure of amorphous graphite remains largely unaltered.A the model of amorphous graphite with a 5% nitrogen impurity concentration is depicted in (e).In (f), the distribution of localized states in this model (in e) is portrayed.The relative size of the atoms visually represents the degree of localization, with non-localized atoms indicated by brown arrows.

Figure
Figure S4: (a) illustrates the radial distribution function (RDF) for the nitrogen-containing impure amorphous graphite.Notably, the nearest neighbor distances for carbon-carbon (C-C at 1.419 Å) and carbon-nitrogen (C-N at 1.404 Å) interactions are highly comparable.Bond angle distributions (BAD) for C-C-C and C-N-C are presented in (b).Given the relatively low nitrogen concentration in the model, a histogram was generated for C-N-C angles, depicting angle frequencies.Interestingly, both the central carbon-nitrogen angle and the nitrogen-carbon angle exhibit peak values around 118 • , closely resembling angles in graphite.The BAD plots were normalized to their maximum values.This normalization suggests that nitrogen substitutions within the layered nanostructure maintain the precise sp 2 layered arrangement, extending (at least) up to the nearest carbon atom.
2. Determine initial model density via NPT simulation, leveraging R o data for the coal (if available).3. Initiate STEAM implementation by selecting appropriate simulation duration for NVT and NPT phases.A recommended runtime should exceed 100 ps, ensuring sufficient time for density convergence during the NPT step.

Table 1 :
Details from the carbonization simulation at 1273 K.The chemical composition of the starting configuration is C 805 H 700 N 10 O 75 S 10 .

Table 2 :
Nearest neighbor bond length of the models compared to experimental data on bond length from other works.The unit of measurement is in Å.