ReaxFF Molecular Dynamics Simulation for the Graphitization of

Subscriber access provided by UNIV OF DURHAM

Dynamics

ReaxFF Molecular Dynamics Simulation for the Graphitization of Amorphous Carbon: a Parametric Study Kejiang Li, Hang Zhang, Guang-Yue Li, Jianliang Zhang, Mohammed Bouhadja, Zhengjian Liu, Adam A Skelton, and Mansoor Barati J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01296 • Publication Date (Web): 19 Apr 2018 Downloaded from http://pubs.acs.org on April 21, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

ReaxFF Molecular Dynamics Simulation for the Graphitization of Amorphous Carbon: a Parametric Study Kejiang Li

a, *

, Hang Zhang b, Guangyue Li c, Jianliang Zhang

a,f

, Mohammed

Bouhadja d, Zhengjian Liu a, Adam Skelton g, h, and Mansoor Barati e a

School of Metallurgical and Ecological Engineering, University of Science and

Technology Beijing, Beijing 100083, P.R. China. b

Modern Technology and Education Centre, North China University of Science and

Technology, Tangshan 063009, P.R. China. c

College of Chemical Engineering, North China University of Science and Technology,

Tangshan, Tangshan 063009, P.R. China. d

Laboratoire de Physique de la Matière Condensée, Université Picardie Jules Verne,

Amiens Cedex, France e

Department of Materials Science and Engineering, University of Toronto, Toronto, ON

M5S 3E4, Canada. f

School of Chemical Engineering, The University of Queensland, St Lucia,QLD 4072,

Australia. g

Department of Chemistry, University of Liverpool, Liverpool, United Kingdom

h

School of Health Sciences, University of KwaZulu-Natal, Durban, South Africa

*Corresponding author: A/Prof. Dr. Kejiang Li E-mail: [email protected] Phone: +86-010-62332550 Address: 30 Xueyuan Road, Haidian District, Beijing, 100083, P. R. China

1 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ABSTRACT A parametric study of ReaxFF for molecular dynamics simulation of graphitization of amorphous carbon was conducted. The responses to different initial amorphous carbon configurations, simulation time steps, simulated temperatures, and ReaxFF parameter sets were investigated. The results showed that time step shorter than 0.2 fs is sufficient for the ReaxFF simulation of carbon using both Chenoweth 2008 and Srinivasan 2015 parameter sets. The amorphous carbon networks produced using both parameter sets at 300 K are similar to each other, with the first peak positions of pair distribution function curves located between the graphite sp2 bond peak position and the diamond sp3 bond peak position. In the graphitization process, the graphene fragment size increases and orientation of graphene layers transforms to be parallel with each other with the increase of temperature and annealing time. This parallel graphene structure is close to the crystalline graphite. Associated with this graphitization is the presence of small voids and pores which arise because of the more efficient atomic packing relative to a disordered structure. For all initial densities, both potential parameter sets exhibit the expected behavior in which the sp2 fraction increases significantly over time. The sp2 fraction increases with increasing temperature. The differences of sp2 fraction at different temperatures are more obvious in lower density at 1.4 g/cc. When density is increased, the gap caused by different temperatures becomes small. This study indicate that both Chenoweth 2008 and Srinivasan 2015 potential sets are appropriate for molecular dynamics simulations in which the growth of graphitic structures is investigated.

Key Words: Molecular dynamics; graphitization; amorphous carbon; ReaxFF;


Page 2 of 33

Page 3 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1. Introduction The structural transformation of carbon (coke or coal char) during heat treatment determines their properties in both ambient and high temperatures. Oberlin et al.

1-2

conducted comprehensive high resolution transmission electron microscopy (HRTEM) studies of carbon structure during heat treatment, demonstrating four stages for the change of micro-texture during carbonization up to 3000 °C

1, 3

. The first stage (< 700 °C)

corresponds to the release of aromatic CH groups. The individual scattering domains [basic structural unit (BSU)] are less than 1 nm in diameter and are isometric. Their thickness is also about 1 nm (N, the number of graphitic layers, varies from 1 to 3). The second stage (700 - 1500 °C) corresponds to the release of interlayer defects between two superimposed BSUs. N increases to 8-10. The third stage (1500 – 2000 °C) corresponds to the release of in-plane defects. A considerable increase in thickness, and then the coalescence of adjacent columns can be observed due to the disappearance of misoriented single BSU. The fourth stage (> 2100 °C) corresponds to the annealing of all distortions. The layers become stiff and perfect. The heteroatoms and both the interlayer and in-plane defects have gradually disappeared, leading to the start of graphite crystal growth. To further understand the carbonization process, molecular dynamics (MD) simulation has been employed to build various carbon models and investigate the structural transformations. Between 1942 and 2010 there were over 134 proposed molecular level representations (models) of coal, typically generated for a specific use 4. In recent years, the Reactive Force Field (ReaxFF) which is a general bond order dependent potential developed by van Duin et al. 5 has been extensively adopted6-13 to investigate the chemical



Page 4 of 33

reaction and structural variations during coal pyrolysis in atomic scale. These works were concentrated on Stages 1 and 2 of the four-stage transformations described by Oberlin et al. 1

, while the structural changes in Stages 3 and 4 have not yet been extensively studied by

MD and understanding about the graphitization process is still incomplete. After the pyrolysis process of carbonaceous materials, heteroatoms have been released and the bulk composition become mainly carbon. Therefore, structural change in stage 3 and 4 are mainly transformation of carbon structural order. Through experiments, it has been concluded that carbonaceous materials may graphitize if they are heat treated further after the carbonization state (semi-coke stage) 2. In industrial blast furnace ironmaking process, graphite crystals were frequently observed in tuyere coke which was extracted from the high temperature zone (1500 – 2700 °C) of blast furnace

14-15

. The graphitization

process of carbon in high temperature has been widely known and confirmed by experiment, but the understanding about graphitization based on atomic simulation is quite insufficient due to the difficulty to simulate pure carbon. Two main reasons are the complexity of the competing hybridizations and long-range effects associated with π electrons. To simulate carbon structure, more than 10 potentials including Tersoff

16

,

17

18

,

reactive empirical bond order (REBO)

, analytic bond-order potential (ABOP)

environmental-dependent interaction potential (EDIP) bound order potential (LCBOP)

20

19

, ReaxFF5, long-range carbon

, modified embedded atom method (MEAM)

charge optimized many body (COMB)

22

21

and

have been put forward since 1988. Recently, de

Tomas et al. 23 performed a comparative study of six common carbon interatomic potentials (Tersoff, REBO-II, ReaxFF, EDIP, LCBOP-I and COMB3) which are implemented in the molecular dynamics package LAMMPS to evaluate their ability to describe amorphous


Page 5 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


carbon and the graphitization process. It was found that ReaxFF and EDIP provide the most accurate picture for the graphitization process, providing a good foundation for the further study of the graphitization process. As ReaxFF potential is already available in LAMMPS and it has the capability to characterize the bond dissociation and bond formation process, more test should be done to accurately simulate the graphitization process by ReaxFF potential. Different ReaxFF parameter sets have been developed for a variety of systems. Parameter sets developed have been shown to accurately describe bond dissociation and formation for systems including oxidation of hydrocarbons

24

, catalytic formation of

nanotubes 25, shock waves in polymers26, and many more carbon-based systems27. However, a variety of ReaxFF force-field parameter sets and simulation variables have been reported in these studies, and it is unclear what combination of parameter sets and simulation conditions should be used for the graphitization process. The objective of this study is to determine the modeling parameters necessary for accurate simulation of graphitization process using the ReaxFF. A parametric study is performed by monitoring the responses of different initial amorphous carbon configurations using various levels of initial densities, simulation time steps, simulated temperatures, and ReaxFF parameter sets. The results of the study identify the maximum usable time step length, beyond which the simulation results spuriously depend on the value chosen. The results also indicate that choice of temperature and parameter set can have a significant effect on the simulation results. This study provides the foundations for future MD simulations on the structural transformation of amorphous or pure carbon using ReaxFF potential, which will promote advances in methodology with ReaxFF application to simulate challenging problems related with carbon structure.



Page 6 of 33

2. Parametric Modelling Procedure Density which is closely related to the fractions of sp3 bonding and sp2 bonding 28 is the key factor determining the initial carbon structure. Based on available experimental data, Merrick et al.

29

proposed a mathematical model to predict the true density (corrected for

the presence of porosity and mineral matter) of coal, semi-coke and coke. The carbon density varies from 1.3 g/cm3 to 2.26 g/cm3 depending on the initial volatile matter and heating rate as well as temperature in the carbonization process. The 2.26 g/cm3 value is the limiting case of pure graphite carbon in this model. To cover a wide range of initial carbon structure, densities of 1.4, 1.8 and 2.2 g/cm3 were selected to perform MD simulations. The ReaxFF potential

5

was originally developed for hydrocarbons and has been

expanded to describe a large number of atomic species including oxygen, nitrogen, sulfur and metal 27. The development philosophy of ReaxFF differs significantly from most other potentials. Components of the potential encompass many chemical possibilities including bond order, charge transfer, van der Waals, under/over-coordination and torsion. The current functional form of the ReaxFF potential, best described in Chenoweth et al.’s

24

hydrocarbon combustion work (herein referred to as Chenoweth 2008 parameter set, abbreviated to Che08), has demonstrated significant transferability across the periodic table 27

. For pure carbon systems simulated at 3000 K, Jensen et al.

30

found that Chenoweth

2008 parameter set yielded the formation of large graphitic structures while Nielson parameter 25 set did not. In a recent study, de Tomas et al. 23 found that a recently improved ReaxFF parameter set released by Srinivasan et al.31 in 2015 (herein referred to as Srinivasan 2015 parameter set, abbreviated to Sri15 ) can describe the graphitization well.


Page 7 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Therefore, both Chenoweth 2008 and Srinivasan 2015 parameter sets were adopted in the present study. The time step required to get accurate results with ReaxFF may be quite different from those used with traditional force field depending on the stiffest portion of the force gradient experienced by the atoms during MD simulation. If the time step is too large, the total energy of the system will not be conserved, typically increasing over the course of the simulation. Jensen et al.30 confirmed that there is a critical time step after which smaller time steps will result in the same response in the system. To encompass the range of typical time steps reported in the literature for a wide range of force fields, three different simulation time steps were selected: 0.1, 0.2 and 0.5 fs. The temperature range was selected based on the range of temperatures that carbonaceous materials can be exposed to during typical graphitization process in synthesis and operating condition. For instance, coke carbon graphitization occurred in the high temperature zone (1500 – 2700 °C) of blast furnace. Three different temperatures were modeled as part of the parametric study: 1500, 2500 and 3500 K. A three dimensional (3D) periodic MD simulation box was constructed for each initial configuration with different density, parameter set, time step and temperature combination. The same general protocol which was extensively tested by de Tomas et al.19, 23 is applied to all simulations and involves four steps. (i) the creation and equilibration of the liquid at 6000 K for 5 ps. The liquid is created from a slightly randomized simple cubic lattice at the desired density which determined the box size. (ii) the quenching of the liquid from 6000 K to 300 K within 1 ps, (1.14×1015K/S) and then equilibration at 300 K for 4 ps to produce



an amorphous solid carbon. (iii) the annealing of the amorphous solid carbon to induce graphitization at desired temperature. The evolution of carbon bond information is calculated by built-in code in the ReaxFF package. An in-house program was made to count the sp2 bond content from the trajectory data and bond information file computed by ReaxFF potential. Detailed information about this calculation method is introduced in the supporting information (session S1). package

32

All simulations were performed with LAMMPS

using periodic boundary conditions. Default Nose-Hoover thermostat in

LAMMPS package is employed to control temperature in this work, and the time constant of the thermostat is set at 10 fs. The equations of motion used are those in Shinoda et al.’s paper33, which combine the hydrostatic equations in Martyna et al.’s paper34 with the strain energy proposed by Parrinello et al.35. The time integration schemes closely follow the time-reversible measure-preserving Verlet and rRESPA integrators derived by Tuckerman et al.36. However, considering that Nose-Hoover chain is non-ergodic to some degree, especially when used in a simulation where a large amount of energy is transferred between the system and the energy bath37. This will not influence structure evolution, but might pose problem in potential future studies based on energies or fluctuations. Therefore, future research is required to validate the results against different thermostats (Andersen, Langevin, Bussi, etc.).

3. Results and Discussion 3.1 Effect of Time Step The effect of time step on the system potential energy evolution at different densities using Sri15 potential and Che08 potential at 3500 K and 1500 K are shown in Fig.S2 and


Page 8 of 33

Page 9 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig.S3 in the supporting information. It shows that the potential energy curves obtained with 0.1 fs and 0.2 fs time steps match quite well. As the potential energy is significantly influenced by the carbon structure, the sp2 bonding fraction is calculated and its evolutions with different time steps at 3500 K and 1500 K are shown in Fig.S4 and Fig.S5 respectively. As expected, there should be a maximum time step below which the sp2 bonding evolution curves show close agreement. For larger time steps, the corresponding sp2 bonding evolution curves drift obviously from those close agreement. From Fig.S4 (A-C) and Fig.S5(A-C) which show the results obtained with Sri15 potential, it can be seen that only at 1500 K with 1.4 g/cc density, the sp2 bonding evolution curve with 0.5 fs time step diverges obviously from the curves obtained with 0.1 fs and 0.2fs. For the other cases, the sp2 bonding evolution curves obtained with 0.1, 0.2 and 0.5 fs intersect with each other with similar values and trends, indicating that time steps below 0.2 fs is short enough for the graphitization simulation using Sri15 potential. For Che08 potential, unexpected computation error occurred using 0.5 fs time step and only results using 0.1 and 0.2 fs time step were shown in Fig.S4 (A’-C’) and Fig.S5 (A’-C’). It can be seen that the sp2 bonding evolution curves obtained using 0.1 and 0.2 fs time steps intersect with each other when Che08 potential is used. Therefore, time step shorter than 0.2 fs is again appropriate for the simulation using Che08 potential. To ensure accuracy, results obtained using 0.1 time step were used for further analysis in the following sections. 3.1 Amorphous carbon structure Snapshots of the amorphous carbon networks (cross-sectional slice of about 1 nm thickness) produced at 1.4, 1.8 and 2.2 g/cm3 at 300 K using Sri15 potential with 0.1 fs time step are shown in Fig. 1 (a)-(c), while the pair distribution function and sp2 fraction



evolution of each system are shown in Fig. 1 (d) and (e), respectively. The amorphous carbon structures were generated via liquid quenching as described earlier. From Fig. 1 (a)1(c), graphitic layers cannot be observed, presenting the characteristics of amorphous carbon. The higher the density, the denser the carbon network. Detailed information about the local structure of amorphous carbon can be obtained by the pair distribution function (PDF) and coordination number, as show in Fig. 1 (d) and (e). The positions of the main peak (~1.46 Å) of the PDF curves were located between the positions of the peaks of graphite (1.421 Å) and diamond (1.547 Å), indicating the local structure falls between the graphite and diamond structures. It is interesting to observe a shoulder peak near the first main peak at around 1.38 Å at the 1.4 g/cc density. The second peak is located at around 2.7 Å, and the peak position moves to a shorter interatomic distance with the increase of density, indicating that the atoms are stacked closer to each other. A small shoulder peak is also observed near the second peak at around 2.1 Å. The plateau value of the coordination number curve is near 3 which is close to that of graphitic sp2 carbon. A slight increase of coordination number is observed with the increase of density. The amorphous carbon networks and their detailed information produced using Che08 potential with 0.1 fs time step are shown in Fig. 2. The overall structure of amorphous carbon (Fig. 2 (a)-(c)) produced from Che08 potential is similar to the one by Sri15 potential. However, the local structure predicted by Che08 potential is different from that of Sri15 potential. As shown in Fig. 2 (d), a clear shoulder peak (~1.2 Å) is observed near the first main peak of the PDF curves. With the increase of density, the shoulder peak intensity decreases. The positions of the first peak in PDF curves are also located between the positions of graphite sp2 diamond sp3 bonding peaks. The position shifted to higher inter-


Page 10 of 33

Page 11 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


atomic distances with the increase of density. The shoulder peak (~2.1 Å) for the second PDF peak using Che08 potential is much smaller than that using Sri15 potential. Similar to the results produced with Sri15 potential, the first plateau value of the coordination number curve is near 3, and a slight increase of coordination number was observed with the increase of density, as shown in Fig. 2 (e). The evolution of sp2 bonding content in the liquid quenching process and the annealing process at 300 K using Sri15 and Che08 potentials are shown in Fig. 3. It can be seen that there is a sharp increase of sp2 content in the liquid quenching process from 6000 K to 300 K. The sp2 content at 300 K remains relatively steady with time. This is understandable that structural transformation at 300 K is extremely small due to the low temperature. From Fig. 3, it can be observed that the final sp2 bonding content increases with the increase of density. At density of 1.4 g/cc, the sp2 bonding content (~57%) obtained with Sri15 potential is similar to that obtained with Che08. At density of 1.8 g/cc, the sp2 bonding content obtained with Sri15 is about 67%, while that from Che08 is approximately 70%. At density 2.2 g/cc, the values increase to 79% and 73% for the Sri15 and Che08 potentials, respectively. Therefore, the sp2 bonding content of amorphous carbon is sensitive to both density and the potential sets. Even though the cooling rate of the liquid in quenching process was thought to influence the final sp2 content, McCulloch et al. [38] concluded that high speed cooling rate at 1015 K/s do not influence the final configuration noticeably.



3.2 Graphitization Process 3.2.1 Validation of model size and simulation time To validate the model size and simulation time, multiple replicates were run to test the variance of predicted values in the graphitization process. The results for the case (1.8 g/cc3500 K-Sri15-0.2) are shown in Fig. 4 as an example. The evolution of potential energy and sp2 bond fraction indicate that no variation was observed for the repeat run, as shown in Fig. 4 (a). Therefore, it is not necessary to repeat the simulation process if no parameter was changed. The evolution of potential energy and carbon bonding obtained using 2000 atoms and 3000 atoms agrees well with each other, indicating that 2000 atoms is sufficiently enough to simulate the graphitization process. Long time simulation till 800 ps indicates that system are very close to equilibrium, even though absolute equilibrium is not reached (Fig. 4 (a)). This is a common characteristic for ReaxFF simulation since the system is allowed to react and new bonds and molecules will form as the evolution of time. Jensen et al.

30

studied the structural transformation of amorphous carbon using a

simulation box with only 900 atoms and the potential energy was not absolutely conserved after 1500 ps simulation. In de Tomas et al.’s study23, a quite large carbon system (32768 atoms) was used to simulate graphitization process and no absolute conservation of potential energy were either observed. However, the potential energy and carbon structure do not change much after 150 ps simulation23. As shown in Fig. 4 (a) B, the sp2 bond content keeps at about 97% after 200

ps simulation. The corresponding carbon atom configurations were shown in Fig. 4 (b). It shows that the main carbon structure is hexatomic ring with is consisted with carbon connected by sp2 bonding, no obvious difference was observed for the configuration


Page 12 of 33

Page 13 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


obtained at 200 ps, 400 ps and 800 ps. Therefore, it is acceptable to simulate the graphitization process for 200 ps using 2000 atoms. 3.2.1 Evolution of carbon bonding In order to evaluate the effectiveness of the potentials to describe the graphitization process, annealing experiments were performed with each potential at 1500, 2500 and 3500 K. Using this approach a tendency towards graphite-like structures is expected at densities below that of crystalline graphite, 2.27 g/cc. Fig. 5 shows the snapshots after the first 200 ps of annealing for each density using the Sri15 potential. As shown in Fig. 5 (a)-(c), hexagons can be clearly observed in all the structures, and assemble to form graphene fragments of varying sizes. With the increase of density, the graphene fragment size increases and orientation of the graphene layers transform to be parallel with each other (Fig. 5 (c)). This parallel graphene structure is close to the crystalline graphite. Associated with this graphitization is the presence of small voids and pores which arise because of the more efficient atomic packing relative to a disordered structure. The lower the density, the larger the voids and pores. The bonding information of those graphitized structure can be obtained from the pair distribution function as shown in Fig. 5 (d). The first peak is positioned at 1.44 Å which is very close to the graphite sp2 bond length (1.42 Å). The plateau of the CN curve is located at 3 which is also the characteristics of sp2 bond, as shown in Fig. 5 (e). The PDF and CN curves are quite stable and do not change with the increase of density, indicating that the graphitization degree is not sensitive to the density. The snapshots after the first 200 ps of annealing for each density using Che08 potential and their corresponding pair distribution functions and CN curves are shown in Fig. 6. Similar to the results produced using Sri15 potential, hexagons and graphene



fragments of varying sizes can be clearly observed. The number of parallel graphene layer produced using Che08 potential seems to be higher than that produced using Sri15 potential. The position of the first peak (1.45 Å) in PDF curves obtained with Che08 potential is also located close to the sp2 bond length, but slightly higher than that obtained with the Sri15 potential. The plateau of CN curves is also stabilized at 3, indicating the characteristics of sp2 bonding. The evolution of sp2 bonding content in the annealing process at 3500 K using Sri15 potential and Che08 potential is shown in Fig. 7. At all densities, both potential parameter sets exhibit the expected behavior in which the sp2 fraction increases significantly over time. For Sri15 potential, it seems that higher density produce higher sp2 fraction, especially in the first 70 ps. After 160 ps annealing using Sri15 potential, the structures are almost entirely sp2-bonded with more than 90% is sp2 bonding. For Che08 potential, the effect of density on the sp2 fraction is not clear as the trend lines intersect with each other with the evolution of time. It is clear that the sp2 fraction produced with Sri15 potential is higher than that of Che08 potential. It should be pointed out that simulations have not yet reached full equilibrium due to the limitation of computation resource. Since ReaxFF is a bond order-based force field and its functional form depends on having all bonds defined explicitly5, calculation of bond information requires more computation resources compared with other potentials which do not consider bond evolution. In de Tomas’s simulation on amorphous carbon with various potentials23, it was found that the computation cost (ms/atom/timestep) of ReaxFF is 184 times more expensive than Tersoff, 122 times more expensive than REBO-II, and 52 times more expensive than LCBOP-I. Therefore, it is


Page 14 of 33

Page 15 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


difficult to reach full equilibrium with ReaxFF potential without a huge amount of computation resources. The sp2 bonding content obtained using Sri15 potential in the present study agrees well with the results obtained by de Tomas et al.23, as shown in Fig. 7. Since 32768 atoms were used in de Tomas’s simulation, while only 2000 atoms were used in the present study, the good agreement of predicted values using 2000 atoms and that using 32768 atoms further indicates that a model size with 2000 atoms is sufficiently enough to simulate the graphitization process. The variation of sp2 bond content caused by different density become less obvious when simulation time exceed 120 ps, indicating that the ideal ultimate structure might be pure graphite structure if sufficient time is provided 3.3 Effect of temperature Temperature is thought to be the dominant factor influencing the graphitization process. By annealing metallurgical coke at 1400 °C– 2000 °C Xing et al.38 found that the graphitic carbon crystalline size Lc of cokes increased significantly from original ~ 20 Å 1400 °C at to ~ 180 Å at 2273 K with increasing heat treatment temperature. Similar graphitization process of coke carbon was observed in blast furnace tuyere coke. Gupta et al.39 found that the width of 002 carbon peak of bosh and the raceway coke is lower than the same peak for coke from other regions, indicating the greatest graphitization degree at the highest temperature region of blast furnace. Interlayer spacing d002 was also strongly affected by the annealing temperature. Increasing annealing temperature resulted in decreased d002 value and a denser structure38. Due to the high temperature annealing the graphitizable ability of coke carbon, pure graphite crystals were observed on the surface of tuyere coke14.



The effect of temperature on the sp2 bonding fraction at different densities using Sri15 and Che08 potentials is shown in Fig. 8. As expected, the sp2 bonding fraction increases with increasing temperature. The differences of sp2 bonding fraction at various temperatures are more obvious in lower density at 1.4 g/cc (Fig. 8 A and A’). With the increase of temperature, the gap caused by different temperatures becomes small. At densities of 1.8 g/cc and 2.2 g/cc, there are regions where the sp2 bonding content curves intersect with each other, showing negligible effect of temperature. However, with the evolution of time and the equilibration of system, higher temperature produced higher sp2 content, indicating a larger degree of graphitization. This is in agreement with Oberlin’s annealing experimental results2 and corresponds to the third and fourth stage during which the graphitic layer thickness increases and the layers become stiff and perfect, leading to the growth of graphite crystal.

4. Conclusions ReaxFF MD simulation was conducted to determine the response of different initial amorphous carbon configurations using various initial densities, simulation time steps, simulated temperatures, and ReaxFF parameter sets. The following conclusions were obtained. (1) The curves for potential energy and sp2 bonding content evolution obtained with 0.1 fs and 0.2 fs time steps overlap, indicating that time step shorter than 0.2 fs is adequate for the ReaxFF simulation of carbon using both potential sets. (2) The amorphous carbon networks produced using Sri15 and Che08 potentials at 300 K are similar to each other. The first peaks in the PDF curves are positioned between graphite sp2 and diamond sp3 bond peaks, while the first peak shifted to higher values with the


Page 16 of 33

Page 17 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


increase of density. The shoulder peak (~2.1 Å) for the second PDF peak using Che08 potential is much smaller than that using Sri15 potential. The plateau value of the coordination number curve is near 3 which is close to the characteristics of graphitic sp2 carbon. A slight increase of coordination number is observed with the increase of density. (3) In the graphitization simulation process, both Sri15 and Che08 potentials produce similar trends. With the increase of density, the graphene fragment size increases and graphene layers become parallel to each other. This parallel graphene structure is close to the crystalline graphite. Associated with this graphitization is the presence of small voids and pores which arise because of the more efficient atomic packing relative to a disordered structure. The positions of the first peak are at 1.44 Å using Sri15 potential and 1.45 Å using Che08 potential which are very close to the graphite sp2 bonding length (1.42 Å). At all densities, both potential parameter sets exhibit the expected behavior in which the sp2 fraction increases significantly over time. For Sri15 potential, it appears that higher density produces higher sp2 fraction, while the influence of density on the sp2 fraction is not clear as the trend lines intersect with each other with the evolution of time when Che08 potential is adopted. (4) The sp2 bonding fraction increases with increasing temperature. The differences of sp2 bonding fraction between various temperatures are more evident in lower density at 1.4 g/cc. With the increase of density, the gap caused by different temperatures becomes small.

Acknowledgements Computations were performed on the GPC supercomputer at the SciNet HPC Consortium. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute



Canada, the Government of Ontario, Ontario Research Fund - Research Excellence, and the University of Toronto. Kejiang Li and Jianliang Zhang acknowledge the support of the Chinese Fundamental Research Funds for the Central Universities (FRF-TP-17-086A1), the National Science Foundation of China (51774032), and the National Key Research and Development Program of China (2017YFB0304300&2017YFB0304303).

Supporting Information The method to calculate sp2 bond content and evolution of potential energy and sp2 bond content are included in the supporting information. This information is available free of charge via the Internet at http://pubs.acs.org

References (1) Rouzaud, J.; Oberlin, A., Structure, microtexture, and optical properties of anthracene and saccharose-based carbons. Carbon 1989, 27 (4), 517-529. (2) Oberlin, A., Carbonization and graphitization. Carbon 1984, 22 (6), 521-541. (3) Feng, B.; Bhatia, S. K.; Barry, J. C., Structural ordering of coal char during heat treatment and its impact on reactivity. Carbon 2002, 40 (4), 481-496. (4) Mathews, J. P.; Chaffee, A. L., The molecular representations of coal–a review. Fuel 2012, 96, 1-14. (5) Van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A., ReaxFF: a reactive force field for hydrocarbons. The Journal of Physical Chemistry A 2001, 105 (41), 9396-9409.


Page 18 of 33

Page 19 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(6) Li, X.; Mo, Z.; Liu, J.; Guo, L., Revealing chemical reactions of coal pyrolysis with GPU-enabled ReaxFF molecular dynamics and cheminformatics analysis. Molecular Simulation 2015, 41 (1-3), 13-27. (7) Liu, J.; Li, X.; Guo, L.; Zheng, M.; Han, J.; Yuan, X.; Nie, F.; Liu, X., Reaction analysis and visualization of ReaxFF molecular dynamics simulations. Journal of Molecular Graphics and Modelling 2014, 53, 13-22. (8) Zheng, M.; Li, X.; Liu, J.; Wang, Z.; Gong, X.; Guo, L.; Song, W., Pyrolysis of Liulin coal simulated by GPU-based ReaxFF MD with cheminformatics analysis. Energy & Fuels 2013, 28 (1), 522-534. (9) Zheng, M.; Li, X.; Liu, J.; Guo, L., Initial chemical reaction simulation of coal pyrolysis via ReaxFF molecular dynamics. Energy & Fuels 2013, 27 (6), 2942-2951. (10) Liang, Y.-H.; Wang, F.; Zhang, H.; Wang, J.-P.; Li, Y.-Y.; Li, G.-Y., A ReaxFF molecular dynamics study on the mechanism of organic sulfur transformation in the hydropyrolysis process of lignite. Fuel Processing Technology 2016, 147, 32-40. (11) Li, G.-Y.; Wang, F.; Wang, J.-P.; Li, Y.-Y.; Li, A.-Q.; Liang, Y.-H., ReaxFF and DFT study on the sulfur transformation mechanism during the oxidation process of lignite. Fuel 2016, 181, 238-247. (12) Li, G.-Y.; Ding, J.-X.; Zhang, H.; Hou, C.-X.; Wang, F.; Li, Y.-Y.; Liang, Y.-H., ReaxFF simulations of hydrothermal treatment of lignite and its impact on chemical structures. Fuel 2015, 154, 243-251. (13) Li, G.-Y.; Xie, Q.-A.; Zhang, H.; Guo, R.; Wang, F.; Liang, Y.-H., Pyrolysis mechanism of metal-ion-exchanged lignite: a combined reactive force field and density functional theory study. Energy & Fuels 2014, 28 (8), 5373-5381.



(14) Gornostayev, S. S.; Härkki, J. J., Graphite crystals in blast furnace coke. Carbon 2007, 45 (6), 1145-1151. (15) Li, K.; Zhang, J.; Liu, Y.; Barati, M.; Liu, Z.; Zhong, J.; Su, B.; Wei, M.; Wang, G.; Yang, T., Graphitization of Coke and Its Interaction with Slag in the Hearth of a Blast Furnace. Metallurgical and Materials Transactions B 2016, 47 (2), 811-818. (16) Tersoff, J., Empirical interatomic potential for carbon, with applications to amorphous carbon. Physical Review Letters 1988, 61 (25), 2879. (17) Brenner, D. W., Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films. Physical Review B 1990, 42 (15), 9458. (18) Pettifor, D.; Oleinik, I., Analytic bond-order potentials beyond Tersoff-Brenner. I. Theory. Physical Review B 1999, 59 (13), 8487. (19) Marks, N., Generalizing the environment-dependent interaction potential for carbon. Physical Review B 2000, 63 (3), 035401. (20) Los, J.; Fasolino, A., Intrinsic long-range bond-order potential for carbon: Performance in Monte Carlo simulations of graphitization. Physical Review B 2003, 68 (2), 024107. (21) Lee, B.-J.; Lee, J. W., A modified embedded atom method interatomic potential for carbon. Calphad 2005, 29 (1), 7-16. (22) Liang, T.; Shan, T.-R.; Cheng, Y.-T.; Devine, B. D.; Noordhoek, M.; Li, Y.; Lu, Z.; Phillpot, S. R.; Sinnott, S. B., Classical atomistic simulations of surfaces and heterogeneous interfaces with the charge-optimized many body (COMB) potentials. Materials Science and Engineering: R: Reports 2013, 74 (9), 255-279.


Page 20 of 33

Page 21 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(23) de Tomas, C.; Suarez-Martinez, I.; Marks, N. A., Graphitization of amorphous carbons: A comparative study of interatomic potentials. Carbon 2016, 109, 681-693. (24) Chenoweth, K.; Van Duin, A. C.; Goddard, W. A., ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation. The Journal of Physical Chemistry A 2008, 112 (5), 1040-1053. (25) Nielson, K. D.; van Duin, A. C.; Oxgaard, J.; Deng, W.-Q.; Goddard, W. A., Development of the ReaxFF reactive force field for describing transition metal catalyzed reactions, with application to the initial stages of the catalytic formation of carbon nanotubes. The Journal of Physical Chemistry A 2005, 109 (3), 493-499. (26) Mattsson, T. R.; Lane, J. M. D.; Cochrane, K. R.; Desjarlais, M. P.; Thompson, A. P.; Pierce, F.; Grest, G. S., First-principles and classical molecular dynamics simulation of shocked polymers. Physical Review B 2010, 81 (5), 054103. (27) Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M., The ReaxFF reactive force-field: development, applications and future directions. npj Computational Materials 2016, 2, 15011. (28) Schwan, J.; Ulrich, S.; Roth, H.; Ehrhardt, H.; Silva, S.; Robertson, J.; Samlenski, R.; Brenn, R., Tetrahedral amorphous carbon films prepared by magnetron sputtering and dc ion plating. Journal of Applied Physics 1996, 79 (3), 1416-1422. (29) Merrick, D., Mathematical models of the thermal decomposition of coal: 3. Density, porosity and contraction behaviour. Fuel 1983, 62 (5), 547-552.



(30) Jensen, B. D.; Bandyopadhyay, A.; Wise, K. E.; Odegard, G. M., Parametric study of ReaxFF simulation parameters for molecular dynamics modeling of reactive carbon gases. Journal of chemical theory and computation 2012, 8 (9), 3003-3008. (31) Srinivasan, S. G.; Van Duin, A. C.; Ganesh, P., Development of a ReaxFF potential for carbon condensed phases and its application to the thermal fragmentation of a large fullerene. The Journal of Physical Chemistry A 2015, 119 (4), 571-580. (32) Plimpton, S., Fast parallel algorithms for short-range molecular dynamics. Journal of computational physics 1995, 117 (1), 1-19. (33) Shinoda, W.; Shiga, M.; Mikami, M., Rapid estimation of elastic constants by molecular dynamics simulation under constant stress. Physical Review B 2004, 69 (13), 134103. (34) Martyna, G. J.; Tobias, D. J.; Klein, M. L., Constant pressure molecular dynamics algorithms. The Journal of Chemical Physics 1994, 101 (5), 4177-4189. (35) Parrinello, M.; Rahman, A., Polymorphic transitions in single crystals: A new molecular dynamics method. Journal of Applied physics 1981, 52 (12), 7182-7190. (36) Tuckerman, M. E.; Alejandre, J.; López-Rendón, R.; Jochim, A. L.; Martyna, G. J., A Liouville-operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal–isobaric ensemble. Journal of Physics A: Mathematical and General 2006, 39 (19), 5629. (37) Cooke, B.; Schmidler, S. C., Preserving the Boltzmann ensemble in replicaexchange molecular dynamics. The Journal of chemical physics 2008, 129 (16), 164112.


Page 22 of 33

Page 23 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(38) Xing, X.; Zhang, G.; Rogers, H.; Zulli, P.; Ostrovski, O., Effects of annealing on microstructure and microstrength of metallurgical coke. Metallurgical and Materials Transactions B 2014, 45 (1), 106-112. (39) Gupta, S.; Ye, Z.; Kanniala, R.; Kerkkonen, O.; Sahajwalla, V., Coke graphitization and degradation across the tuyere regions in a blast furnace. Fuel 2013, 113, 77-85.



Figure Captions Fig. 1 Snapshots of the amorphous carbon networks (cross-sectional slice of about 1 nm thickness) produced at 1.4 (a), 1.8 (b) and 2.2 g/cm3 (c) at 300 K using Sri15 potential with 0.1 fs time step and their corresponding pair distribution functions (d) and coordination number (e). Fig. 2 Snapshots of the amorphous carbon networks (cross-sectional slice of about 1 nm thickness) produced at 1.4 (a), 1.8 (b) and 2.2 g/cm3 (c) at 300 K with Che08 potential using 0.1 fs time step and their corresponding pair distribution functions (d) and coordination number (e). Fig. 3 The evolution of sp2 bonding content in the liquid quenching annealing processes at 300 K using Sri15 potential (a) and Che08 potential (b). Fig. 5 Snapshots at 200 ps of the graphitized carbon networks (cross-sectional slice of about 1 nm thickness) produced at 1.4 (a), 1.8 (b) and 2.2 g/cm3 (c) at 3500 K with Sri15 potential using 0.1 fs time step and their corresponding pair distribution functions (d) and coordination number (e). Fig. 6 Snapshots of the graphitized carbon networks (cross-sectional slice of about 1 nm thickness) produced at 1.4 (a), 1.8 (b) and 2.2 g/cm3 (c) at 3500 K with Che08 potential using 0.1 fs time step and their corresponding pair distribution functions (d) and coordination number (e). Fig. 7 The evolution of sp2 bonding content in the annealing process at 3500 K using Sri15 potential and Che08 potential Fig. 8 Effect of temperature on the sp2 bonding fraction at different densities using Sri15 and Che08 potentials.


Page 24 of 33

Page 25 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 1 Snapshots of the amorphous carbon networks (cross-sectional slice of about 1 nm thickness) produced at 1.4 (a), 1.8 (b) and 2.2 g/cm3 (c) at 300 K using Sri15 potential with 0.1 fs time step and their corresponding pair distribution functions (d) and coordination number (e).



Fig. 2 Snapshots of the amorphous carbon networks (cross-sectional slice of about 1 nm thickness) produced at 1.4 (a), 1.8 (b) and 2.2 g/cm3 (c) at 300 K with Che08 potential using 0.1 fs time step and their corresponding pair distribution functions (d) and coordination number (e).


Page 26 of 33

Page 27 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 3 The evolution of sp2 bonding content in the liquid quenching annealing processes at 300 K using Sri15 potential (a) and Che08 potential (b).



Fig. 4 Validation of simulation model size and simulation time for the graphitization process: (a) evolution of potential energy and sp2 bond content; (b) evolution of atomic configuration


Page 28 of 33

Page 29 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 5 Snapshots at 200 ps of the graphitized carbon networks (cross-sectional slice of about 1 nm thickness) produced at 1.4 (a), 1.8 (b) and 2.2 g/cm3 (c) at 3500 K with Sri15 potential using 0.1 fs time step and their corresponding pair distribution functions (d) and coordination number (e).



Fig. 6 Snapshots of the graphitized carbon networks (cross-sectional slice of about 1 nm thickness) produced at 1.4 (a), 1.8 (b) and 2.2 g/cm3 (c) at 3500 K with Che08 potential using 0.1 fs time step and their corresponding pair distribution functions (d) and coordination number (e).


Page 30 of 33

Page 31 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Fig. 7 The evolution of sp2 bonding content in the annealing process at 3500 K using Sri15 potential and Che08 potential and the comparison with results obtained using Sri15 by de Tomas et al.23



Fig. 8 Effect of temperature on the sp2 bonding fraction at different densities using Sri15 and Che08 potentials.


Page 32 of 33

Page 33 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents (TOC) graphic


ReaxFF Molecular Dynamics Simulation for the Graphitization of

Recommend Documents