Porous Characteristics of Three-Dimensional Disordered Graphene Networks

: The porous characteristics of disordered carbons are critical factors to their performance on hydrogen storage and electrochemical capacitors. Even though the porous information can be estimated indirectly by gas adsorption experiments, it is still hard to directly characterize the porous morphology considering the complex 3D connectivity. To this end, we construct full-atom disordered graphene networks (DGNs) by mimicking the chlorination process of carbide-derived carbons using annealing-MD simulations, which could model the structure of disordered carbons at the atomic scale. The porous characteristics, including pore volume, pore size distribution (PSD), and speciﬁc surface area (SSA), were then computed from the coordinates of carbon atoms. From the evolution of structural features, pores grow dramatically during the formation of polyaromatic fragments and sequent disordered framework. Then structure is further graphitized while the PSD shows little change. For the obtained DGNs, the porosity, pore size, and SSA increase with decreasing density. Furthermore, SSA tends to saturate in the low-density range. The DGNs annealed at low temperatures exhibit larger SSA than high-temperature DGNs because of the abundant free edges.


Introduction
The large variety of carbon materials with different chemical and structural forms are used in many fields, such as hydrogen storage devices [1], filtering membranes [2], catalysis [3], and supercapacitors [4]. Even though low-dimensional carbon-based materials, including graphene [5], carbon nanotubes [6], and fullerene [7] demonstrate unique chemical and physical properties because of their ordering and reduced size [8], disordered carbons are more widely used because of less cost and widespread resources. Disordered carbons, including activated carbons [9], carbide-derived carbons (CDCs) [10], flash graphene [11], pyrolytic carbons [12], etc., can be produced using different strategies, and the products exhibit different porosity and structures. Active carbons are usually made from carbonaceous source materials (wood, bamboo, coconut husk, etc.) through an activation process by hot gases or certain chemicals [13]. The CDCs have been synthesized from different metal carbides by etching noncarbon species in various temperatures [14,15]. Recently, flash graphene (FG) are obtained by flash joule heating of inexpensive carbon sources with low energy consumption [11,16]. These disordered carbon materials usually exhibit local ordering of a honeycomb in-plane structure, but they are noncrystalline on a large scale [17][18][19]. So, they can be regarded as 3D disordered graphene networks (DGNs) with varying graphitization levels and densities.
The porosity and specific surface area (SSA) are the critical factors that influence the performance of DGNs in gas adsorption applications [20,21]. Cabria and coworkers predicted that spherical pores' optimal radius lies in the range between 8.5 and 10.7 Å [22]. In experiments, argon sorption analysis is often used to calculate SSA, pore size distributions (PSD), and pore volumes [23]. Many groups reported that the pore size

Annealing-MD Simulations
Inspired by the formation process of CDCs [15], we used the annealing-MD simulations to construct DGNs at different conditions. During the experimental preparation of CDCs, metal atoms are removed from the carbide lattice by high-temperature chlorination treatment, leaving a pure carbon material. In the simulated reconstruction, we ignored the influence of metal atoms and only considered the carbon atoms in the lattice. The starting model containing 108,000 atoms was obtained by replicating an FCC lattice of 1.0 g/cm 3 . The atoms were assigned initial velocities corresponding to 300 K, and the temperature was maintained in the canonical (NVT) ensemble for 10 ps to reduce the energy. During this procedure, the density can be adjusted by changing the box size and remapping the atoms. After that, the temperature was raised to the target annealing temperature in another 10 ps, and annealing simulations were performed for a very long time of 2000 ps. The chlorination temperature used in experiments ranges from 600 to 1200 • C, and the process lasts for several hours [35]. However, the time scale of MD simulations is much smaller than that of experiments. To achieve a similar graphitization level with experiments, we should raise the simulated temperature to accelerate the annealing process. Tomas and coworkers suggested an Arrhenius framework to correlate the experimental temperatures and simulation temperatures [31]. We followed this framework and used target temperatures of 1000 K to 4000 K to mimic the experimental conditions (from 600 to 1200 • C).
After a long time of annealing simulations under NVT ensemble, DNGs formed, but the temperature and pressure were not suitable for further characterization. So, we lowered the system temperature to 300 K and changed the box size to meet zero pressure under the isothermal-isobaric (NPT) ensemble. After annealing, there were free atoms and gaseous clusters left in the system. Cluster analysis were performed, and we can easily distinguish the small clusters from those spanning the simulation box. These gaseous clusters were far smaller, containing less than ten atoms. So, these atoms were not considered part of the adsorbent. We followed the treatment in the literature and removed them.
In the annealing simulations, the carbon-carbon interactions were simulated by the EDIP potential provided by Marks [33]. This potential has been proved reliable to reproduce the graphitization behavior of carbons in a wide range of densities [34]. However, a drawback is that the long-range term was not included in the EDIP potential. To reduce the impact of this factor, we performed further relaxation under NPT ensemble for 200 ps, using the Adaptive Intermolecular Reactive Empirical Bond Order potential (AIREBO) [36], which includes an L-J term describing the attractive forces between layers [37][38][39]. After that, the system's energy minimization was performed to get the optimized atomic configuration of DGNs [40].
The MD simulations were iterated using the code of LAMMPS [41]. We obtained a series of DGNs with varying densities and annealing temperatures. As a representative, the DGN of 1.0 g/cm 3 annealed at 4000 K was selected to investigate the formation of structures and pores during the annealing process. The coordination, from which sp 2 and sp 3 atoms can be identified, was computed by the embedded command in LAMMPS. The carbon rings were characterized using homemade code based on the "shortest-path" algorithm [42]. To calculate the edge length, we first identified edge atoms. Second, cluster analysis was applied to these atoms. Third, bond lengths were computed in each cluster and were summed up as the total edge length. The atomic configurations were rendered in OVITO [43].

Porosity Characterization
POREBLAZER is a set of computational tools for the porosity characterization of porous materials. Its code has been released by Prof. Sarkisov and coworkers [44]. We have improved some algorithms in the code to make it efficient to calculate the porous properties of materials with a side length of more than twenty nanometers. The computation started with gridding the simulation cell into cubes of 0.5 Å in length. Then each cube was tested whether it was occupied by the carbon atoms with a diameter σ C = 3.4 Å. For cubes that were not occupied by carbon atoms, we used a probe of a hydrogen atom with a diameter σ p = 2.65 Å to test whether the probe can occupy the cube without overlapping with the carbon atom. We chose the L-J parameter for H-H interactions as the probe diameter because this parameter has been used to simulate hydrogen adsorption in carbon materials [45]. If the cube could be occupied by the probe, it was defined as a portion of pores. Furthermore, the total volume of these cubes can be computed, representing the total pore volume in the material. We performed percolation analysis on the pores. Pores were replicated by 3 × 3 × 3 along the periodic boundary. After that, cluster analysis was performed in the nonperiodic boundary condition. If a cluster can connect two opposite surfaces of the replicated box, the prereplicating cubes in the cluster are considered to belong to the open pores. The centers of porous cubes were imported to OVITO to construct the surface mesh, which can be further rendered as the morphology of pores.
The surface area of a DGN is the sum of accessible surface area for all carbon atoms. That is For each carbon atom, to calculate A i , a sphere was created at the center with a radius R = (σ C + σ p )/2. Then five hundred points were randomly selected on the sphere surface, and we tested whether the point was located in one of the accessible cubes. From the number of the accessible points N i , the accessible surface area of the carbon atom can be estimated as where S is the total surface area of the sphere.
To calculate the pore size distribution of the DGNs, first, the estimated range of pore size was divided into bins in 2 Å. Second, ten thousand points were randomly selected in the simulation box. Third, for each point, if it was located in one of the accessible cubes, the largest empty sphere that contains this point was detected. Then, the values of the bin corresponding to the size of this sphere and the smaller bins were added by one. This procedure resulted in the cumulative pore size function. Finally, the PSD can be calculated by numerical differentiation of this function.

Results and Discussion
We simulated the evolution of nanoporous disorders graphene networks from a precursor with FCC carbon lattice. In the beginning, the carbon atoms were set an ensemble of velocities corresponding to 300 K, followed by a thermal relaxation of 10 ps. The initial configuration of individual carbon atoms was highly unstable, and it collapsed to a network of carbon chains quickly. Then the temperature was raised to the target annealing temperature in another 10 ps and be kept constant until 2 ns. We performed percolation analysis on pores to test whether the pore is open or closed. Furthermore, surface area is calculated only for the open pores. As a representative, the DGN of 1.0 g/cm 3 annealed at 4000 K was investigated in detail to obtain insights into the formation and evolution of porous structures. According to structural evolution, three stages can be identified.
In the first stage (from 10 ps to 20 ps), as the temperature began to rise, the number of carbon rings increased dramatically. We identified the edges of the fragments composed of more than three rings to help to understand how the pores formed. At 10 ps, the start point of the annealing process, the structure was a network of carbon chains. The atomic structure in Figure 1c shows a small amount of separated small fragments of rings. Overall, carbon atoms distributed sparsely and evenly without certain gathering, preventing the probe atom from passing through. So, at this point, there were no open pores, and the accessible surface area was zero. As time goes, carbon atoms gathered to form more rings, and the size of sp 2 dominated fragments increased, resulting in the formation of local pores. At 17 ps, several local pores connected to form open pores, and the volume of open pores increased sharply during the following annealing process. Benefit from the increasing open pore volume and edges, the accessible surface area also increased during this stage. As shown in Figure 1d, at the end of this stage, the structure transformed into a network of polyaromatic fragments with abundant edges.
In the second stage (from 20 ps to 100 ps), the polyaromatic fragments merged into a disordered framework with interconnected sp 2 dominated layers. The volume of open pores continued to increase as the aggregation of small fragments allowed more pores connected and accessible. This assembling process involves the merging of edges, so the total edge length declined. The accessible surface area, which was affected by both the open pore volume and the edges, increased until 30 ps. Then, it decreased slightly with annealing time. At the end of this stage, a framework of graphene-like layers formed, which were mainly composed of hexagonal rings. There were also some nonhexagonal rings as abundant defects in the disordered framework. In the second stage (from 20 ps to 100 ps), the polyaromatic fragments merged into a disordered framework with interconnected sp 2 dominated layers. The volume of open pores continued to increase as the aggregation of small fragments allowed more pores connected and accessible. This assembling process involves the merging of edges, so the total edge length declined. The accessible surface area, which was affected by both the open pore volume and the edges, increased until 30 ps. Then, it decreased slightly with annealing time. At the end of this stage, a framework of graphene-like layers formed, which were mainly composed of hexagonal rings. There were also some nonhexagonal rings as abundant defects in the disordered framework.
In the last stage (from 100 ps to 2000 ps), the framework underwent further graphitization. A large portion of nonhexagonal rings transformed into 6-membered rings. In the final structure, the number of hexagonal rings is 0.443 per atom, near 0.5 in infinite graphene. According to the Arrhenius framework, the model simulated at 4000 K for 2 ns has the same reaction probability with the one annealed at 1200 °C for 10 4 s. The sp 2 fractions in the model annealed at 4000 K accounted for 97.1%, which is close to that of CDCs prepared at 1200 °C for hours in experiments [46], indicating the strength of the simulation approach to mimic experimental results. We can tell that the structure finally transformed into a network of disordered graphene layers. During this stage, the free edges continued to decline because of the merging of layers and healing of defects. The change of pore structure was not significant since the main framework was already established in the former stage. In the last stage (from 100 ps to 2000 ps), the framework underwent further graphitization. A large portion of nonhexagonal rings transformed into 6-membered rings. In the final structure, the number of hexagonal rings is 0.443 per atom, near 0.5 in infinite graphene. According to the Arrhenius framework, the model simulated at 4000 K for 2 ns has the same reaction probability with the one annealed at 1200 • C for 10 4 s. The sp 2 fractions in the model annealed at 4000 K accounted for 97.1%, which is close to that of CDCs prepared at 1200 • C for hours in experiments [46], indicating the strength of the simulation approach to mimic experimental results. We can tell that the structure finally transformed into a network of disordered graphene layers. During this stage, the free edges continued to decline because of the merging of layers and healing of defects. The change of pore structure was not significant since the main framework was already established in the former stage.
The surface area and pore volume do not include all the information we need about pore structure. The pore size is also one of the main factors influencing adsorption and mass diffusion in the material. In experiments, pore size can be estimated indirectly from experimental isotherms for gas adsorption. From the computational point of view, we can compute the PSD function from the coordinates of carbon atoms. The computed PSD at several annealing times are presented in Figure 2a. We found that the PSD became wider over annealing time. At 20 ps, the maximum diameter was 1.15 nm. It increased to 1.71 nm at 100 ps. The average pore size, calculated from the PSD, also increased dramatically from 20 to 100 ps. From the rendered pore pictures in Figure 2b-d, we observed the formation and growth of pores. At the beginning, the open pores formed through the connecting of initially segregated closed pores during the first stage, and the size was even and small. As the graphene-like framework developed during the second stage, the pore size became bigger, and slit pores were found in the interlayer space of stacked structures. The PSD and rendered pore structure in Figure 2e show little change in the last stage, during which the average pore size changed slightly from 0.93 nm to 0.96 nm. After annealing, there are some free atoms and gaseous triangles, involving about 0.16% of the initial atoms. These atoms were removed from the system. After energy minimization, the average pore size of the final model is 1.09 nm, and the maximum pore size is 2.19 nm. Experimentally, the PSD is calculated from adsorption isotherms. The experimental average pore size of the DGN annealed at 1200 • C is 1.28 nm [46], slightly larger than our results. A minor portion of 2.0-4.0 nm sized mesopores were found in experiments. To describe the minor portion of mesopores in simulations, a much larger system size is needed, which is beyond our computing capability. pore structure. The pore size is also one of the main factors influencing adsorption and mass diffusion in the material. In experiments, pore size can be estimated indirectly from experimental isotherms for gas adsorption. From the computational point of view, we can compute the PSD function from the coordinates of carbon atoms. The computed PSD at several annealing times are presented in Figure 2a. We found that the PSD became wider over annealing time. At 20 ps, the maximum diameter was 1.15 nm. It increased to 1.71 nm at 100 ps. The average pore size, calculated from the PSD, also increased dramatically from 20 to 100 ps. From the rendered pore pictures in Figure 2b-d, we observed the formation and growth of pores. At the beginning, the open pores formed through the connecting of initially segregated closed pores during the first stage, and the size was even and small. As the graphene-like framework developed during the second stage, the pore size became bigger, and slit pores were found in the interlayer space of stacked structures. The PSD and rendered pore structure in Figure 2e show little change in the last stage, during which the average pore size changed slightly from 0.93 nm to 0.96 nm. After annealing, there are some free atoms and gaseous triangles, involving about 0.16% of the initial atoms. These atoms were removed from the system. After energy minimization, the average pore size of the final model is 1.09 nm, and the maximum pore size is 2.19 nm. Experimentally, the PSD is calculated from adsorption isotherms. The experimental average pore size of the DGN annealed at 1200 °C is 1.28 nm [46], slightly larger than our results. A minor portion of 2.0-4.0 nm sized mesopores were found in experiments. To describe the minor portion of mesopores in simulations, a much larger system size is needed, which is beyond our computing capability. From the above results, we know that the pores in the annealed disordered graphene networks have complicated shapes because of the curving and stacking of graphene layers. Furthermore, the porous structures are determined at the very early annealing stage. During the most prolonged stage, the pore size and volume show little change. We can estimate the corresponding experimental time of each stage. The simulation at 4000 K for 100 ps can be related to 1200 • C for 500 s based on the Arrhenius framework [31]. However, the real experimental time needed for forming the framework will be longer than 500 s because it takes additional time to remove noncarbon atoms from the matrix.
Then we investigate the influence of density on the porous properties of DGNs. The initial density of carbon atoms ranging from 0.1 g/cm 3 to 2.3 g/cm 3 was controlled by changing the box size during the first 10 ps when the number of carbon atoms was 108,000. Because of the changes in box size during relaxation, the final density is different from the initial density. So, we only refer to the final density when mentioned.
The level of graphitization can be characterized by the fraction of sp 2 atoms. For infinite graphene, sp 2 atoms account for 100%. In our simulations, the models at all densities show the same high degree of graphitization when the annealing temperature was 4000 K. As shown in Table 1, the fractions of sp 2 atoms are around 98% from the lowest density 0.10 g/cm 3 to the highest density 1.64 g/cm 3 . The carbon rings in all models are dominated by hexagons, and each atom shares 0.45 hexagons on average, which is close to the limit of infinite graphene, 0.5. The difference comes from a small amount of nonhexagonal rings and edges of limited-size graphene layers. We also calculated the total edge length of all models. The specific length sizes are similar, around 20 × 10 10 m/g. In 2011, López studied the porosity structure of a carbon system based on Tersoff potential [27]. Their models annealed at high temperature have similar specific edge lengths with our results. To compare the porous properties at different densities, we calculated the ratio of porosity. As shown in Figure 3a, the total porosity decreases almost linearly with the increase of density. This tendency is easy to explain. From the definition of total porosity, it can be related to a function of the density of disordered graphene networks, ρ, and the density of dense carbon material, ρ 0 :  We also calculated the ratio of open porosity and closed porosity by percolation analysis. At the very low density of 0.1 g/cm 3 , nearly all pores are open because there is no stacking of graphene layers, enabling the connection and opening of all void space. When the density increases, stacking of graphene layers forms, and some interlayer voids become inaccessible. So, the open porosity decreases faster than total porosity. At the same time, closed porosity increases. When density is larger than 1.0 g/cm 3 , graphene layers are large enough to divide the pores into segregated spaces, and some spaces are wrapped This equation directly shows the linear relationship between the total porosity and the density of the material. By linearly fitting the data of total porosity in Figure 3a, the intercept is 99.15%, which agrees well with the equation, indicating the rationality of the method to calculate the pore volume. The fitted ρ 0 is 2.50 g/cm 3 , corresponding to the density of single-layered graphene. It should be noted that this density is based on the sphere model of atoms, which includes less volume than the membrane model with a thickness of 3.4 Å. So, the fitted ρ 0 is larger than the often used density of crystalline graphite, 2.267 g/cm 3 .
We also calculated the ratio of open porosity and closed porosity by percolation analysis. At the very low density of 0.1 g/cm 3 , nearly all pores are open because there is no stacking of graphene layers, enabling the connection and opening of all void space. When the density increases, stacking of graphene layers forms, and some interlayer voids become inaccessible. So, the open porosity decreases faster than total porosity. At the same time, closed porosity increases. When density is larger than 1.0 g/cm 3 , graphene layers are large enough to divide the pores into segregated spaces, and some spaces are wrapped and isolated from other open pores. When density reached 1.35 g/cm 3 , open porosity drops to 3.9%, even though the total porosity is 33.4%. So, the permeability of these disordered graphene networks is weak. As density continues to increase from 1.35 g/cm 3 , all types of porosity decline, and nearly all of the pores are closed.
In porous carbon materials, the SSA is highly dependent on the porosity. Figure 3b shows that the SSA of DGNs declines monotonically with the increase of density. When density is lower than 0.5 g/cm 3 , the SSA is above 2100 m 2 /g, and we find a saturation of the SSA. Given the same mass of disordered graphene networks and the same graphitization level, when the density is low enough, nearly all pores are open, so the accessible surface areas of different models are saturated to the theoretical SSA of single-layered graphene. Further reduction of density does not result in a significant increase in the SSA. At densities higher than 1.35 g/cm 3 , even though there are pores in the material, the SSA is near zero because few pores are accessible.
The rendered images of pores, including total pores and open pores, are presented in Figure 4c. For low-density models, pores are connected in the whole space, and open pores have the same shape with total pores. While at higher densities, the pores are segregated by walls of graphene layers, and only part of them are open. The intensity and distribution range of PSDs for different densities show a huge difference, as shown in Figure 4a. We obtained the maximum pore size and average pore size from the PSD value. They both decreased with increasing density. The average pore size can range from 1 nm to more than 10 nm. So, the pore size is tunable by adjusting the density of disordered graphene networks.
Another critical factor that influences the structure of DGNs is the annealing temperature. We performed a series of annealing-MD simulations by changing the target annealing temperature from 1000 to 4000 K while kept the initial densities as 1.0 g/cm 3 . The temperature settings can be related to experimental temperatures from 400 to 1200 • C. After the same relaxation and energy minimization process, we get DNGs annealed at different temperatures. We find that annealing temperature has a significant influence on the graphitization level. The sp 2 fraction increases from 74% at 1000 K to 98% at 4000 K, indicating the DGNs annealed at low-temperatures are not highly graphitized. In low-temperature models, carbon atoms form a network of small fragments with many free edges, which is different from the graphenelike layers in high-temperature models. Figure 5a shows that all DGNs have a similar total porosity, while the open porosity for 2500 K (corresponding to 800~1000 • C in experiments) is higher than other models. This result is consistent with Dash's experimental study on CDCs from Titanium carbide, where they reported that the micropore volumes at 800 and 1000 • C of are larger than other temperatures [23]. In Figure 5b, the SSA decreases dramatically with the increasing annealing temperature, and the tendency is different from the porosity, indicating that free edges have a more significant influence than the porosity in these models. In Gläsel's work on CDCs [46], the experimental SSA also decreases with increasing temperature.  Another critical factor that influences the structure of DGNs is the annealing temperature. We performed a series of annealing-MD simulations by changing the target annealing temperature from 1000 to 4000 K while kept the initial densities as 1.0 g/cm 3 . The temperature settings can be related to experimental temperatures from 400 to 1200 °C. After the same relaxation and energy minimization process, we get DNGs annealed at different temperatures. We find that annealing temperature has a significant influence on the graphitization level. The sp 2 fraction increases from 74% at 1000 K to 98% at 4000 K, indicating the DGNs annealed at low-temperatures are not highly graphitized. In low-temperature models, carbon atoms form a network of small fragments with many free edges, which is different from the graphene-like layers in high-temperature models. Figure 5a shows that all DGNs have a similar total porosity, while the open porosity for 2500 K (corresponding to 800~1000 °C in experiments) is higher than other models. This result is consistent with Dash's experimental study on CDCs from Titanium carbide, where they reported that the micropore volumes at 800 and 1000 °C of are larger than other temperatures [23]. In Figure 5b, the SSA decreases dramatically with the increasing annealing temperature, and the tendency is different from the porosity, indicating that free edges have a more significant influence than the porosity in these models. In Gläsel's work on CDCs [46], the experimental SSA also decreases with increasing temperature. The calculated PSD for DGNs annealed at different temperatures are shown in Figure  6a. When the annealing temperature rises, the distribution of pore size becomes broader. At the same time, the position of the peak, which represents the most accessible pore size, becomes larger. The average pore size is calculated and presented in each panel. It increased from 0.71 nm to 1.09 nm. The increasing pore size can be seen from the rendered images of each DGNs, as is shown in Figure 6b-e. At low temperatures, the pores are uniform, separated by walls of small fragments. While at high temperatures, space is divided by graphene layers, and big slit pores are found. The PSD of our simulations agrees with the experimental results of CDCs [15], where the PSD of ZrC-CDCs shows increasing pore size with the increase of temperature from 600 to 1200 • C. Experimental results suggested that highly graphitized carbon materials adsorbed less hydrogen than less graphitized carbons [24]. Our simulations reveal that smaller pore size results in larger SSA given the same pore volume. Therefore, one can adjust the pore size to achieve better gas adsorption performance by changing the annealing temperature during material processing. Crystals 2021, 11, x FOR PEER REVIEW 10 of 14 The calculated PSD for DGNs annealed at different temperatures are shown in Figure 6a. When the annealing temperature rises, the distribution of pore size becomes broader. At the same time, the position of the peak, which represents the most accessible pore size, becomes larger. The average pore size is calculated and presented in each panel. It increased from 0.71 nm to 1.09 nm. The increasing pore size can be seen from the rendered images of each DGNs, as is shown in Figure 6b-e. At low temperatures, the pores are uniform, separated by walls of small fragments. While at high temperatures, space is divided by graphene layers, and big slit pores are found. The PSD of our simulations agrees with the experimental results of CDCs [15], where the PSD of ZrC-CDCs shows increasing pore size with the increase of temperature from 600 to 1200 °C. Experimental results suggested that highly graphitized carbon materials adsorbed less hydrogen than less graphitized carbons [24]. Our simulations reveal that smaller pore size results in larger SSA given the same pore volume. Therefore, one can adjust the pore size to achieve better gas adsorption performance by changing the annealing temperature during material processing.

Conclusions
In this work, we have applied MD simulations to investigate the porous characteristics of DGNs using large models containing over 100,000 atoms. The DGNs of varying densities were annealed at different temperatures for 2 ns, and the annealing temperature and time can be related to experimental conditions using an Arrhenius framework. To understand how porous structures form and evolve, we chose the DGN of 1.0 g/cm 3 annealed at 4000 K as the representative and calculated the porous properties at different annealing times. In the formation of polyaromatic fragments, the volume of small open pores and SSA increased dramatically. These polyaromatic fragments then merge into larger layers and form a framework for DGNs, during which open pores connect and evolve into larger pores with a slit feature, and the SSA decrease because of the reduction of free edges. In the rest time of annealing, the porous properties show little changes.
Then, the influence of densities on the porous structures of DGNs was investigated. The ratio of porosity, pore size, and SSA increase with decreasing density. For highdensity DGNs, the pores are confined in the interlayer spaces because of the stacking of graphene layers, resulting in closed and slit pores. For low-density DGNs, there is no stacking structure, and the pores are well interconnected. We also observed a saturation phenomenon of SSA with decreasing density, which comes from a similar graphitization level and edge lengths. Finally, we investigated the porous structures of DGNs annealed at different temperatures. At low temperatures, the pores are small, and the walls are small polyaromatic fragments, which results in higher SSA than highly graphitized DGNs at high temperatures. Our results suggest that the pore size in DNGs can be adjusted by carefully choosing the density and annealing temperature, which would guide the design and procession of nanoporous carbon materials.

Conclusions
In this work, we have applied MD simulations to investigate the porous characteristics of DGNs using large models containing over 100,000 atoms. The DGNs of varying densities were annealed at different temperatures for 2 ns, and the annealing temperature and time can be related to experimental conditions using an Arrhenius framework. To understand how porous structures form and evolve, we chose the DGN of 1.0 g/cm 3 annealed at 4000 K as the representative and calculated the porous properties at different annealing times. In the formation of polyaromatic fragments, the volume of small open pores and SSA increased dramatically. These polyaromatic fragments then merge into larger layers and form a framework for DGNs, during which open pores connect and evolve into larger pores with a slit feature, and the SSA decrease because of the reduction of free edges. In the rest time of annealing, the porous properties show little changes.
Then, the influence of densities on the porous structures of DGNs was investigated. The ratio of porosity, pore size, and SSA increase with decreasing density. For high-density DGNs, the pores are confined in the interlayer spaces because of the stacking of graphene layers, resulting in closed and slit pores. For low-density DGNs, there is no stacking