A ReaXFF carbon potential for radiation damage studies

Although molecular dynamics simulations of energetic impacts and collision cascades in graphite have been investigated for over 25 years, recent investigations have shown a difference between the types of defects predicted by the commonly used empirical potentials compared to ab initio calculations. As a result a new ReaXFF potential has been ﬁtted which reproduces the formation energies of many of the defects predicted by the ab initio calculations and the energy pathways between different defect states, important for investigating long term defect evolution. The data sets in the ﬁtting have been added to the existing data sets used for modelling hydrocarbons and fullerenes. The elastic properties of the potential are less well modelled than the point defect structures with the elastic constants c33 being too high and c44 too low compared to experiment. Preliminary results of low energy collision cascades show many point defect structures develop that are in agreement with those predicted from the ab initio results. (cid:1) 2016 The Authors. Published


Introduction
The advent of three-body and bond order potentials of the Tersoff form in the late 1980s [1,2] meant that large scale molecular dynamics computer simulations could be carried out on covalent materials for the first time. Simulations involving sputtering of both silicon and carbon [3,4] were some of the first to be performed with these potentials. Such simulations were successful in explaining, for example, ejection patterns and sputtering yields, but for graphite, the Tersoff model did not include any interlayer bonding terms. Later, a model was produced that included such terms between atoms that were not joined via covalent bonds [5]. In the same paper displacement energy thresholds in graphite were also evaluated, but the bonding terms required a different parameterisation for graphite compared to fullerene. Later more sophisticated potentials such as AIREBO and EDIP were developed [6,7] to overcome some of these deficiencies, and such potentials have been recently used to investigate collision cascades in perfect and polycrystalline graphite [8,9]. These potential models have been mainly developed for use in materials science applications, but carbon potentials for chemical processes and fullerenes have also been developed using a reactive force field approach.
Over recent years many ab initio calculations have also been done on defect structures in graphite, and a detailed comparison of the various potentials was carried out by Latham et al. [10]. These showed some agreement but there were also some structures that the potential models predicted that did not agree with the density functional theory, (DFT) calculations.
The Reactive Force Field (ReaxFF) potential, developed by van Duin et al. [11] also depends on the bond order of atoms in a system but has many more free parameters than the Tersoff type. This potential has proven to be highly transferable, with an ability to describe covalent, ceramic and metallic materials and their interfaces, and in particular has also been used to model hydrocarbon reactions [11] and fullerenes [12]. The advantage of this type of potential is that the parameter training process can be carried out by the addition of new data sets to those that have been determined for other systems.
This paper reports on a new ReaxFF model which fits the energies of small defect structures in graphite more exactly to the DFT with the local density approximation (LDA) results. In addition since the potentials should also be able to model the long time scale evolution of defects after irradiation damage, we also include  transition pathways between different defect structures in the training process.

The defect structures and energy pathways used in the fitting process
The most stable interstitial defect in LDA [10] is the 'spiro' interstitial, where a carbon atom between two graphite planes is bonded to two atoms in each plane twisted through 60°. This and all the stable (in LDA) interstitial, di-interstitial, vacancy and divacancy defects together with their associated formation energies given in [10] are used in the fitting process. The formation energy of another important defect, the Dienes defect [13], sometimes also known as the Stone-Thrower-Wales defect after its later rediscovery, is also fitted together with the energy barrier to form this defect from a perfect graphite structure. Fig. 1 gives the 3 main stable interstitial defects predicted by LDA and Fig. 2, a comparison of the main defect formation energies calculated by LDA and reaxFF. There is generally good agreement between them.   2. A comparison between DFT (LDA) and reaxFF for the formation energies of various defects using the notation described in [10].
The terminology for the names of the defects is defined in [10] based on their geometrical appearance.
Energy pathways in LDA were calculated using the nudged elastic band method and besides the Dienes transformation, consist of pathways from the 'spiro' to grafted interstitial, the grafted to grafted, the a-split to b-split, the a-split to grafted, the b-split to grafted and the b-split to Frenkel pair.
The formation energies from the parameterisation are shown in Fig. 2 and the comparison with some of the pathways in Fig. 3. There is relatively good agreement for all cases. The terminology used in Fig. 2 is consistent with the definitions of the defect structures given in [10]. The previous reaxFF parameterisation [12] captured the main trends for the formation energies but the current parameterisation is far superior. The reaction pathways in Fig. 3 are described in the figure caption. Note that in Fig. 3a, the barrier to form the Dienes defect is very high at over 9 eV and the defect itself has an energy of over 5 eV greater than perfect lattice. The original paper by Dienes in 1952 [14] proposed this as a diffusion mechanism, however this defect is highly unlikely to form in any annealing process. The values of the elastic constants c33 and  c44 were not fitted but tested. c44 was determined as 0.1 GPa while c33 was 85 GPa. The c44 value is too low, as experiment gives a value of around 5 GPa and c33 is higher than the experimental value of 35 GPa.

Low energy displacement cascades
In order to test the potential, some low energy collision cascades were undertaken using PKA's with energies of 500 eV. The final interstitial states of 4 cascades are shown in Fig. 4. A PKA (with 500 eV of kinetic energy) was directed in a random direction in a pristine graphite lattice at 0 K. These results confirm the observation of Christie et al. [8] that it is very easy for the PKA to channel through the lattice and when it does it loses around 18 eV per layer.
Examination of the defect structures produced by the cascades, shows the existence of the main interstitials predicted by LDA. The b split interstitial state shown in Fig. 4a, initially formed in a slightly distorted state with one bond break, however, after annealing at 1000 K for 3 ps, the system relaxed to the expected state. The 'spiro' interstitial state was also observed, and is shown in Fig. 4b. An interstitial state that is two bond breaks away from the b split interstitial also formed, as can be seen in Fig. 4c. This structure remained stable at 1000 K for over 25 ps. Fig. 4d shows an interstitial state similar to a split interstitial. Here, one bond is in a different location, compared to the split interstitial case. The interstitial atom is not directly above any atom on the underlying graphene lattice (as in the b split interstitial case), but instead distorts the underlying lattice to form 4 bonds with atoms that are part of the same hexagon.

Conclusion
Preliminary results indicate that a new reaxFF potential could be used to model the formation and evolution over long time scales of point defects produced by collision cascades in graphite but further tests are required to determine under which circumstances it performs better than the more established models of Tersoff, AIREBO and EDIP.

Acknowledgments
This work was supported by United Kingdom EPSRC grant EP/M018822/1, UNIGRAF: Understanding and Improving Graphite for Nuclear Fission.
Appendix A. Force field parameterisation