On the validity of empirical potentials for simulating radiation damage in graphite: a benchmark

In this work, the ability of methods based on empirical potentials to simulate the effects of radiation damage in graphite is examined by comparing results for point defects, found using ab initio calculations based on density functional theory (DFT), with those given by two state of the art potentials: the Environment-Dependent Interatomic Potential (EDIP) and the Adaptive Intermolecular Reactive Empirical Bond Order potential (AIREBO). Formation energies for the interstitial, the vacancy and the Stone–Wales (5775) defect are all reasonably close to DFT values. Both EDIP and AIREBO can thus be suitable for the prompt defects in a cascade, for example. Both potentials suffer from arefacts. One is the pinch defect, where two α-atoms adopt a fourfold-coordinated sp3 configuration, that forms a cross-link between neighbouring graphene sheets. Another, for AIREBO only, is that its ground state vacancy structure is close to the transition state found by DFT for migration. The EDIP fails to reproduce the ground state self-interstitial structure given by DFT, but has nearly the same formation energy. Also, for both potentials, the energy barriers that control diffusion and the evolution of a damage cascade, are not well reproduced. In particular the EDIP gives a barrier to removal of the Stone–Wales defect as 0.9 eV against DFT's 4.5 eV. The suite of defect structures used is provided as supplementary information as a benchmark set for future potentials.


Introduction
The response of graphite to damaging radiation is extremely complex, owing to its anisotropic structure, and the ability of carbon atoms to adopt a great variety of metastable configurations, yet the material can remain recognisably graphitic (with modified properties) even after exposure that Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. is sufficient to displace every carbon atom in a specimen tens of times. This resilience, together with several other properties of graphite, make it attractive for applications in the nuclear industry. Indeed, the world's first artificial nuclear fission reactor used blocks of graphite for the moderator and reflector material. It is still used in present-day reactors for this purpose, and appears often in designs proposed for the future. Thus, since the modification of properties that does occur upon exposure to damaging radiation may affect the performance of operating nuclear reactors, radiation damage in graphite has remained a subject of interest since the 1940s.
The process of radiation damage begins in a fission reactor with the displacement of a carbon atom by a neutron. A chaotic cascade then follows, which eventually settles into a population of Frenkel defects and defect complexes, accompanied by the formation of dislocation loops [1]. These defects then affect the properties of the irradiated material compared with the pristine state.
Insight into the nature of defects created by radiation damage at the atomic scale has in recent times been facilitated by the discovery of graphene, and by progress in atomicresolution imaging techniques. High-resolution electron micrographs can show clearly how the atoms are arranged [2,3], and their dynamical behaviour [4,5]. Graphite, however, is a three-dimensional solid, and much less amenable to experimental investigation, and even less so for specimens which have become contaminated with radioactive matter after long-duration exposure to radiation inside nuclear reactors. Also, since the observations are made after irradiation, it is not possible to determine how the atoms reached their destinations. To do this, it is necessary to model radiation damage events using molecular dynamics. Except for the small size and short time duration events near the displacement threshold [6], the number of atoms, timescales, and lengths involved in radiation damage events is far too great to use ab initio methods to simulate these processes. Empirical potentials must be employed, instead. Their use inevitably raises questions about fidelity and transferability, owing to the vast range of possibilities for hybridization and coordination that can occur during a radiation damage event in graphite. Nevertheless, it has been seen previously that empirical potentials are able to provide good models of carbon in a wide variety of situations [7][8][9][10].
The tests used to validate empirical potentials for carbon in the past have mainly been based on statistics of structural parameters, typically using model systems that are initially prepared in a high-temperature liquid-like state, with a set density, which is subsequently quenched using molecular dynamics [7][8][9]. These are compared with equivalent ab initio simulations. This random, sampled approach covers a wide range of different bonding configurations without any obvious bias. The complementary defect-physics methodology, used in the present work, is to prepare specific models of defects and defect complexes, then compare results with different methods. Appropriate defects in this instance are the fundamental, intrinsic defects of the graphite lattice, namely self-interstitials and lattice vacancies, and stable pairings of them. This type of approach provides a more tightly controlled, targeted analysis of the relative performance of different potentials versus ab initio benchmarks, since particular structures can be examined on a case-by-case basis. One example where this has been done to a limited extent for a few selected defects is in [10]; however, this is not the main part of that work, and there is insufficient information to use it as a benchmark.

Methods
The ab initio calculations follow the same methodology as described in [11,12].
These are based on density functional theory (DFT) implemented by the AIMPRO program package [13][14][15] to calculate the total energies of supercells in which model defects have been constructed. The supercells are orthorhombic, comprising 6 × 3 × 2 eight-atom unit cells when no defect is present. Fixed lattice parameters are used for the DFT calculations, following the same protocol adopted in [12]. The dimensions of the unit cells are a = 2.445 821 Å, b = a √ 3, and c = 6.518 6297 Å, which represent the values optimized using the local density approximation with the AIMPRO program package. It is necessary here to specify these lengths with more digits than is physically significant, in order to keep the lattice vectors numerically consistent with the atomic coordinates given in the supplementary data, since even a small mismatch can make calculations fail in some instances. Results are given for both the PW92 local density approximation (LDA) by Perdew and Wang [16], and the PBE96 generalized gradient approximation (GGA) by Perdew, Burke, and Ernzerhof [17]. Core electrons are replaced by norm-conserving pseudopotentials [18].
A wide variety of potentials are available for carbonbased systems. For the purposes of the present work, we have selected two that are frequently used, and are applicable to graphitic systems: AIREBO (Adaptive Intermolecular Reactive Empirical Bond Order) [19] and EDIP (Environment-Dependent Interatomic Potential) [20,21] reformulated for carbon [7]. This choice should not be seen as any form of advocacy, nor does the non-selection of other possibilities mean we view them unfavourably; they are only meant as representative examples.
Each of the structures optimized with the AIMPRO program package using the LDA (the structures optimized with GGA are nearly the same) is reoptimized by the LAMMPS [22] and EDIP [7] packages. The LAMMPS [22] package includes the AIREBO potential as an option, while the EDIP package is specific to the EDIP. Some additional calculations were also done using the GULP package [23,24]. In all these calculations, we use lattice parameters that are appropriate to each potential. For the AIREBO potential calculations, the orthorhombic unit cell dimensions are a = 2.418 Å, b = a √ 3 Å, and c = 6.738 Å, which are the values that give the minimum total energy. In the case of the EDIP calculations a and b are the same as the DFT values stated previously, and c = 6.4 Å exactly, which is double the cutoff parameter c 0 , beyond which carbon atoms are assumed not to interact. The lattice parameters measured by x-ray crystallography of hexagonal graphite are a = 2.4612 ± 0.0001 Å, and c = 6.7079 ± 0.0007 Å [25].

Results
A total of 33 point defects and defect pairs are considered in this work. These include 25 defects that have been described previously, plus six metastable structures for bound pairs of interstitial atoms, and two for a pair of vacancies that have not been examined with DFT before. Also, the level of detail provided in earlier descriptions varies. In the present work, all defects are described in full, in order to facilitate comparison and future work. Illustrations of the defects, and their atomic coordinates (as individual plain text xyz-format files) are provided in supplementary information.

Self-interstitial atoms
Isolated, single self-interstitial atoms are predicted by DFT to exist in four distinct forms: spiro, grafted, α-split, and β-split, reviewed previously in [26]. These are illustrated in figure 1.
The two types of split-interstitial have a pair of carbon atoms sharing an α or β site in a graphite crystal, with the bond between the pair aligned in the prismatic direction [27]. They both possess D 3h symmetry. According to Lieb's theorem, substitutional defects spontaneously induce a magnetic moment in a π-electron system on a hexagonal bipartite lattice [28]. DFT shows this is true for both types of split-interstitial, and predicts M ≈ 2µ B .
Grafted interstitial atoms form a triangular structure bridging the bond between two neighbouring atoms in the same graphene sheet of the host, and have C s symmetry. A similar grafted adatom defect is also observed on the surface of graphite [29]. The length of the bond that is bridged in the host sheet is calculated to be 1.54 Å (LDA) or 1.53 Å (GGA), compared with its unperturbed length of 1.42 Å. This corresponds to the bond order being single, as in diamond, where the C-C bond length is also 1.54 Å [30]. In contrast to earlier DFT results for adatoms on graphene [31], the present calculations find that this defect in graphite has no magnetic moment.
The spiro state comprises four nearest neighbouring atoms from the host, two in each of two adjacent layers, plus the interstitial atom near their centre of mass, with six bonds linking the five atoms [27]. It has C 2 symmetry. Spinpolarized DFT confirms the expected result that this defect is not magnetic.
The results of our calculations are summarized in table 1, where it can be seen that the extent to which both potentials are able to reproduce the DFT structures and energies is in most cases at least qualitatively acceptable, and in a few cases quantitatively good. These are all challenging structures for potentials, owing to their triangular bonding arrangements, which are known from previous studies on models of disordered carbon to be underrepresented in the statistics [8]. In particular, the EDIP fails to reproduce the spiro structure, and instead finds a minimum for the cantedinterstitial structure, which has C 2h symmetry (figure 1). According to DFT calculations, this structure has an energy more than 3 eV above the spiro, and is highly unstable. The result of optimizing a canted-interstitial by DFT with the conjugate-gradient method is a spiro, while the AIREBO method finds a local minimum in energy in between the canted and the spiro, with E f = 6.61 eV. Nevertheless, in EDIP simulations the canted-interstitial can represent a partial substitute for a true spiro-interstitial structure owing to its energy being close to that given by DFT. The EDIP and AIREBO potential energy surfaces have another energy minimum in between the canted and grafted states, that is unstable for DFT. This state can be viewed as being a distorted grafted interstitial, where the interstitial atom spans a gap of 2.2 Å between the α and β atoms in the neighbouring graphene sheet, instead of there being a bond. Its formation energy, calculated with the EDIP program package, is 8.64 eV, while the AIREBO formation energy is 6.15 eV. The defect is illustrated in the supplementary material. Both potentials overestimate the formation energies of the split interstitials, and their energies are much higher than for the grafted interstitial, contrary to the results with DFT. Nevertheless, the grafted structure is reproduced well by both potentials.
Collision cascades and annealing processes can be simulated using molecular dynamics. This covers a much larger configuration space than that occupied by point defects.
A systematic way to sample this space for validation of empirical potentials is to compare the results of calculations for selected reactions using either the nudged elastic band method, or constrained optimization, with DFT benchmarks. In graphite, the process of self-interstitial migration provides several suitable examples.
A spiro-interstitial defect can reorient its C 2 symmetry axis to point along one of six 1 12 0 directions via a transition state known as a 'Y-lid' [32], which has a Y-shaped bonding configuration, illustrated in figure 1. This process exchanges one pair of the interstitial atom's β-atom bonds, while both α-atom bonds and one β-atom bond remain intact. A simple constraint, R 2 β1 −R 2 β2 = 0, where R β1 and R β2 are the distances between the interstitial and the two exchanged β-atoms, can keep the atoms in the Y-lid state configuration during structural optimization. Both the LDA and GGA yield very similar energies for this state relative to the spiro-interstitial, i.e. E a = 1.12 eV and 1.14 eV, respectively, according to our calculations. Since the spiro state is not a stable minimum for the EDIP, the reorientation process is not applicable in this instance. Spiro reorientation for the AIREBO potential is possible, at least; however, the Y-lid structure is a minimum in energy 0.07 eV lower than the spiro. Thus, its formation energy is E f = 5.32 eV.
Reorientation via a Y-lid yields no net migration, since the interstitial atom always remains attached to one pair of α-atoms. To reach a neighbouring pair of α atoms, a migrating interstitial atom must first escape from the spiro state to a neighbouring grafted state, then move to a split-interstitial on a β-site, followed by the next neighbouring grafted state, before reaching a new spiro state.
The first migration step, spiro to grafted, has a high activation energy: our results using DFT are E a = 2.00 eV for the LDA and E a = 1.87 eV for the GGA. This is very close to the results of Zhang et al who report that E a = 1.88 or 2.12 eV for spiro to grafted, using a GGA [33].
It can be seen that both types of split-interstitial state provide a route for self-diffusion in the prismatic direction. Both carbon atoms in these defects are equivalent; thus, both of them have an equal opportunity to move into a neighbouring grafted state, regardless of which of the pair was formerly in a grafted state. This also means that the activation energy for self-diffusion by this process will be the same in the prismatic direction as it is in the basal plane, since escape from the spiro state dominates the overall process. This is not entirely obvious for such an anisotropic material as graphite. Nevertheless, if the atoms are labelled isotopically, then the diffusion rate in the prismatic direction will appear to be very much less than it is in the basal plane, because for a particular atom to pass through a graphene sheet, after losing its partner in a split state, Figure 1. Interstitial atoms in graphite. From top to bottom: the spiro interstitial, β-split interstitial, grafted interstitial, Y-lid transition state, and canted interstitial. The α-split interstitial is similar to the β. The canted interstitial is predicted to be the most stable state by the EDIP, but DFT finds it to be highly unstable.
it must dwell as a host atom until another interstitial atom is able to displace it via the split state.
For the EDIP, the process of migration necessarily begins from the canted-interstitial, which is its version of the spiro state. The activation energy to reach the grafted state then is 3.02 eV, which makes it slightly over 1 eV too large. AIREBO, however, performs well: the result for spiro to grafted is E a = 2.14 eV. The next step is grafted to β-split. DFT predicts that this barrier is low: E a = 0.23 eV using the LDA and E a = 0.20 eV for the GGA, calculated with AIMPRO. The barrier from grafted to α-split is slightly higher: E a = 0.31 eV using the LDA and E a = 0.27 eV for the GGA. Neither the EDIP nor the AIREBO potential are able to reproduce these results. First, for both potentials, the split states are about 2 eV higher than the grafted, while for DFT split is slightly lower than grafted. This Table 1. Calculated properties of isolated self-interstitial atoms in graphite: E f is the formation energy, R α is the distance from the interstitial atom to the nearest α-atom or partner sharing the α site for the α-split state, R β is likewise for β-atoms, and M is the magnetic moment (DFT only).

Self interstitial pairs
In previous work it was shown that self interstitial atoms in graphite can form pairs with binding energies of over 3 eV, relative to isolated spiro interstitial atoms [11]. In all, seven structures were identified and described. The present work adds six more. A description of each defect follows, given in order of descending binding energy. The first five structures are shown in figure 2; illustrations of the others are provided in the supplementary information, together with their atomic coordinates.
Bipentagon grafted intralayer bridge. The interstitial pair forms a bond that is about 1.40 Å long, aligned along a [1 12 0] direction, and which lies parallel with the basal plane [11]. This bond divides a hexagon in the host graphene sheet into two distorted pentagonal rings of atoms, so that they share a common side along the bond. Two separated, distorted heptagonal rings are created by the insertion of each interstitial atom between the pairs of atoms on the opposite sides of the hexagonal ring in the host that is divided by the interstitial pair bond. The overall structure accommodates the two additional atoms by creating a bulge in the graphene sheet to which they are attached ( figure 2). This defect has been described as an 'inverse Stone-Wales' defect [34]; however, the nomenclature is misleading because a Stone-Wales is generated by the displacement of atoms only, without any additional atoms. 'Inverse' has a strict meaning in group theory, which is not applicable in this case. Bent twin-triangle interlayer bridge. This defect is formed from two grafted interstitial atoms occupying the common interstitial space between neighbouring graphene sheets.
The bond between the two interstitial atoms is about 1.29 Å long, and has its axis at an angle to the prismatic direction [11]. The two pairs of host atoms to which each interstitial atom is bonded lie on different, parallel (1 12 0) planes ( figure 2).
Flat twin-triangle interlayer bridge. Two α-atoms, two β-atoms, and the two interstitial atoms form a pair of coplanar triangles on a (1 12 0) plane in this defect. The interstitial atoms are separated by a distance of about 1.28 Å, with their bond tilted at an angle to the prismatic direction, forming a bridge between two neighbouring graphene sheets ( figure 2). This defect is not stable for both potentials: optimization opens the triangles resulting in different structures in each case. αβ double interlayer bridge. One of the interstitial atoms is bonded to two α-atoms in neighbouring graphene sheets, and the other interstitial atom is bonded to two β-atoms, such that all six atoms lie on the same (1 12 0) plane. Four of these atoms, two each from the host and the interstitial pair, form a square on one side of the defect. The length of the bond between the interstitial atoms is about 1.33 Å (figure 2). ββ bent interlayer bridge. In this defect the interstitial pair forms a bent structure bridging two β-atoms in adjacent sheets, with all four atoms lying on the same (1 12 0) plane [11]. The bond between the interstitial pair is about 1.21 Å long, and oriented at an angle to the prismatic direction (figure 2). Bipentagon interlayer bridge. For this defect, four αatoms and the two interstitial atoms lie on the same (1 01 0) plane, with the axis of the bond between the interstitial pair perpendicular to the prismatic direction, forming the shared side of two pentagonal rings of atoms [11]. The two interstitial atoms are separated by a distance of about 1.35 Å. αβ bent interlayer bridge. Similar to the ββ bent interlayer bridge, this defect is a bent structure bridging two adjacent sheets, except that one of the interstitial atoms is bonded to an α-atom in the host instead of a β-atom. The two host atoms and the interstitial pair lie on the same (1 12 0) plane, and the bond length for the interstitial pair is about 1.21 Å. Twisted twin-triangle interlayer bridge. This defect is formed from two grafted interstitial atoms on neighbouring graphene sheets. The length of the bond between the pair is about 1.28 Å and lies in a (1 12 0) plane, with its axis oriented at an angle to the prismatic direction. The two pairs of host atoms to which each interstitial atom is bonded lie on (1 12 0) planes at 30 • to each other, giving the defect a twisted form. This defect is not stable for both potentials: optimization yields the αβ bent interlayer bridge.
αα arch bridge. The axis of the bond between the two interstitial atoms is aligned along a [1 12 0] direction, and oriented perpendicular to the prismatic direction for this defect. The interstitial pair form a bridge beween two second-neighbour α-atoms in the same graphene sheet, passing opposite a β-atom, with no bonds to the adjacent graphene sheet. The bond length for the interstitial pair is about 1.30 Å. ββ arch bridge. This defect takes the form of a bridge between two second-neighbour β-atoms in the same graphene sheet, with no bonds to the adjacent graphene sheet. The bond between the interstitial pair is aligned along a [1 12 0] direction, and oriented perpendicular to the prismatic direction, with the two atoms separated by about 1.31 Å. Skew bipentagon interlayer bridge. The bond between the interstitial pair forms the common side of two pentagonal rings, with each interstitial bonded to two β-atoms in neighbouring graphene sheets of the host [11]. These six atoms lie in the same plane, with its normal vector oriented at an angle to the prismatic direction. The axis of the bond between the interstitial atoms is perpendicular to the prismatic direction, and points along a [1 12 0] direction. This bond is about 1.35 Å long. αβ double split pair. This defect comprises two split interstitial pairs on neighbouring α and β sites, separated by about 1.51 Å. The four atoms all lie on the same (1 12 0) plane, in a nearly square configuration.

Isolated pentagons intralayer defect.
The isolated pentagons defect is obtained by performing a 90 • Stone-Wales rotation on one of the pentagon-heptagon bonds in a bipentagon grafted intralayer bridge, leaving both pentagons and both heptagons isolated from one another by a hexagonal ring [11,35]. The distance between the two interstitial atoms is about 1.43 Å.
Binding energies for all thirteen interstitial pairs are given in table 2. The binding energy is defined here as being the energy needed to separate the pair of the complex into two isolated interstitial atoms in their lowest energy state (i.e. spiro for DFT and AIREBO and canted for EDIP). The extent to which systematic errors cancel is expected to be better for this type of calculation than it is for calculating formation energies. It is apparent that the LDA and GGA binding energies are in good agreement with each other (the standard deviation is ±0.2 eV and largest difference is 0.37 eV). However, the binding energies calculated for each defect using the two potentials deviate by an order of magnitude more than these amounts from the the DFT results. Specifically, the EDIP and AIREBO results both have standard deviations of about ±1.7 eV for the binding energies of each defect from the individual mean DFT values (i.e. half the sum of the LDA and GGA binding energy of each defect). The largest deviations are 3.25 eV (EDIP bent twin-triangle interlayer bridge) and 3.90 eV (AIREBO αβ double split pair). Both potentials severely overestimate the binding energies of the two bipentagon interlayer bridges, and both underestimate the binding energy of the bipentagon grafted intralayer bridge. This defect, which is a dislocation dipole, must have lower energy (meaning higher binding energy) than the isolated pentagons defect, where the atoms are rearranged so as to increase the separation of the two dislocations. The difference in energy between the two states given by both potentials is close to the DFT result: LDA/GGA/EDIP/AIREBO give 2. 36 Not all the structures for the interstitial pairs are reproduced well by the two potentials. They do not work well when the defects contain triangular arrangements of atoms. Typically, optimization opens one of the three bonds forming a triangle, giving a strongly distorted form of the defect, or results in another defect, but the result is not necessarily the same for both potentials. Overall, it is quite a mixed picture when looking at the defects on an individual basis, but on average over all these defects the behaviour resembles the results of DFT calculations.

Vacancies
Both EDIP and AIREBO yield values for the formation energies of isolated vacancies that are reasonably close to those given by DFT calculations. Details of the results are given in table 3, where it can be seen that the structures appear to be a close match as well, including the Jahn-Teller distortion that lowers the symmetry of the defect from D 3h to C 2v on both the α and β sites. This is noteworthy owing to the fact that the effect has a quantum-mechanical origin [36], which is not included in the potentials. However, there is a serious problem with the AIREBO potential, which predicts that the splitvacancy structure (where an atom is located midway between two unoccupied lattice sites) has an energy that is 2.32 eV lower than the normal C 2v structure.
A crucial test for the potentials is the formation energy of a Frenkel pair, since its creation is a fundamental process in any collision cascade. Experiments on electron-irradiated graphite have concluded that the energy released by the recombination of a Frenkel defect (or needed to form one) is close to 13.5 eV [37]. The LDA and GGA predict this energy to be 13.8 and 13.6 eV, respectively. The formation energy of a Frenkel pair calculated using the EDIP is 12.2 eV, which is within 10% of its measured value. The AIREBO potential presents us with two choices for which type of a vacancy to use when calculating this quantity. If the split state is disregarded as being the actual ground state, then the result is E f (V I ) = 13.6 eV; however, when the split-vacancy is used, then E f (V I ) = 11.3 eV. DFT calculations have shown that vacancy migration and coalescence in graphite (and graphene) is a much more dynamic process than previously believed [12,38,39]. Both theoretical [12,40] and experimental [41] estimates for the activation energy for migration, including our own, suggest that E a ≈ 1 eV, and can be lower in the right circumstances [39,41]. The EDIP gives E a ≈ 0.79 eV, while the AIREBO potential fails completely because it has a minimum in energy the split-vacancy structure, which is close to the true transition state. Thus, the EDIP simulates the behaviour of isolated lattice vacancies in graphite (or graphene) rather well.

Vacancy pairs
Bound pairs of vacancies fall into two classes: coplanar and cross-layer. Many coplanar divacancy structures have been identified in previous work on graphene [5,40,[42][43][44][45][46], and graphite [12], of which possibly seven could be regarded as being 'close', in the sense that they have significant binding energy, and are not more than third-neighbour separation. When close vacancies are created in neighbouring layers, they reconstruct by forming bonds that bridge the gap. These cross-layer divacancies come in four forms with different symmetries and binding energies up to about 3 eV [12,32]. It has been suggested [32] that the trapping of vacancies in this manner may explain the discrepancy between the longheld view, based on an analysis of experimental data, that the activation energy for vacancy migration appears to be quite large (E a = 3.1 ± 0.2 eV [47]), while more recent studies suggest it is much smaller. Extended structures may also grow from cross-layer divacancies by trapping additional vacancies [38].
A description of each defect follows; selected examples are shown in figure 3, and the remainder are illustrated in the supplementary information.
Haeckelite structure divacancy. This defect comprises three fivefold and three sevenfold coplanar rings of carbon atoms, arranged alternately around its principal axis of symmetry [40,48]. In graphite, the central atom on the axis is collinear with a line of β atoms in the host crystal [12]. It can be constructed from a nearest neighbour αβ divacancy by a 90 • rotation of one pair of two pairs of atoms that are equivalent in the central ring, and adjusting the bond lengths and angles so that they are compatible with graphite. Butterfly defect divacancy. Starting from a haeckelite structure divacancy, the butterfly defect is constructed by making a second 90 • rotation of the pair of atoms on the opposite side of the original eightfold ring in a nearest neighbour αβ divacancy. This produces a hexagonal ring of carbon atoms at its centre, oriented 30 • to the hexagons of the host, surrounded by four fivefold and four sevenfold rings [42,44]. In graphite, the mirror plane of this defect is the same (1 01 0) plane as the nearest neighbour αβ divacancy from which the defect is constructed. Nearest neighbour αβ divacancy. Reconstruction of this defect yields a structure with a flattened eightfold ring with two fivefold rings on each side. The two pentagonal rings are equivalent in graphene, but not in graphite, owing to the presence of the host crystal. 'Second neighbour' divacancy. This is a complex defect, intermediate in structure between the third neighbour trans divacancy and the nearest neighbour αβ divacancy, hence the choice of name, even though it is not possible to construct a true second neighbour divacancy without extensive rearrangement of the structure. The most characteristic feature of this defect is a nearly square ring of atoms, neighboured by a flattened sevenfold ring, and a pentagon [46]. Twin pentagon-heptagon divacancy.
As its name suggests, this defect comprises two fivefold and two sevenfold rings arranged as a pair. Observed in reduced graphene-oxide, this defect represents a dislocation dipole [49]. It can be constructed from a nearest neighbour αβ divacancy by a 90 • rotation of one pair of carbon atoms in its central eightfold ring [5]. In graphene the twin pentagon-heptagon structures are identical giving C 2h symmetry; in graphite they are not, and there is no symmetry. Offset trans third neighbour divacancy. There are two ways in which coplanar third neighbour divacancies can be constructed. In this case both missing atoms lie on a (1 12 0) plane. This leaves a pair of atoms forming a bridge across the centre of the stucture, which are offset slightly from the plane of the neighbouring atoms, lowering the symmetry of the defect from C s to C 1 [12]. Cis third neighbour divacancy.
In this case, two third neighbour atoms on a (1 01 0) plane are missing. The defect does not reconstruct. In fact all atoms around this defect remain very close to their ideal lattice positions [12].
Single-cross second-neighbour ββ divacancy, V ββ * 2,1 . The lowest energy cross-layer divacancy defect has a single reconstruction bond between a pair of α atoms next to two empty second-neighbour β sites, all four of which lie on the mirror plane of the defect [12,32]. Double-cross close ββ divacancy, V ββ 2,2 . When two β atoms are removed from neighbouring layers that are closest to each other, then one or two reconstruction bonds between α atoms can form that bridge the gap. This defect is the version with two cross-layer bonds [12]. The crosslayer reconstruction nearly vanishes when using the EDIP. Single-cross close ββ divacancy, V ββ 2,1 . This version of the close ββ cross-layer divacancy has only one bond between two α atoms, which lowers the symmetry of the structure from C 2h to C 2 [12,32]. The cross-layer reconstruction nearly vanishes when using the EDIP. Single-cross αβ divacancy, V αβ 2,1 . Our calculations have confirmed earlier predictions [32] that a pair of αvacancies in neighbouring sheets is not stable; however, an α-vacancy can form a bound pair with a nearest β-vacancy in a neighbouring sheet [12]. The orientation of both vacancies is the same, locked in place by a bond between the two atoms that would be unpaired if the vacancies were isolated. The vacancies and atoms forming the cross-layer bond all lie on the mirror plane of the defect. The crosslayer reconstruction is weak when using the EDIP.
Compared with the benchmark DFT calculations, AIREBO underestimates binding energies for the coplanar divacancies, and overestimates them for cross-layer divacancies. Nevertheless, the overall level of agreement is quite good: the standard deviation from the mean DFT binding energies is about ±1.3 eV. EDIP underestimates all the binding energies, to such a large extent that none of the cross-layer states are bound, so the results cannot be considered a success. This is reflected in the large standard deviation of about ±2.7 eV from the mean DFT binding energies. The interlayer reconstruction of the cross-layer divacancies is also weak or non-existent. Thus, while both potentials may work well for graphene, only the AIREBO potential seems suitable for simulating the formation of extended interlayer structures of the type described in [38]. Further details of the results are given in table 4.

Displacement defects
Displacement defects are expected to be among the defects formed in the wake of a collision cascade. These are defects where the number of atoms is conserved over their extent. Specifically, there are three such defects that are relevant here: the Stone-Wales defect and the α and β forms of intimate Frenkel defects. The creation of a Stone-Wales defect is a pericyclic reaction, i.e. one that involves a concerted rearrangement of bonds arranged in a cyclic manner [50]. These reactions are governed by the Woodward-Hoffmann rules [51]. Owing to long-ranged elastic effects, the formation energies of Stone-Wales defects are sensitive to the size of model used, and the orientation of the defects with respect to the lattice vectors of the model system [52]. Thus, it is essential to use models that have identical size and geometry when making comparisons, which is the case in the present work. In terms of both formation energy and structure, it can be seen from the results presented in table 5 that both potentials apparently perform well with respect to DFT, with EDIP being the winner by a small margin.
Although it might not be obvious, a Stone-Wales defect can also be formed from a Frenkel defect [53]. Thus, the energy of a Stone-Wales defect can be expressed as a 'binding energy' E b , with respect to an isolated lattice vacancy and a self-interstitial. The estimate for E b given by the AIREBO potential about 1 eV lower than the DFT values, but noticeably closer than the EDIP, which underestimates this energy by about 2 eV (table 5).
The AIREBO potential also performs well when we consider the activation energies for formation and removal of a Stone-Wales defect, giving a result that is fairly close to DFT, while it is found that the EDIP underestimates the height of the barrier by slightly more than AIREBO (table 5).
The existence of energetic barriers to the collapse of Frenkel defects has been demonstrated by experiments on graphite exposed to damaging radiation [37,54,55]. Calculations based on DFT support this idea [12,53]. Metastable states exist where an interstitial atom and vacancy are in intimate proximity, that may form either by them encountering each other by migration, or are formed directly during a collision cascade. In one form, an α-atom is displaced from its lattice site, while in the other form, a β-atom is displaced. These states have a formation energy that is about 3 eV lower than a widely separated Frenkel pair, according to DFT calculations [12,53]. Both potentials reproduce the DFT results quite closely for the β intimate Frenkel defect (table 5). However, only AIREBO gives a satisfactory result for the α variant, which, in contrast to its sibling, has barely any barrier to collapse, according to DFT [12]. EDIP yields a so-called 'pinch' defect, instead. This is an unphysical state, where two atoms in neighbouring layers adopt an sp 3 configuration, and are joined by a bond. The pinch defect structure is also a metastable energy minimum with AIREBO as well, similar to the EDIP version, except that its depth is much smaller.

Summary and conclusions
Owing to the approximate nature of empirical potentials, and the compromises involved in their construction, it is necessary to identify their limitations in order to justify their use. The purpose of the present work is to provide a benchmark for carbon potentials, that is applicable to simulations of radiation damage in graphite at the atomic scale. As such, it can both test the validity of potentials, and guide their construction.
This study finds that, for 33 model defects, the EDIP and AIREBO empirical potentials are mostly able to reproduce the results of calculations based on DFT, but that there are some notable failures. These are not necessarily fatal flaws, however. In some cases when using a potential, it may be possible to accept a substitute for one of the DFT states. For example, the canted interstitial can be used to represent a spiro interstitial when doing simulations with the EDIP: the structure is wrong, but the energy is close to the DFT result for the spiro. In other instances this approach sometimes works less well, such as with single vacancies for the AIREBO potential, where there is the complication that the split-vacancy is the ground state, and its energy is more than 2 eV too low.
The potentials work less well for self-interstitial and vacancy pairs, and fail to reproduce activation energies for defect transformations. This limits their usefulness for modelling processes such as annealing, but the average properties of these defects are close enough to those given by DFT to be used in some simulations that involve a population of defects. The structures and energies of displacement defects fare much better, except for the spurious, unphysical pinch defect. Other spurious local minima found in the present work are the bridge-interstitial defect (both potentials) and the Y-lid (AIREBO). Nevertheless, at least isolated vacancies and interstitials within the empirical approximations are represented with sufficient fidelity to reproduce the primary events occuring in radiation damage, after their shortcomings are taken into account. Typically, these calculations might involve molecular dynamics simulations of cascade events, which themselves employ further approximations and assumptions that may affect the outcome more than the total energy calculation.