# Influence of screw dislocations on stacking order in graphite

Gabriel R Francas, Anouar-Akacha Delenda, Jacob W Martin, Colin Bousige, Irene Suarez-Martinez, Nigel A Marks, Chris Ewels

► To cite this version:

Gabriel R Francas, Anouar-Akacha Delenda, Jacob W Martin, Colin Bousige, Irene Suarez-Martinez, et al.. Influence of screw dislocations on stacking order in graphite. Carbon, 2026, 247, pp.120995. 10.1016/j.carbon.2025.120995. hal-05343529

HAL Id: hal-05343529

https://hal.science/hal-05343529v1

Submitted on 3 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research paper

# Influence of screw dislocations on stacking order in graphite

Gabriel R. Francas<sup>a, ID</sup>, Anouar-Akacha Delenda<sup>b, ID</sup>, Jacob W. Martin<sup>a, ID</sup>, Colin Bousige<sup>c, ID</sup>,
Irene Suarez-Martinez<sup>a, ID</sup>, Nigel A. Marks<sup>a, ID</sup>, Chris Ewels<sup>b, ID, *</sup>

## ARTICLE INFO

### Keywords:

Screw dislocation

Graphite

AB graphite

ABC graphite

Stacking

Bernal

Rhombohedral graphite

Graphitisation

## ABSTRACT

This work examines how screw dislocations disrupt ideal stacking in graphenic materials, using a machine learning interatomic potential to model screw dislocation dipoles within different periodic cells. A novel tool is developed to assess local interlayer registration to quantify and visualise regions exhibiting stacking order. Using molecular dynamics simulations, we demonstrate that single screws exhibit greater stability in rhombohedral stacking compared to AA stacking, whereas double screws exhibit greater stability in Bernal AB stacking. The investigation reveals several mechanisms through which the ideal stacking configuration is achieved despite the presence of screw dislocations, including shearing, bond length distortion, and buckling. Furthermore, an upper threshold for the density of screw dislocations is calculated, beyond which ideal stacking cannot be realised, potentially offering an explanation for certain forms of turbostratic carbon. The findings also indicate that single screw dislocations hinder Bernal AB stacking in graphite, whereas double screw dislocations support Bernal stacking, which is significant for understanding the formation processes of graphite.

## 1. Introduction

Graphite is a crucial material for many modern industries, from electrodes for steel recycling to lubrication, composites and as a moderator material in nuclear reactors. For over thirty years it has been the primary anode material in lithium-ion batteries [1], typically comprising around 15–20% of modern car battery cells by mass [2]. As global industries increasingly shift towards renewable energy sources to tackle climate change, interest and investment into improving the reliability, speed and efficiency of energy storage is growing [3]. Graphite is likely to continue to play an important role in this area for the foreseeable future [4,5], and so improving our understanding of the nano-structure of graphite is of increasing relevance in the search for improvement in these technologies.

Since the 1960's, defects in graphite have largely been studied in the context of nuclear power moderation [6–8]. While point defects [9], basal (in-plane) dislocations [10] and stacking faults [11] have been well documented, there is less literature on extended out-of-plane line defects such as non-basal screw dislocations. The ideal structure of graphite consists of stacked graphene layers which are only weakly bonded, however, non-basal screw dislocations connect these layers through strong covalent bonding. This has been theorised to radically change thermal and electrical conductivity in graphite [12,13]

as well as influence mechanical properties and thermal expansion in pyrocarbons [14–16]. Screw dislocations may also have a strong effect on the stacking order in graphite, with single screws potentially inhibiting ideal AB stacking [13]. In addition, our recent work [17] has demonstrated the important role screws play during graphitisation. For these reasons, improving our understanding of screw dislocations in graphite is crucial for improving both synthetic graphite production and properties relevant to modern technologies involving graphite such as energy storage.

Describing dislocations in highly anisotropic layered materials has historically been a challenge, as classical dislocation theory was originally developed for isotropic metals, limiting its transferability. In recent decades, the advancement of computational methods such as density-functional theory (DFT) and molecular dynamics (MD) simulations have aided the study of nano-structures in materials science. Applying these methods to study dislocation behaviour has remained a challenge due to the need to capture a multitude of relevant effects such as the extended strain fields and the strong intra-layer and weaker inter-layer bonding, all of which influence the various stacking arrangements in graphite. The recent advancement of machine learning interatomic potentials (MLIP) however has promised the attainment of DFT level accuracy at a fraction of the compute cost [18], enabling the study of larger systems with extended dislocations.

<footnote><sup>a</sup> Dept. of Physics and Astronomy, Curtin University, Perth, Australia</footnote>

<footnote><sup>b</sup> Nantes Université, CNRS, Institut des Matériaux de Nantes Jean Rouxel, IMN, Nantes, F-44000, France</footnote>

<footnote><sup>c</sup> Université Claude Bernard Lyon 1, CNRS, LMI UMR 5615, Villeurbanne, F-69100, Lyon, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: chris.ewels@cnrs-imn.fr (C. Ewels).</footnote>
This paper begins with a background on non-basal screw dislocations in graphite. Following benchmarking of various interatomic potentials, the atomic cluster expansion for carbon (ACE) MLIP [19, 20] is chosen to investigate these screws during relaxation under various constraints. A novel algorithm is presented to quantify local stacking registration quality and type. Using this, several mechanisms are demonstrated by which preferential stacking arrangements are achieved despite the presence of screws. We also calculate an upper limit on screw dislocation density above which ideal stacking cannot be restored, suggesting that graphite dominated by single screws may be a candidate for the structure of turbostratic carbon.

## 2. Non-basal screw dislocations in graphite

### 2.1. Theory

Dislocations are linear defects found in all crystalline materials, and are defined by a Burgers vector and line vector which may be at arbitrary angles to one another. Dislocations are labelled screw, edge or mixed when these vectors are parallel, perpendicular or somewhere between these respectively. Due to the anisotropy of graphite where in-plane interactions are much stronger than out-of-plane, dislocations are further categorised by the direction of their Burgers vector [21]: those directed within the basal plane are called basal dislocations, and those out of the plane non-basal (also called prismatic where the Burgers vector lies in the c-direction). Non-basal screws appear like spiral staircases which wind through the graphene layers, connecting them into a single sheet. Depending on the size of the Burgers vector, a single wind around the screw core may involve moving through one (single-screws), two (double-screws) or many tens or thousands of layers in the case of super-screws.

In a perfect dislocation the Burgers vector is equal to a lattice vector of the crystal, leading to restrictions on the types of screw dislocations possible in different graphite stackings. Fig. 1 illustrates the structure of perfect screw dislocations within the three high-symmetry graphite stacking types. In thermodynamically unstable AA (simple hexagonal) graphite, the prismatic lattice parameter *c* is equal to one interlayer spacing and is directed normal to the basal planes, allowing perfect single screws as shown in Fig. 1(a). In ABC (rhombohedral) graphite, the prismatic lattice parameter *c* is not normal to the basal plane but rather at an angle of ~67° to incorporate the interlayer stacking. This creates perfect single screws which are on a slant relative to the basal planes as shown in Fig. 1(b). In AB (Bernal) graphite, the prismatic lattice parameter is normal to the basal planes and equal to two interlayer spacings. This necessitates a perfect double screw which appears like a double helix connecting alternate planes which wind up two layers at a time as shown in Fig. 1(c). The screws shown in Fig. 1 all have full cores, although the possibility of hollow screws where the core atoms are removed have also been explored elsewhere [22].

### 2.2. Experiment

Experimental evidence for non-basal screw dislocations appears in numerous graphitic materials. Super-screws with Burgers vector between 450 and 2000 Å were some of the earliest to be observed in natural graphite at densities of $10^6-10^8$ m⁻² using oxidation techniques combined with optical and electron microscopy [25–29]. Since then, atomic force microscopy (AFM) has been used to study screws with Burgers vectors ranging from 4 µm down to two interlayer spacings (6.7 Å) known as double screws, an example of which is shown in Fig. 2(b) [24]. The current work explores these, and screws with Burgers vector equal to one interlayer spacing (3.35 Å) known as single or elementary screws. These were first observed by Hennig [27] in pyrolytic graphite at a density of $5 \times 10^{12}$ m⁻² using a gold decoration technique in 1965 as shown in Fig. 2(a). Subsequent work using the same technique observed a dominance of single screws in pyrolytic

<figure><img src="image_3.png" /><figcaption>Fig. 1. Perfect non-basal screw dislocation models according to elastic theory in (a) AA graphite, (b) ABC graphite and (c) AB graphite. Perfect screw dislocations have a Burgers vector equal to the *c* lattice vector as shown by the arrows. The core of the screws are emphasised with balls representing the atoms, while only the bonds are shown for the surrounding structure.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. (a) Etched section of pyrolytic graphite showing screw spirals of small Burgers vector, reproduced from [23]. (b) Atomic force microscopy image of a double screw dislocation spiral in natural graphite, reproduced from [24]. (c) High resolution transmission electron micrograph of polyvinyl chloride during graphitising heat treatment at 2400 °C from previous work by the authors [17]. Interdigitated zig-zag fringes are single non-basal screw dislocations highlighted by boxes. An example is shown in the magnified inset.</figcaption></figure>

graphite, although at lower concentration of $10^7-10^9$ m⁻² [30,31]. More recently, it has been suggested that by viewing the graphenic planes of anthracite [32,33] and coking coal [34] along the basal plane using high-resolution transmission electron microscopy (HRTEM), screws can be identified as interdigitated zig-zag fringes which connect misaligned layers. Applying this reasoning, there may in fact be further examples in previous transmission electron microscopy (TEM) images from the literature which were not at the time identified as screws [35–37]. We recently confirmed these features to be indicative of screws by gathering HRTEM images, like that shown in Fig. 2(c), of polyvinyl chloride during graphitising heat treatment [17]. Screw dislocations
are also an important mechanism for crystal growth in layered materials, as recently reviewed by Zhao and Jin [38]. In the context of graphite, spiral growth along the *c*-axis is a suggested mechanism for the development of nodular graphite in carbon alloys such as steel and cast iron [39–42]. Fig. 2(b) shows an example of a fundamental double screw growth spiral in natural graphite formed within calcite [24]. The potential for controlling this growth has recently sparked interest in single and double screws for the assembly of novel spintronics devices involving twisted few-layer graphene (FLG) [43–45]. The promise for these devices lie not only in the exotic properties of ‘magic-angle’ graphene [46] but also potentially in the strong magnetic properties of the screws themselves [22,47–49].

## 2.3. Modelling

There have been relatively few computational studies of non-basal screws in bulk graphite. We previously calculated the energetics and optimised structures of the three perfect dislocations shown in Fig. 1 using DFT [13,50,51]. We suggested that upon consideration of elastic energy alone the most stable structure should be in AA graphite, however, this is complicated by the stacking fault energy of AA graphite being much higher than AB or ABC, as highlighted in Appendix A. Other works have used DFT to explore highly screw-dislocated graphite as a metamaterial with interesting mechanical [52] and magnetic [22] properties. More recently the mechanics and energetics of non-basal screws in alternating twisted multilayer graphene (ATMG) was analysed using elastic membrane theory and MD [53]. Non-basal screws have also appeared in MD-based works which develop atomistic models of graphitic carbons. Leyssale et al. [54,55] developed an image guided atomic reconstruction (IGAR) algorithm to generate models of pyrocarbon based on HRTEM images. They found that non-basal screws appeared as the structural feature which resolves misaligned graphenic fringes. A similar observation was made by Chartier et al. [56] who observed non-basal screws at the grain boundaries between crystallites following annealing of their model of nano-crystalline graphite. Using MD simulations, we found non-basal screws arise in our self-assembled model of pre-graphitic carbon [57], and studied how they form loops which are removed during simulated annealing at graphitising temperatures [58]. Aside from graphitic carbons, non-basal screws have also been modelled in other layered carbon materials such as multi-walled carbon nanotubes (MWCNTs) [59,60] and carbon onions [61].

# 3. Methodology

## 3.1. Relaxation and annealing of screwed graphite

Initial structures are created by inserting a perfect non-basal screw dislocation dipole into a hexagonal unit cell of AA or AB graphite using elastic theory with the program Atomsk [62]. The screws have full cores with the minimal radius of 1.42 Å as indicated by the red markings on Fig. 3(a). We examine cells with standard graphite lattice parameters of $a = b = 2.46$ Å and an interlayer spacing of 3.355 Å. Various layer numbers (1 to 8) and sizes (18 to 28,800 atoms/layer) are explored, however the relative placement of the screws is maintained at the $\alpha$ (1/3, 1/3) and $\beta$ (2/3, 2/3) positions in the unit cell as shown by the red marks in Fig. 3(a). We use periodic boundary conditions in all directions, which makes the line vectors of the screws effectively infinite as suggested in Fig. 3(c) and creates an ordered hexagonal dipole array in the periodic images as shown in Fig. 3(d). This placement fixes the shortest distance between neighbouring screws along the [1010] armchair direction. These serve as the initial structures which we relax and in some cases anneal using MD.

As discussed in Section 2.1, there is an important interplay between screw dislocations and the different stacking types available in graphite. To ensure the validity of the following investigation, we perform benchmarking of several MD potentials against our own DFT

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Initial system structure example of AA graphite containing a perfect non-basal screw dipole. (a) view along *c*-axis with screw positions marked in red. (b) orthographic view. (c) side-on view showing relative spacing of the screws, with line vectors marked by red-arrows. Translucent cells above and below are shown to indicate periodic boundary conditions. (d) view along *c*-axis including some nearest periodic images to show the regular hexagonal array of screw dipoles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

calculations and broader literature values. Based on these results, which are presented in detail in the Appendices, we opt for the ACE potential for carbon [19,20] as available in LAMMPS [63], which crucially captures the energy differences between AA, AB and ABC graphite. For the relaxations, the method outlined for the potentials in Appendix B is used, except with the damped dynamics fire minimisation style [64,65] being used for the atom position relaxation in the first stage of the loop, as we find this is more efficient than the cg method at finding a reliable minimum. Low temperature annealing is performed prior to relaxation for the 8-layer systems where we found a simple relaxation was not reliable in finding the fully relaxed structures. A ramped temperature profile from 1 K up to and down from 300 K is used in an NPT ensemble with target pressures set to 0 GPa. The systems are taken up to temperature over 30 ps and back down over the same time. To accelerate the re-structuring of the largest systems over 90,000 atoms a higher temperature of 1000 K is required with ramp times of 100 ps. This is still well below the threshold of ~2000 K for screw glide. A time-step of 0.2 fs is used for both the annealing and relaxation where the fire style requires a time-step. All simulation figures are rendered using Ovito [66].

## 3.2. Stacking registration measure algorithm

We develop a new algorithm to quantify the degree of registration and type of stacking within our graphitic systems. A geometric approach is used which shares some similarities with one previously developed by Hod and others [67–70] to assess the stacking registry in hexagonal boron nitride, FLG and transition metal dichalcogenides. This previous method however is only valid for a pair of flat, parallel layers. Our approach differs in that it uses a mesh representation, allowing for the assessment of local registration for sheets which may be curved. This also allows consideration of the interactions between three, rather than two layers that is important where ABC stacking is of relevance. This leads to a bond-wise definition we call registration measure. Although what follows is a geometric definition, the registration measure metric aligns well with energy surface maps between various stacking configurations as shown in the Supplementary Materials.

A mesh is generated from our atomic systems by first using an efficient implementation [71] of the ring finding algorithm described in Ref. [72] to find rings of up to eight sides. Following the method outlined in our previous works [58,73], the mesh is constructed by placing a vertex at the site of each 3-fold coordinated atom ($V_{atom}$) and
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Schematic of registration measure procedure as explained in the text. (a) Mesh composed of atom ($V_{atom}$ shown in grey) and ring-centre ($V_{ring}$ shown in blue) vertices. Mean unit normal vector $\hat{n}$ is calculated for each three-fold coordinated atom vertex. (b) Nearest vertices are identified along $\pm 3.4\hat{n}$, labelled type $\alpha$ for atom vertices, or $\beta$ for ring-centre vertices. (c) Interlayer distances $d_{002}$ are calculated as well as the offset from ideal positions for each vertex found up and down. (d) The four possible ideal stacking arrangements in graphite for three layers, assessed for one layer relative to the layers above and below. Registration measure is assigned on a bond-wise basis due to the need to consider the $\alpha$ or $\beta$ status of each pair of atoms in the up and down directions to determine the stacking type. A registration measure less than 0.5 is deemed turbostratic while 1 denotes perfect stacking. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

at the centre of each ring ($V_{ring}$). This mesh representation is necessary because in a 3D crystalline graphite system, each atom is situated above and below either another atom or the centre of a hexagonal ring. Thus, calculating the registration relies on the knowledge of both atom and ring centre positions, and entails measuring how close each atom is to overlapping with these ideal positions when projected into the layers above and below.

To measure this we first calculate a unit normal vector $\hat{n}$ for each $V_{atom}$ by finding the mean of the cross products of the vectors between itself and its atomic neighbours, as in Fig. 4(a). We then take this unit normal vector $\hat{n}$ and multiply by a sensible guess of $3.4 \AA$ for the interlayer distance and find the nearest vertex. This is shown in Fig. 4(b) where in the upwards direction we find a $V_{ring}$ (small blue vertex circled in black) while repeating the process downwards by multiplying $\hat{n}$ by $-3.4 \AA$, we find a $V_{atom}$ (grey vertex circled in black). We use the notation of $\alpha$ and $\beta$ to denote whether the nearest vertex found in the neighbouring layer is of type $V_{atom}$ or $V_{ring}$ respectively. Next, the vector connecting the central atom vertex with the closest vertex found in the neighbouring layer, shown as a dashed purple arrow in Fig. 4(c), is projected onto $\hat{n}$ to find the true local interlayer distance $d_{002}$. Finally, we use these two vectors to find the offset distance of the projected atom from its ideal position in both the layer above $o_{up}$ and below $o_{down}$ and sum these for a total offset value for each atom $o_{tot}$.

In order to evaluate the type of stacking which best describes each region of the system, we consider the $\alpha$ and $\beta$ categorisations made for each pair of neighbouring atoms, leading to the final bond-wise definition. Considering three layers, there are four potential ideal stacking configurations where each atom is aligned with either another atom or a ring centre both above and below, as shown in Fig. 4(d). In AAA stacked graphite, each atom in a pair is sandwiched by atoms and thus assigned type $\alpha$ in the up and down direction while for any pair of neighbouring atoms in ABA, one will be type $\alpha$ in both directions, and the other type $\beta$ (sandwiched by ring centres). We call the ABB stacking configuration mixed as we have essentially AA stacking in one direction and AB stacking in the other. ABC stacking is meanwhile characterised by pairs of atoms which are type $\alpha$ in one direction while $\beta$ in the other with each atom in the pair being the inverse of the other. Due to the necessity of defining stacking type by atomic pairs, we calculate the ultimate registration measure for each pair of nearest-neighbour atoms as the mean of the total offset $o_{tot}$ for each atom, normalised within the range [0,1]. The closer the value is to 1 for each bond, the more ideally stacked the atoms are, while for registration measures less than

~0.5 we label the poor registration as turbostratic. Finally we sort each pair into one of the four relevant stacking types based on their $\alpha$ and $\beta$ values and each bond can then be coloured accordingly as shown by the scale bar in Fig. 4(d). The stacking registration code can be accessed at doi.org/10.5281/zenodo.17338926.

# 4. Results

## 4.1. Perfect prismatic screws

We begin with the idealised cases of perfect prismatic single and double screws outlined in Section 2.1 as a benchmark to compare to prior works and demonstrate the new registration measure algorithm. Starting with single screws in AA graphite, we allow a 6-layer system of 21,168 atoms to relax with almost complete freedom. Both the atom coordinates and cell dimensions are allowed to vary including the angle of the c-axis. Fig. 5(a,b) presents the initial system prior to relaxation, with bonds coloured by registration measure. As expected, the entire system is coloured yellow, identifying it as AA stacked which is confirmed by the magnified inset showing all layers stacked on top of one another. The only grey regions indicating poor stacking, are around each screw core. One of these AA screws has been extracted and magnified in the figure up to a radius of two rings where the curvature limits the degree of registration.

Following relaxation, there is a significant shift as the c-axis tilts, allowing the majority of the system to change to an ABC registration as shown by the red colouration in Fig. 5(c,d). This change in stacking order is confirmed by the characteristic patterning shown in the magnified red inset, and results in an energy decrease of 34 meV/atom down to -9.162 eV/atom. Once again, there is a region of poor registration coloured grey around each screw, however here it also extends along the ramp region between the screws. The poor registration is confirmed by viewing a slice of the system in a direction normal to the ramp which is shown in the grey inset. The poor stacking here could be because the steepness of the ramp has increased significantly relative to the surrounding flat layers. This is evidenced by the structure of the screw cores, one of which has been extracted and magnified in the figure. It has become significantly slanted with its Burgers vector shifting in line with the tilt of the c-axis. In order to achieve ABC stacking, this tilt occurs along the armchair direction which is normal to the b-axis, as shown by the two black arrows in Fig. 5(d). Thus by viewing the system along the b-axis in Fig. 5(c) we see the true tilt of the screws, which is
<figure><img src="image_7.png" /><figcaption>Fig. 5. (a) Initial 6-layer, 21,168 atom structure of AA graphite containing a perfect non-basal screw dipole according to elastic theory viewed along the *b*-axis. Screw structure extracted in box. (b) Same as (a) but viewed along *c*-axis, with magnified inset confirming AA stacking as indicated by yellow colouring of registration measure. Stacking around screw cores is shown to be disrupted by the grey colouring. (c,d) Same system as (a,b) following relaxation with full freedom where the *c*-axis shifts, allowing the stacking to largely become ABC. This is reflected by the extracted screw which is now tilted, with Burgers vector aligning with the *c*-axis. The stacking is confirmed by the red colouring of the registration measure and the magnified insets in (d), which also shows a region of poor grey stacking between the screw cores. (e,f) Perfect double screw dipole in 6-layer 21,168 atom AB graphite following relaxation where AB stacking is maintained and the screw cores are reconstructed to include two *sp*³ atoms per bilayer. All bonds are coloured by calculated stacking registration measure and associated scale bar as defined in Fig. 4. Cohesive energies of each structure are also given for comparison. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

at ~67° to the basal plane. This demonstrates that perfect single screws are more stable in ABC graphite than AA graphite.

A perfect screw in AB stacked graphite is a double screw. To model this, a 6-layer AB graphite system of 21,168 atoms is constructed in which a dipole of double screws is inserted as predicted by elastic theory, as in Fig. 1(c). Fig. 5(e,f) shows the system following relaxation where AB stacking is maintained throughout the majority of the system as confirmed by the blue colouring and characteristic layering pattern in the inset. As with the previous systems, AB stacking is not achieved in the neighbourhood of the screws, and is also somewhat compromised between the cores, although not to the extent of the ABC case. The structure is just 1 meV/atom less stable than the relaxed single screw structure in Fig. 5(c,d). The extracted core structure changes significantly from the original predicted by elastic theory to adopt a pair of four-fold coordinated atoms per bilayer, in agreement with our previous DFT calculations [13].

## 4.2. 8-Layer orthogonal system

The most idealised cases above show that for two aligned single screws, if given the freedom, there is a preference to tilt into an ABC stacking. This is not, however, what is typically observed experimentally for example in Fig. 2(c) where the fringe patterns show vertical screws. This could be due in part to imaging difficulties as a slanted screw would generally not be contained to a single focal plane. It is also likely that in a real system there will be constraints such as grain boundaries or misaligned screws that keep the screws more vertical. As a first step to modelling a more tightly constrained system, we impose the condition that the *c*-axis must remain orthogonal during the relaxation. Applying this restriction first to an 8-layer system of 28,224 atoms we found that it was necessary to perform a low temperature (300 K) anneal to consistently find the minimum energy structure (see methodology Section 3.1 for details). Fig. 6(a,b) shows the structure that results from these constraints where once again the system shifts from AA stacking to a more stable stacking order. This time however, due to the constraint on the *c*-axis and number of layers, the system cannot become wholly ABC and rather has two types of layers - the

blue AB layers labelled 1 and 5, and the other red layers which are ABC. This layering pattern is explained by the structure of the screws which again become tilted to facilitate the ABC stacking, but which in this case turn back on themselves. As shown by the orange arrows in Fig. 6(b), the AB layers are those where the screws pivot direction. Again, there is a region of poor stacking in between the screw cores, however here it is much larger and asymmetrical, while also showing some small region of mixed stacking order ABB (coloured green) near one of the screw cores. The asymmetry is reflected in the opposite AB layer, caused by the symmetry broken by the screws to preserve the periodic boundary conditions.

## 4.3. Bilayer orthogonal system

In order to test a limiting case where ABC stacking is discouraged, we next take a similar initial system of 7056 atoms comprising a 2-layer AA graphite with a screw dipole. The *c*-axis is again constrained to remain orthogonal. Following the relaxation procedure outlined in Section 3.1, we see this system shifts to a largely AB stacking structure as shown in Fig. 6(c,d). Note that for consistency in presentation and to better visualise the structure we have replicated the system four times in the *c*-direction to match the 8-layer system. Because the system is only composed of two layers with a fixed *c*-axis direction, the system cannot achieve ABC stacking. Once again the region of poorest registration is between the screw cores along the ramp, this time in a symmetrical fashion. We also see more clearly here that the registration is subpar between each screw and its two other nearest screw periodic images. The region of best stacking is as far away from the screws as possible at the corners of the cell, suggesting a radial dependence in registration development.

To examine how the structure changed from AA to AB stacking during the relaxation, we plot the atomic displacements during relaxation in the *a*–*b* plane for both layers in Fig. 7(a). We see an overall shearing between the atoms in each layer having primarily moved parallel to the armchair direction between the screws, in the opposite direction in each layer. The magnitude of this shift is largest at the corners of the cell where the material is able to achieve the best
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** (a) Orthographic view of 8-layer orthorhombic system following low-temperature annealing and relaxation showing the development of a mixture of AB (blue) and ABC (red) stacking. (b) Same structure as (a) viewed from the side to show the structure of the screws which are tilted through the ABC regions. (c) 2-layer orthorhombic system, showing the development of AB stacking, with a region of poorly stacked material around the screw cores. (d) Same structure as (c) viewed from the side to show the vertical screw structure. (e,f) 1-layer orthorhombic system showing a poorly stacked, turbostratic structure which nonetheless buckles to escape the initial AA stacking, forming small patches of ABC stacking in red. The 2- and 1-layer systems are repeated up to eight layers to better visualise the structures. All bonds are coloured by calculated stacking registration measure and associated scale bar as defined in Fig. 4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

registration, while smallest between the screw cores. Around the screw cores there is evidence of a slight twisting as the arrows are angled into the core. The minimal shearing close to the screws is due to the periodic boundaries which constrain the screws to remain vertical. Fig. 7(b) shows that the most significant bond distortions occur at the regions around the screws. This is expected and reflective of the higher strain closer to the cores. The bonds around each screw core alternately expand and contract, both within and between layers. The extent of the most significant distortion aligns quite well with the extent of the poorly stacked grey region in Fig. 6(c), suggesting that the bond distortions achieve a basal dislocation loop of sorts which demarcates this transition. These mechanisms also dominate the 8-layer orthogonal system (Supplementary Materials).

## 4.4. Monolayer orthogonal system

Having observed shearing as a dominant mechanism for improving stacking registration in these systems, we increase the constraints to probe what other mechanisms may arise where shearing is not possible. This is achieved by reducing our system to a single layer of graphene, connected to itself across the periodic boundaries by the two screws. Maintaining a perpendicular c-axis prevents the layers from shearing and also prevents in-plane bond distortions from impacting stacking order. Nonetheless, following relaxation the result is again a clear move away from the initial AA stacking, but this time less successfully. As shown in Fig. 6(e,f) the system is unable to develop the long-range stacking seen in the previous examples, but still largely moves off AA stacking and becomes turbostratic, with some small pockets of ABC ordering as shown by the red patches. This is achieved through buckling, as the rippled layers now have local normals which do not point straight up at the atom above but rather at an angle, allowing it to mostly escape AA stacking. The exception is some AA wrinkles coloured in yellow, which wind through the structure where the slope of the buckles change direction. We can examine the curvature of this buckling by calculating the angular defect throughout the structure as we described in a previous work [73]. This measure quantifies the degree and type of local Gaussian curvature. Fig. 7(c) shows that the most buckled regions are those lines of AA stacking which are negatively curved ridges, as well as the screw cores themselves.

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** (a) Map of atomic displacements in the basal-plane during relaxation of the 2-layer orthorhombic system. For clarity, only 2% of atom displacements are shown and displacements are scaled by 20x. (b) 2-layer orthorhombic system coloured by bond length distortion. The colour bar is centred on the ideal bond length for AB graphite of 1.4225 Å. In both (a) and (b) the dashed pink line indicates a discontinuity in the sheet where the ramp sits between the two screws. (c) 1-layer orthorhombic system coloured by angular defect which gives a measure of the Gaussian curvature. The buckled ridges and screw cores are coloured red indicating negative curvature. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** 8-layer orthorhombic system across a variety of screw separations for the two main layer types in the supercell. The equivalent dislocation density is also shown. Stacking quality reduces with the screw separation, until at spacings of 3 nm and below the structure becomes buckled and turbostratic, similar to the 1-layer orthorhombic system. All bonds are coloured by calculated stacking registration measure and associated scale bar as defined in Fig. 4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 4.5. Registration changes with screw density

We next investigate the effect of screw dislocation density on the stacking order of the relaxed structures. In the 8- and 2-layer systems we observe a radius around the screw cores, within which it is not possible to develop quality stacking. This suggests that as the dislocation density rises (or the screw separation shrinks), it should become harder to achieve good registration between the layers. This is first tested with the 8-layer orthorhombic system across a range of dislocation densities, achieved by varying the initial cell sizes while maintaining the relative placement of the screws. By having the same number of screws in a smaller or larger cell the overall screw density is adjusted, although we could have equally maintained cell size and varied screw density directly (see Supplementary Materials). Five of these structures are shown in Fig. 8 following low-temperature annealing and relaxation. There is a clear trend where stacking quality is reduced as the screw separation decreases or as the dislocation density is increased. While the systems with 5.5 or 7.2 nm between the screws have developed significant regions of ABC and AB stacking, at 4.3 nm this is only just starting to occur. Below this, the systems are mostly grey indicating a lack of registration, i.e. turbostratic carbon throughout. There is a notable step-change in system structure between 3.4 and 3.0 nm as the layers in the 3.0 nm system are buckled, forming a structure very similar to the 1-layer case presented above. Meanwhile, the 3.4 nm has managed to flatten out but does not have sufficient distance between the screws to develop a proper stacking order. Similar trends are seen for the 2-layer system (see Supplementary Materials). In the buckled 3.0 nm system, we observe a mean interlayer distance calculated through the registration measure procedure of 3.48 Å, which reduces to 3.43 Å for the screw separation of 7.2 nm. Comparing these to the ACE values for AB and ABC graphite in Table A.1, the buckled turbostratic system shows an interlayer distance increase of 2.3% from AB. The larger system meanwhile has a spacing slightly less than that for ABC graphite which follows from it containing two layers of AB graphite.

The observed screw density trend combined with the region of poor stacking around the screw cores leads us to investigate this quantitatively using the registration measure. Fig. 9 presents plots of the registration measure for both an ABA and ABC layer of the 8-layer system. A radial measurement is taken along the armchair direction which passes through the screws for half the diagonal of the cell (the yellow arrow in the inset). Four cell sizes are presented where the screws are separated by varying distances showing that for all separations and

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Radial plots of registration measure for various screw separations in the 8-layer orthorhombic system as a measure of the distance from a single screw core. Results are shown for layers 5 and 7 as labelled in Fig. 6(b). The direction measured is shown by the yellow arrow in the inset. Regions and insets are shaded based on the type of stacking present as defined by the registration measure and associated scale bar in Fig. 4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

both the ABA and ABC layers, there is an approximately constant radius of ~0.8 nm where the registration measure is <0.5, indicating that ideal stacking is not achieved. This region is coloured in grey to reflect this. Although the threshold of 0.5 for poor stacking is somewhat arbitrary, it appears to align in this case with a sharp jump in registration measure. For the ABA layer, beyond this core on the left side of the screw we see a progressive increase in the registration measure into ABA stacking which is achieved more completely and rapidly for cells with larger screw separations. On the other side within the ramp region between the screw cores there is a slight development of mixed ABB stacking, hence the green shading. This is also the region which first develops registration measure >0.5 in Fig. 8, but we see in this plot that regardless of screw separation the registration measure remains low in this region. In contrast, both sides of the screw in the ABC layer develop rather high registration measure, which again increases with screw separation. Radial registration measure plots for the other previously
explored systems are available in the Supplementary Materials and show similar trends in turbostratic core size and increasing stacking order with screw separation.

## 4.6. Stability & Energetics

To confirm this notion of a dislocation density or screw separation threshold to achieve stacking, we investigate the energetic stability of the previously discussed systems over a range of screw separations. The most stable system across a large portion of the range is the 6-layer rhombohedral system (i.e. perfect ABC screws seen in Fig. 5(c)). Taking these energies as the baseline, Fig. 10(a) presents the energy per atom difference between this system and each of the other systems as a function of screw separation. By examining these differences, two notable separation thresholds can be identified which divide the range of screw separations into three regions illustrated by grey, red and blue wash.

The first region is found with screw separations of less than 3–4 nm, where all the orthorhombic single screw systems have nearly identical energies as they form similarly buckled structures like the 3 nm system in Fig. 8. This region is coloured grey to indicate that, due to the high screw density, no proper stacking registration can be achieved and almost all systems resort to a buckled structure with low registration measure. The only exception is the 6-layer rhombohedral, which although it is the most stable and less buckled, still fails to achieve significant ideal ABC stacking (see Supplementary Materials). The least stable system in this region is the double screw structure in AB graphite for which the data points diverge beyond the limit of the plot at the smallest screw separations due to the high strain in the cores.

The second region is marked by a divergence of the 2- and 8-layer systems away from the buckled structure of the least stable 1-layer system, and a rapid improvement in the stability of the double screw system. This is confirmation of the structural difference noted between the 3 and 3.4 nm systems in Fig. 8; in this region almost all systems are now able to flatten and achieve more ideal stacking. This region is coloured red because the perfect ABC screws are the most stable and develop more significant ABC stacking. Despite this the 1-layer system remains buckled and turbostratic while the 2-layer orthorhombic and 6-layer double screw systems continue to develop ABA stacking throughout this region. The final region begins for dislocation spacings above ~10 nm and is marked by the ABA stacked double screw system becoming more stable than the ABC stacked single screws. The implications of these different stacking regions is discussed in the following section.

## 5. Discussion

The exploration of screw density and core size presented here has relevance to stacking in graphite in general. The structure and explanation of turbostratic carbon has remained an open question within carbon science. It is commonly understood to be graphitic carbon where the layers are at random orientations to one another, thus not forming any periodic stacking sequence. The interlayer spacing in turbostratic carbon was historically determined to be 3.44 Å [74,75] although more recent studies have established that a broader range of metastable spacings around this value are possible [76,77]. In this work, below the threshold screw separation of 3–4 nm all systems converge into similar buckled structures with poor registration. We thus label this grey region in Fig. 10(a) turbostratic with the most stable structures here represented by the wavy grey diagram with many nearby screws in Fig. 10(b). This region corresponds to a screw dislocation density of between 4–9 ×10¹⁶ m⁻² which must be considered a limiting case as our model is rather idealised. For a polycrystalline system for example containing a disordered array of screws or non-hexagonal rings, this screw density threshold would likely be lower. These systems also show a 2.3% increase in interlayer spacing relative to AB graphite that is

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** (a) Energy/atom difference between the various relaxed screw dipole structures as compared with the 6-layer rhombohedral structure with perfect ABC screws. Positive energy difference indicates a less stable structure. Three regions of different stacking preferability are coloured accordingly as discussed in the text. (b) Diagrams showing the structure of the most stable configurations in each region as coloured in (a). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

very similar to the 2.5% increase observed experimentally between AB (3.355 Å [78]) graphite and turbostratic carbon (~3.44 Å). This suggests that some turbostratic carbons may be explained by an abundance of non-basal screws, for example pyrolytic graphite which can often be turbostratic [79,80]. This is supported by the early experimental observations of single non-basal screw dislocations in pyrolytic graphite [27,30,31] as well as more recently in the IGAR models of pyrocarbons [54,55] which contain varying densities of screws depending on their heat treatment.

These findings also provides insight into the role of screws during graphitising heat treatment. Our previous works [17,57] have shown that single screws form during the secondary carbonisation or pregraphitic stages. Within the model developed by Oberlin, Monthioux and others [35,81,82], small stacks of misoriented graphenic fragments known as basic structural units (BSUs) are connected into wavy columns during secondary carbonisation. We believe this is achieved initially through a high-density network of elementary screws. In many of the historic TEM images of this stage (such as Figure 2(a) of Ref. [83] or Figure 7(a)-1300 °C of Ref. [35]), the misoriented BSUs are initially separated by <2 nm, forming the wavy buckled texture we see in the grey region for our simulations of similar screw separation (see Supplementary Materials). Following further heat treatment the distance between zig-zag fringes is reduced as some screws are removed and wavy columns transform into distorted layers as in Figure 7(1)-1600 °C of Ref. [35]. This matches what is seen in our simulations, as once screw separation is beyond ~3 nm, the layers can begin to flatten as in Fig. 2(c). We have previously modelled [58] this screw removal process which can occur once they are glissile above ~2000 °C and marks the beginning of the final stages of graphitisation. We have also shown however that double screws do support AB stacking, becoming more
stable than single screws beyond a separation of ~10 nm. This raises questions around whether some single screws, rather than annihilate, may accumulate into double screws and how such a transition may occur. Unlike in the models presented here, in a physical sample even once the screws have been removed and the layers are flat, they require further heating to remove non-hexagonal rings and allow the final stacking order to develop. There remains work to be done to understand why some screws persist during the later stages of graphitisation or how they may accumulate into super-screws and inhibit crystal growth due to interactions with vacancies and interstitials as first suggested by Roscoe and Baker in 1969 [30].

The results found here are enabled by the MLIP ACE [20] which allows us to approach DFT-level accuracy at MD scale. Previous computational studies of non-basal screws in graphite with DFT have relied on systems containing the minimum number of layers to describe each dislocation, i.e. one layer for single screws, two layers for double screws [13,22,50]. While this was an adequate starting place for exploring core structure, the current work has shown that using the minimum number of layers will not necessarily capture the intermediate to long-range structural changes. The registration measure has clarified various stacking outcomes depending on layer number, as well as a radial dependence on stacking which shows that large cells are required to better model stacking behaviour around screw cores. We previously used DFT to calculate a core radius defined by where elastic theory fails of between 1.6–1.9 Å [13] in small hydrogenated graphene clusters containing single screw dislocations. We suggest the larger radius of ~8 Å calculated here is a secondary stacking core radius to the traditional core radius. This radius is still much smaller than the typical core widths found for basal dislocations in graphite [84,85], but may be related to the 20 Å core defined by in-plane strain recently calculated for a screw dipole in ATMG [53]. In fact, the stacking core radius found here seems to align with the region of most significant bond distortion as shown in Fig. 7(b), suggesting the elastic strain field as the major driver for size of this radius. Previous DFT work by ourselves also suggested that single screws could encourage AA stacking locally [51], however we have seen here that this is not the case. Instead due to the high stacking fault energy of AA stacking, the layers in the systems will shear where possible to achieve more ideal AB or ABC stacking. We have also shown that above a threshold dislocation density, the only available mechanism to improve the stacking is for the layers to become buckled. These mechanisms and the energetics of the screw cores remains an area to explore in further work.

Although we have attempted to differentiate between more and less constrained dislocation networks by varying the dislocation density, layer number and c-axis freedom, the screws themselves are all nonetheless ordered in a hexagonal array. In particular, all the systems presented here contain screws separated along the armchair direction. Some testing was done for screws separated along the zig-zag direction where we did not find a significant change in behaviour (see Supplementary Materials). To model a more realistic system in future the screws could be distributed randomly and with Burgers vectors in various directions. In addition, by using periodic boundary conditions, the screws here are effectively infinite, thereby modelling very long screws in bulk graphite. However, previous works [32] including our own self-assembly simulations [58] suggest that non-basal screws often form loops with non-basal edge dislocations. All these effects could further constrain the way stacking is resolved in screw-dislocated graphite. In addition, although this work has focused on single and double screws as a starting point, super-screws which are prevalent in natural graphite would be a natural extension.

## 6. Conclusions

The impact of non-basal screw dislocations on stacking in graphite is examined under a range of constraints. Although it was previously suggested that single screws are most stable in AA stacked graphite,

this does not outweigh the large stacking fault energy. Instead when possible, single screws shear into ABC stacking while double screws are stable in AB graphite. Below a screw separation threshold of 3–4 nm ideal stacking is not possible and all systems resort to a buckled, turbostratic structure. This has wider implications, suggesting screw-dislocated graphite as a possible explanation for turbostratic carbons as well as helping to explain the wavy column phase seen during graphitisation where screw density must be reduced before layers can flatten. These findings are made possible by the new MLIP carbon ACE, which describes the stacking energy differences between AA, AB and ABC graphite allowing for MD-scale simulations approaching DFT level accuracy. The analysis is aided by the development of a new registration measure metric to visualise and quantify the stacking order within these systems. The techniques and approach used here could in future be applied to testing the impact of screw loops, disordered dislocations arrays and super-screws.

## CRediT authorship contribution statement

**Gabriel R. Francas:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Data curation. **Anouar-Akacha Delenda:** Writing – review & editing, Methodology, Investigation, Data curation. **Jacob W. Martin:** Writing – review & editing, Validation, Supervision, Project administration, Funding acquisition, Formal analysis, Data curation. **Colin Bousige:** Writing – review & editing, Validation, Supervision. **Irene Suarez-Martinez:** Writing – review & editing. **Nigel A. Marks:** Writing – review & editing, Validation, Supervision, Project administration, Funding acquisition, Conceptualization. **Chris Ewels:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Funding acquisition, Formal analysis, Data curation, Conceptualization.

## Declaration of competing interest

All authors declare there have no competing interests to declare for this work.

## Acknowledgements

G.R.F. acknowledges support through an Australian Government Research Training Program Scholarship. J.W.M. acknowledges the support of the Forrest Research Foundation, Australia. N.A.M. acknowledges funding from the Australian Research Council, Australia through DP190101438. A.A.D, C.B. and C.E. acknowledge support from the French National Research Agency, France grant (ANR-23-CE09-0005-03). Computational resources are provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia, Australia. This work was partially supported by the “PHC FASIC” programme (project number: 51462SL), funded by the French Ministry for Europe and Foreign Affairs, France, the French Ministry for Higher Education and Research, France and the Australian academic and institutional partners, Australia.

## Appendix A. Stacking in graphite in the literature

There is an important interplay between screw dislocations and the different stacking types available in graphite. This means accurately modelling screw dislocations requires an approach that captures the various stacking behaviours in graphite.

To ground our assessment of different modelling techniques, we begin by summarising the experimental literature on the different high-symmetry stacking modes of graphite being AB (Bernal or 2H), ABC (rhombohedral or 3R) and AA (simple hexagonal), as presented in Table A.1. The most stable form of graphite is AB stacked while ABC is metastable. ABC is found in quantities ranging up to 30% [102,103] in natural graphite samples or following milling, and mostly transforms
<caption>**Table A.1**</caption>

Values from experimental and computational literature and the current work for the relative stacking energies (meV/atom) and c-axis interlayer spacings (Å) for different high-symmetry bulk graphite stacking orders. Positive energies denote stackings less stable than AB graphite.

<table><thead><tr><th rowspan="2">Method</th><th>AB</th><th colspan="2">ABC</th><th colspan="2">AA</th></tr><tr><th>d<sub>AB</sub><br>(Å)</th><th>ΔE<sub>AB→ABC</sub><br>(meV/atom)</th><th>d<sub>ABC</sub><br>(Å)</th><th>ΔE<sub>AB→AA</sub><br>(meV/atom)</th><th>d<sub>AA</sub><br>(Å)</th></tr></thead><tbody><tr><td><strong>Experimental</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>XRD [78]</td><td>3.355</td><td></td><td></td><td></td><td></td></tr><tr><td>HRTEM [86]</td><td></td><td></td><td></td><td></td><td>3.55</td></tr><tr><td>Calorimetry [87]</td><td></td><td>6.1</td><td></td><td></td><td></td></tr><tr><td>AFM [88]<sup>a</sup></td><td></td><td>0.15</td><td></td><td></td><td></td></tr><tr><td><strong>Benchmark calculations</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ACFDT-RPA [89]</td><td>3.34</td><td></td><td></td><td>10</td><td>3.42</td></tr><tr><td>RPA [90]</td><td>3.33</td><td></td><td></td><td>11</td><td>3.48</td></tr><tr><td>DMC [91,92]</td><td>3.426</td><td></td><td></td><td>20</td><td>3.63</td></tr><tr><td><strong>DFT</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LDA [86,93–96]</td><td>3.3–3.35</td><td>0.1</td><td>0.327–0.34</td><td>9.3–20</td><td>3.5–3.66</td></tr><tr><td>GGA-D2 [97]</td><td>3.347<sup>b</sup></td><td>0.02</td><td>3.347<sup>b</sup></td><td></td><td></td></tr><tr><td>GGA-D3 [98]</td><td>3.36</td><td></td><td></td><td>10</td><td>3.64</td></tr><tr><td>vdW-DF2 [99]<sup>a</sup></td><td>3.55</td><td>–0.085</td><td>3.55</td><td></td><td></td></tr><tr><td>vdW-DF1-optB88 [100]</td><td>3.343</td><td></td><td></td><td>10.5</td><td>3.528</td></tr><tr><td>vdW-DF1-optB86b [100]</td><td>3.316</td><td></td><td></td><td>11.3</td><td>3.522</td></tr><tr><td>HSE-TS [101]</td><td>3.3<sup>b</sup></td><td></td><td></td><td>18.2</td><td>3.3<sup>b</sup></td></tr><tr><td>HSE-MBD [101]</td><td>3.3<sup>b</sup></td><td></td><td></td><td>17.5</td><td>3.3<sup>b</sup></td></tr><tr><td><strong>MD</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>MLBOP [98]</td><td>3.27</td><td></td><td></td><td>5</td><td>3.66</td></tr><tr><td><em>This work below</em></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><strong>DFT</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>LDA (AIMPRO)</td><td>3.325</td><td>–1.18</td><td>3.314</td><td>10.22</td><td>3.585</td></tr><tr><td>GGA-D2 (AIMPRO)</td><td>3.221</td><td>–0.25</td><td>3.215</td><td>13.33</td><td>3.466</td></tr><tr><td>GGA-D2 (VASP)</td><td>3.193</td><td>0.25</td><td>3.195</td><td>13.64</td><td>3.417</td></tr><tr><td>GGA-TS (VASP)</td><td>3.324</td><td>0.15</td><td>3.325</td><td>13.37</td><td>3.450</td></tr><tr><td><strong>Empirical potentials</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>AIREBO</td><td>3.363</td><td>0.02</td><td>3.363</td><td>0.94</td><td>3.388</td></tr><tr><td>ReaxFF-C2013</td><td>3.260</td><td>0.0</td><td>3.26</td><td>0.53</td><td>3.270</td></tr><tr><td>COMB3</td><td>3.365</td><td>0.0</td><td>3.365</td><td>1.07</td><td>3.393</td></tr><tr><td>LCBOP-I</td><td>3.346</td><td>0.0</td><td>3.346</td><td>0.94</td><td>3.364</td></tr><tr><td><strong>ML potentials</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>GAP20</td><td>3.353</td><td>–4.58</td><td>3.226</td><td>10.01</td><td>3.200</td></tr><tr><td>ACE-D2</td><td>3.402</td><td>0.09</td><td>3.437</td><td>13.90</td><td>3.454</td></tr></tbody></table>

to AB upon heating above 1300 °C [103,104]. The interlayer spacing value for AB graphite ($d_{AB}$) is 3.355 Å [78] and thought to be the same for ABC graphite [102]. The energy difference between the two has been sparsely investigated experimentally but should be quite small given the metastability of ABC graphite. Early work by Boehm and Coughlin [87] calculated an energy difference of 6.1 meV/atom in favour of AB graphite but the authors note that this value is unexpectedly high. More recently, an investigation of FLG by Li et al. [88] found that the preferred stacking order is tuneable through doping and applying an electric field, with the energy difference between the two stacking orders varying by at most 0.15 meV/atom. Pure AA stacked graphite is unstable [95] and is rarely found experimentally except for in intercalation compounds [105], although one work has reported growing pure AA graphite on diamond, measuring an interlayer separation of 3.55 Å [86].

The interlayer distances and relative energies have been calculated for all stacking types using a variety of computational means and are summarised in Table A.1. Beginning with the highest level theory, two studies using the random phase approximation (RPA) method report very similar findings with rather accurate $d_{AB}$ of 3.33–3.34 Å and a significant $\Delta E_{AB→AA}$ of 10–11 meV/atom [89,90] reflecting the instability of AA graphite. Diffusion quantum Monte Carlo (DMC) methods applied to AB [91] and AA [92] graphite, yield interlayer spacings 2% larger than experiment and an energy difference ($\Delta E_{AB→AA}$) of 20 meV/atom, twice that of RPA.

Density functional theory in various forms has been widely applied to these systems. The earliest work published in 1994 [93] uses the

local density approximation (LDA) and reports interlayer distances impressively close to those found experimentally with $d_{AB}$ of 3.35 Å, $d_{ABC}$ of 3.34 Å and $d_{AA}$ of 3.66 Å. The energy differences too are in line with what we would expect with a small difference of 0.1 meV/atom for $\Delta E_{AB→ABC}$, and a much larger difference of 17.3 meV/atom for $\Delta E_{AB→AA}$. Other calculations using LDA have found values which roughly agree with these, though the $\Delta E_{AB→AA}$ values vary from 9.3 to 20 meV/atom [86,94–96]. Nery et al. [97] compared the energy differences between AB and ABC using several functionals and electronic temperatures, finding similar results between LDA and GGA-D2 which both indicated that ABC is the more stable structure for low electronic temperature (284 K) while AB is more stable for higher temperatures. The preference of ABC stacking over AB was also reported in tri-layer graphene by a recent DFT study [99] which utilised a van der Waals functional, vdW-DF2. Two related functionals were used in another study [100] to compare AB and AA graphite, finding similar results to those previously with $\Delta E_{AB→AA}$ around 11 meV/atom. Finally Ouyang et al. [101] assessed the many-body dispersion (MBD) and Tkatchenko-Scheffler (TS) methods, finding $\Delta E_{AB→AA}$ values at the higher end of the range with 17.5 and 18.2 meV/atom respectively.

At the next scale of simulation lies MD where there are several carbon potentials which seek to include the long-range interactions involved in stacking. The Adaptive Intermolecular Reactive Empirical Bond Order potential (AIREBO) [106], Long-range Carbon Bond-Order Potential (LCBOP) [107], the reactive force-field (ReaxFF) [108] and Charge-Optimised-Many-Body (COMB) [109] potentials all improve on the interlayer spacing for graphite from preceding carbon potentials,

<footnote>a Some few-layer graphene results are included where bulk graphite results are not available.</footnote>

<footnote>b Denoted interlayer spacings have been set and not optimised.</footnote>
however they are not able to capture the differences between AB and AA stacking [110]. There are however several Kolmogorov-Crespi [94] based potentials [111,112] which are designed for this problem. Unfortunately for our purposes, these require assignment of each atom to a particular layer which is not possible for systems containing screws which connect the layers. More recently there has been a surge of MLIPs such as the hybrid neural network potential for multilayer graphene (hNN-Grₓ) [110] which was used by McHugh et al. [113] to study ripplocations in graphite where they found hNN-Grₓ performs similarly to the Kolmogorov-Crespi based potentials. A recent pre-print [98] has compared AA and AB graphite using a machine-learning bond-order potential (MLBOP) against DFT calculations using a GGA with D3 correction, finding good agreement in interlayer spacings but an $\Delta E_{AB→AA}$ of around half that of the DFT calculations. Due to a lack of explicit benchmarking of the stacking behaviours of most other potentials for bulk graphites in the literature, these and other MLIPs such as the Gaussian approximation potential (GAP20) [114] and the atomic cluster expansion (ACE) [19,20] will be explored further in the following Appendix B.

The anisotropy of graphite makes the calculation of interlayer interactions challenging even for the most accurate electronic structure theory calculations [115]. Despite this, by comparing with experimental results there are several consistent features we can observe in Table A.1. The most robust methods point to a $\Delta E_{AB→AA}$ of 10–20 meV/atom while $\Delta E_{AB→ABC}$ is two orders of magnitude smaller, at around 0.1 meV/atom. The interlayer spacing for AA graphite is generally found to be larger than that of turbostratic carbon, in agreement with the one experimental value of 3.55 Å. These are the key features required of a MD potential for it to adequately describe the various stacking structures in graphite.

## Appendix B. Benchmarking

To find a suitable potential for modelling dislocations while accurately describing the stacking features, we begin by benchmarking several traditional empirical potentials and comparing these to some newer MLIPs as well as our own DFT calculations. Of the traditional potentials we test AIREBO, ReaxFF with the C2013 reparametrisation [116], COMB and LCBOP-1. These are compared to two MLIPs: GAP20 and ACE with a D2 dispersion correction [19,20]. All potentials are used as implemented in the package LAMMPS [63] version 02Aug23.

For each potential we test the relaxation of three initial systems: AA, AB and ABC graphite, with the AA and AB systems being composed of eight layers totalling 2304 atoms while the ABC system contains nine layers totalling 2592 atoms. Periodic boundary conditions are used in all directions and the relaxations are performed with a looping algorithm that makes use of the conjugate gradient (cg) and Hessian-free truncated Newton (hftn) styles available in LAMMPS. We minimise the energy by adjusting both the atomic coordinates and the unit cell dimensions, aiming to achieve a pressure of 0 GPa on the box. To achieve this we begin with the cg style, alternating between relaxing just the atom positions and the atoms and box. The slower hftn style then replaces the cg method for the atom position relaxation to achieve final small improvements in the energy minimisation. The results are presented in the lower part of Table A.1 and show that although most of the empirical potentials achieve reasonable interlayer spacings for AB graphite, $\Delta E_{AB→ABC}$ are almost all non-existent while $\Delta E_{AB→AA}$ is severely underestimated. Clearly, these traditional potentials are not capable of adequately describing the stacking differences that exist in graphite.

On the other hand, the two MLIPs considered present $\Delta E_{AB→AA}$ within the range suggested by the higher-level theories, although $\Delta E_{AB→ABC}$ for GAP20 suggests ABC is more stable than AB. The interlayer spacings are mixed with GAP20 having a rather accurate $d_{AB}$ while severely underestimating the spacings for ABC and AA, while

<figure><img src="image_13.png" /><figcaption>Fig. B.11. Energy surfaces for different graphitic stacking calculated using GGA-D2 and ACE-D2. Energies are given relative to the AB values for each method. (a) $\gamma$-surface energy for a periodic two-layer unit cell. For each point on the plot the first layer is fixed, the second layer translated in-plane by a fraction of the two $a$ unit cell vectors, then the $c$-axis vector is allowed to relax. (b) A two-atom periodic single-layer unit cell, tilting the $c$-axis by varying $\alpha$ and $\beta$, the angles between $c$ and the two $a$ axis vectors. For each combination of $\alpha$ and $\beta$, the magnitude of $c$ is varied to find the energetic minimum. The idealised configurations (AA, SP, ABC and AB) are coloured according to the stacking registration scale depicted in Fig. 4(b). For the full registration measure projection of the ACE-D2 maps see Supplementary Materials.</figcaption></figure>

ACE has interlayer spacings which are too large, but do exhibit the correct relative differences between stacking types. As discussed by Qamar et al. [20], ACE was trained on DFT data computed with the Perdew–Burke–Ernzerhof (PBE) GGA without an empirical correction, allowing the flexibility to be used with various empirical dispersion corrections. They showed ACE-D2 reproduces PBE-D2 results and also yields a more accurate interlayer separation than using either a D3 or D4 correction. Despite these interlayer discrepancies, these are of less relevance compared to the key importance of computing the energy surfaces between AB, ABC and AA stacking accurately. Qian et al. reports (in Figure 7(e) of Ref. [117]) that the shape of the $\gamma$-surface energy for GAP20 sliding a bilayer between AA and AB stacking does not agree with DFT. To confirm ACE as a suitable potential for our investigation we calculate the $\gamma$-surface energy for translation and tilting (to achieve ABC stacking) and compare this to DFT calculations.

DFT calculations are performed using the AIMPRO code [118]. Carbon atoms are modelled using HGH pseudopotentials [119] and 38 independent Gaussian-based functions per atom, with a 14 × 14 × 4 k-point mesh [120] for the unit cell. Both, generalised gradient approximation (GGA) [121] with D2 dispersion [122] and the LDA [123] were tested for the exchange correlation functional. The GGA energetics and interlayer spacing are highly sensitive to the calculation setup. As observed by Nery et al. [97], we also find ABC to be very slightly ($\Delta E_{AB→ABC} = -0.25$ meV/atom) more stable than AB-stacking with the low electron temperature used here (kT = 0.02 eV); raising kT to 0.2 eV inverts the ordering with AB more stable by 0.07 meV. The interlayer spacing is highly dependent on the choice of dispersion correction used; a series of tests using Vienna Ab initio Simulation Package
(VASP) [124–127] (see Table A.1) show the closest interlayer spacings to experiment are found with the TS dispersion correction [128], while Grimme-D2 typically underestimates the spacing by ~0.19 Å (see Supplementary Materials).

The γ-surface energy landscape is mapped for GGA-D2 calculated with AIMPRO starting from a 4-atom AB bilayer graphite cell, then translating one of the two layers in a series of steps along the basal directions [10̅10] and [01̅0]. The *a*-lattice parameter is held fixed at the AB optimised value (a = 2.463 Å) and the *c*-spacing allowed to relax as shown in Fig. B.11(a). In order to map the shift to ABC graphite, a 2-atom AA-stacked cell is used; this time the *c*-vector is tilted by varying α and β, the two angles between the *c*- and *a*-vectors. For each fixed pair of angles, the length of *c* is varied to find the energetic minimum as shown in Fig. B.11(b). Repeating the same procedure with ACE-D2 shows remarkably good reproduction of both energy surfaces, confirming it as a suitable potential to accurately reproduce the energetics of interlayer stacking.

## Appendix C. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.carbon.2025.120995.

## References

[1] M. Li, J. Lu, Z. Chen, K. Amine, 30 years of lithium-ion batteries, Adv. Mater. 30 (2018) 1800561.

[2] Q. Dai, J.C. Kelly, J. Dunn, P.T. Benavides, Update of Bill-of-Materials and Cathode Materials Production for Lithium-ion Batteries in the GREET® Model, Argonne National Laboratory, 2018, https://greet.anl.gov/files/update_bom_cm.

[3] R.L. Moss, E. Tzimas, P. Willis, J. Arendorf, L. Tercero Espinoza, et al., Critical Metals in the Path Towards the Decarbonisation of the EU Energy Sector: Assessing Rare Metals as Supply-Chain Bottlenecks in Low-Carbon Energy Technologies, European Commission, Publications Office of the European Union, 2013, https://op.europa.eu/en/publication-detail/-/publication/505c089c-7655-4546-bd17-83f9d581190/language-en.

[4] D.W. Olson, 2018 Minerals Handbook - Graphite, National Minerals Information Center U.S. Geological Survey, 2022.

[5] J. Zhang, C. Liang, J.B. Dunn, Graphite flows in the US: Insights into a key ingredient of energy transition, Environ. Sci. Technol. 57 (2023) 3402–3414.

[6] W. Bollmann, Electron-microscopic observations on radiation damage in graphite, Phil. Mag. 5 (54) (1960) 621–624.

[7] P. Thrower, Interstitial and vacancy loops in graphite irradiated at high temperatures, Br. J. Appl. Phys. 15 (1964) 1153.

[8] R. Telling, M. Heggie, Radiation defects in graphite, Phil. Mag. 87 (31) (2007) 4797–4846.

[9] P. Thrower, R. Mayer, Point defects and self-diffusion in graphite, Phys. Status Solidi, A, Appl. Res. 47 (1) (1978) 11–37.

[10] S. Snyder, W.W. Gerberich, H. White, Scanning-tunneling-microscopy study of tip-induced transitions of dislocation-network structures on the surface of highly oriented pyrolytic graphite, Phys. Rev. B 47 (16) (1993) 10823.

[11] S. Amelinckz, P. Delavignette, M. Heerschap, Dislocations and stacking faults in graphite, in: Chemistry and Physics of Carbon vol. 1, Marcel Dekker, 1965, pp. 1–69.

[12] J. Rothstein, Possible macroscopic effects of single lattice defects, Phys. Rev. 95 (2) (1954) 370.

[13] I. Suarez-Martinez, Theory of Diffusion and Plasticity in Layered Carbon Materials (Ph.D. thesis), University of Sussex, 2007.

[14] M. Chen, B. Wu, L. Zhou, Y. Zhu, H. Wu, Micromechanical properties of pyrolytic carbon with interlayer crosslink, Carbon 159 (2020) 549–560.

[15] J.-M. Leyssale, G. Couégnat, S. Jouannigot, G.L. Vignoles, Mechanisms of elastic softening in highly anisotropic carbons under in-plane compression/indentation, Carbon 197 (2022) 425–434.

[16] F. Polewczyk, P. Lafourcade, N. Pineau, C. Denoual, G.L. Vignoles, J.-M. Leyssale, A molecular dynamics investigation of laminar pyrocarbons elasticity up to high-temperatures, Carbon Trends 16 (2024) 100376.

[17] J.W. Martin, J.L. Fogg, G.R. Francas, K.J. Putman, E.P. Turner, I. Suarez-Martinez, N.A. Marks, Graphite rapidly forms via annihilation of screw dislocations, Carbon 215 (2023) 118386.

[18] V. Meunier, G. Bepete, M.-S. Cao, Y. Chen, C. de Tomas, J. Di, C. Ewels, N. Koratkar, Q. Li, C. Liu, et al., Carbon science perspective in 2024: Current research and future challenges, Carbon 229 (2024) 119488.

[19] R. Drautz, Atomic cluster expansion for accurate and transferable interatomic potentials, Phys. Rev. B 99 (1) (2019) 014104.

[20] M. Qamar, M. Mrovec, Y. Lysogorskiy, A. Bochkarev, R. Drautz, Atomic cluster expansion for quantum-accurate large-scale simulations of carbon, J. Chem. Theory Comput. 19 (15) (2023) 5151–5167.

[21] F.E. Fujita, K. Izui, Observation of lattice defects in graphite by electron microscopy, part I, J. Phys. Soc. Japan 16 (1961) 214–227.

[22] T. Mitran, G. Nemnes, Helical graphite metamaterials for intense and locally controllable magnetic fields, RSC Adv. 7 (2017) 49041–49047.

[23] G. Hennig, Reactivity changes near lattice defects, in: Chemistry and Physics of Carbon vol. 2, Marcel Dekker, 1966, pp. 1–50.

[24] J. Rakovan, J.A. Jaszczak, Multiple length scale growth spirals on metamorphic graphite {001} surfaces studied by atomic force microscopy, Am. Mineral. 87 (2002) 17–24.

[25] F.H. Horn, Spiral growth on graphite, Nature 170 (4327) (1952) 581–581.

[26] A. Patel, O. Bahl, Evidence of screw dislocations in graphite, Br. J. Appl. Phys. 16 (2) (1965) 169.

[27] G. Hennig, Screw dislocations in graphite, Science 147 (1965) 733–734.

[28] C. Roscoe, J.M. Thomas, The identification and some physico-chemical consequences of non-basal edge and screw dislocations in graphite, Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 297 (1967) 397–407.

[29] H. Badenhorst, Microstructure of natural graphite flakes revealed by oxidation: Limitations of XRD and Raman techniques for crystallinity estimates, Carbon 66 (2014) 674–690.

[30] C. Roscoe, J. Baker, Electron microscope studies of the deviation of a stress- annealed pyrolytic graphite matrix from a single crystal, J. Appl. Phys. 40 (4) (1969) 1665–1673.

[31] J.J. Baker, Microstructural Defects and Their Importance in Highly Crystalline Pyrolytic Graphites (Ph.D. thesis), The Pennsylvania State University, 1970.

[32] Y. Sun, L.B. Alemany, W. Billups, J. Lu, B.I. Yakobson, Structural dislocations in anthracite, J. Phys. Chem. Lett. 2 (2011) 2521–2524.

[33] K. Li, Q. Liu, D. Hou, Z. Wang, S. Zhang, Quantitative investigation on the structural characteristics and evolution of high-rank coals from Xinhua, Hunan Province, China, Fuel 289 (2021) 119945.

[34] X. Xu, Y. Qi, L. Tian, J. Dou, J. Appiah, J. Wang, R. Zhu, X. Chen, J. Yu, Evolution of carbon structure and cross-linking structure during coking process of coking coal, J. Anal. Appl. Pyrolysis 182 (2024) 106670.

[35] A. Oberlin, Carbonization and graphitization, Carbon 22 (1984) 521–541.

[36] M. Ghazinejad, S. Holmberg, O. Pilloni, L. Oropeza-Ramos, M. Madou, Graphitizing non-graphitizable carbons by stress-induced routes, Sci. Rep. 7 (1) (2017) 16551.

[37] S. Tulić, T. Waitz, M. Caplovicova, G. Habler, M. Varga, M. Kotlár, V. Vretenar, O. Romanyuk, A. Kromka, B. Rezek, et al., Covalent diamond–graphite bonding: mechanism of catalytic transformation, ACS Nano 13 (4) (2019) 4621–4630.

[38] Y. Zhao, S. Jin, Stacking and twisting of layered materials enabled by screw dislocations and non-Euclidean surfaces, Acc. Mater. Res. 3 (3) (2022) 369–378.

[39] I. Minkoff, W. Nixon, Scanning electron microscopy of graphite growth in iron and nickel alloys, J. Appl. Phys. 37 (13) (1966) 4848–4855.

[40] D. Double, A. Hellawell, Cone-helix growth forms of graphite, Acta Metall. 22 (4) (1974) 481–487.

[41] J. Lacaze, K. Theuwissen, L. Laffont, M. Véron, Misorientations in spheroidal graphite: some new insights about spheroidal graphite growth in cast irons, in: IOP Conference Series: Materials Science and Engineering, 117, (1) IOP Publishing, 2016, 012024.

[42] D.M. Stefanescu, G. Alonso, R. Suarez, Recent developments in understanding nucleation and crystallization of spheroidal graphite in iron-carbon-silicon alloys, Metals 10 (2) (2020) 221.

[43] R.Y. Tay, H.J. Park, J. Lin, Z.K. Ng, L. Jing, H. Li, M. Zhu, S.H. Tsang, Z. Lee, E.H.T. Teo, Concentric and spiral few-layer graphene: growth driven by interfacial nucleation vs screw dislocation, Chem. Mater. 30 (2018) 6858–6866.

[44] P. Pochet, H.T. Johnson, A double-helix dislocation in graphene, Nat. Mater. 23 (3) (2024) 306–307.

[45] Z.-J. Wang, X. Kong, Y. Huang, J. Li, L. Bao, K. Cao, Y. Hu, J. Cai, L. Wang, H. Chen, et al., Conversion of chirality to twisting via sequential one-dimensional and two-dimensional growth of graphene spirals, Nat. Mater. 23 (3) (2024) 331–338.

[46] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556 (7699) (2018) 43–50.

[47] F. Xu, H. Yu, A. Sadrzadeh, B.I. Yakobson, Riemann surfaces of carbon as graphene nanosolenoids, Nano Lett. 16 (2016) 34–39.

[48] F.S. Boi, S. Wang, W. Li, O. Odunmbaku, S. Gao, J. Guo, J. Song, J. Wen, H. Zhang, Interplay of edge/screw dislocations and enhanced magnetism in exfoliated pyrolytic graphite with distorted hexagonal Moiré superlattices, Carbon Trends 5 (2021) 100106.

[49] Y. Zhang, B. Xie, Y. Yang, Y. Wu, X. Lu, Y. Hu, Y. Ding, J. He, P. Dong, J. Wang, et al., Extremely large magnetoresistance in twisted intertwined graphene spirals, Nat. Commun. 15 (1) (2024) 6120.

[50] I. Suarez-Martinez, G. Savini, M. Heggie, First principles modelling of scroll-to-nanotube defect: Screw-type dislocation, in: Materials Science Forum, vol. 527, Trans Tech Publ, 2006, pp. 1583–1586.
[51] I. Suarez-Martinez, G. Savini, G. Haffenden, J.-M. Campanera, M. Heggie, Dislocations of Burgers vector c/2 in graphite, Phys. Status Solidi C 4 (2007) 2958-2962.

[52] A. Kuc, G. Seifert, Hexagon-preserving carbon foams: Properties of hypothetical carbon allotropes, Phys. Rev. B—Condensed Matter Mater. Phys. 74 (21) (2006) 214104.

[53] M.-k. Choi, H.T. Johnson, Mechanics of out-of-plane screw dislocation in a 2D material, Extrem. Mech. Lett. 77 (2025) 102320.

[54] J.-M. Leyssale, J.-P. Da Costa, C. Germain, P. Weisbecker, G. Vignoles, Structural features of pyrocarbon atomistic models constructed from transmission electron microscopy images, Carbon 50 (2012) 4388-4400.

[55] F. Polewczyk, P. Lafourcade, J.-P. Da Costa, G. Vignoles, J.-M. Leyssale, Polygranular image guided atomistic reconstruction: A parametric model of pyrocarbon nanostructure, Carbon 212 (2023) 118109.

[56] A. Chartier, L. Van Brutzel, J. Pageot, Irradiation damage in nuclear graphite at the atomic scale, Carbon 133 (2018) 224-231.

[57] G.R. Francas, J.W. Martin, I. Suarez-Martinez, N.A. Marks, Topological defects and anisotropic development during pre-graphitization, Carbon (2023) 118251.

[58] G.R. Francas, J.W. Martin, K.J. Putman, I. Suarez-Martinez, N.A. Marks, Mechanisms and kinetics of martensic dislocation loop removal during graphitization, Carbon 229 (2024) 119400.

[59] I. Suarez-Martinez, G. Savini, A. Zobell, M. Heggie, Dislocations in carbon nanotube walls, J. Nanosci. Nanotechnol. 7 (10) (2007) 3417-3420.

[60] J.Y. Huang, F. Ding, B. Yakobson, Dislocation dynamics in multiwalled carbon nanotubes at high temperatures, Phys. Rev. Lett. 100 (2008) 035503.

[61] E. Akatyeva, J. Huang, T. Dumitrică, Edge-mediated dislocation processes in multishell carbon nano-onions? Phys. Rev. Lett. 105 (2010) 106102.

[62] P. Hirel, Atomsk: A tool for manipulating and converting atomic data files, Comput. Phys. Comm. 197 (2015) 212-219.

[63] A.P. Thompson, H.M. Aktulga, R. Berger, D.S. Bolintineanu, W.M. Brown, P.S. Crozier, P.J. in't Veld, A. Kohlmeier, S.G. Moore, T.D. Nguyen, R. Shan, M.J. Stevens, J. Tranchida, C. Trott, S.J. Plimton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comp. Phys. Comm. 271 (2022) 108171, http://dx.doi.org/10.1016/j.cpc.2021.108171.

[64] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, P. Gumbsch, Structural relaxation made simple, Phys. Rev. Lett. 97 (17) (2006) 170201.

[65] J. Guénolé, W.G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash, E. Bitzek, Assessment and optimization of the fast inertial relaxation engine (fire) for energy minimization in atomistic simulations and its implementation in lammps, Comput. Mater. Sci. 175 (2020) 109584.

[66] A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO-the open visualization tool, Model. Simul. Mater. Sci. Eng. 18 (1) (2010) http://dx.doi.org/10.1088/0965-0393/18/1/015012.

[67] N. Marom, J. Bernstein, J. Garel, A. Tkatchenko, E. Joselevich, L. Kronik, O. Hod, Stacking and registry effects in layered materials: the case of hexagonal boron nitride, Phys. Rev. Lett. 105 (4) (2010) 046801.

[68] O. Hod, Quantifying the stacking registry matching in layered materials, Isr. J. Chem. 50 (4) (2010) 506-514.

[69] O. Hod, The registry index: A quantitative measure of materials interfacial commensurability, ChemPhysChem 14 (11) (2013) 2376-2391.

[70] W. Cao, O. Hod, M. Urbakh, Interlayer registry index of layered transition metal dichalcogenides, J. Phys. Chem. Lett. 13 (15) (2022) 3353-3359.

[71] J.D. Morrow, RingsStatisticsMatter, 2023, https://github.com/MorrowChem/RingsStatisticsMatter.jl.

[72] X. Yuan, A. Cormack, Efficient algorithm for primitive ring statistics in topological networks, Comput. Mater. Sci. 24 (3) (2002) 343-360.

[73] J. Martin, C. de Tomas, I. Suarez-Martinez, M. Kraft, N. Marks, Topology of disordered 3D graphene networks, Phys. Rev. Lett. 123 (2019) 116105.

[74] R.E. Franklin, Crystallite growth in graphitizing and non-graphitizing carbons, Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 209 (1951) 196-218.

[75] R.E. Franklin, The structure of graphitic carbons, Acta Crystallogr. 4 (3) (1951) 253-261.

[76] P. Kokmat, P. Surinlert, A. Ruammaitree, Growth of high-purity and high-quality turbostratic graphene with different interlayer spacings, ACS Omega 8 (4) (2023) 4010-4018.

[77] J. Sha, S. Wang, X. Zhang, J. Zhang, Evidence of a connection between turbostratic structure and twisted graphene, Cryst. Growth Des. (2025).

[78] J.Y. Howe, C.J. Rawn, L. Jones, H. Ow, Improved crystallographic data for graphite, Powder Diffr. 18 (2) (2003) 150-154.

[79] J. Pappis, S. Blum, Properties of pyrolytic graphite, J. Am. Ceram. Soc. 44 (12) (1961) 592-597.

[80] C. Klein, M. Holland, Thermal conductivity of pyrolytic graphite at low temperatures. I. Turbostratic structures, Phys. Rev. 136 (2A) (1964) A575.

[81] M. Monthiox, M. Oberlin, A. Oberlin, X. Bourrat, R. Boulet, Heavy petroleum products: microtexture and ability to graphitize, Carbon 20 (3) (1982) 167-176.

[82] M. Monthiox, 4. Structure, texture, and thermal behaviour of polyaromatic solids, in: Carbon Molecules and Materials, CRC Press, 2002, pp. 127-178.

[83] M. Bonijoly, M. Oberlin, A. Oberlin, A possible mechanism for natural graphite formation, Int. J. Coal Geol. 1 (4) (1982) 283-312.

[84] B. Yang, M. Barsoum, R. Rethinam, Nanoscale continuum calculation of basal dislocation core structures in graphite, Phil. Mag. 91 (10) (2011) 1441-1463.

[85] P. Mouratidis, J.G. McHugh, K. Jolley, Modelling of partial basal dislocation dipoles in bilayer graphene and graphite, Carbon 217 (2024) 118613.

[86] J.-K. Lee, S.-C. Lee, J.-P. Ahn, S.-C. Kim, J.I. Wilson, P. John, The growth of AA graphite on (111) diamond, J. Chem. Phys. 129 (23) (2008).

[87] H. Boehm, R. Coughlin, Enthalpy difference of hexagonal and rhombohedral graphite, Carbon 2 (1) (1964) 1-6.

[88] H. Li, M.I.B. Utama, S. Wang, W. Zhao, S. Zhao, X. Xiao, Y. Jiang, L. Jiang, T. Taniguchi, K. Watanabe, et al., Global control of stacking-order phase transition by doping and electric field in few-layer graphene, Nano Lett. 20 (5) (2020) 3106-3112.

[89] S. Lebègue, J. Harl, T. Gould, J. Ángyán, G. Kresse, J. Dobson, Cohesive properties and asymptotics of the dispersion interaction in graphite by the random phase approximation, Phys. Rev. Lett. 105 (19) (2010) 196401.

[90] N. Leconte, J. Jung, S. Lebègue, T. Gould, Moiré-pattern interlayer potentials in van der Waals materials in the random-phase approximation, Phys. Rev. B 96 (19) (2017) 195431.

[91] L. Spanu, S. Sorella, G. Galli, Nature and strength of interlayer binding in graphite, Phys. Rev. Lett. 103 (19) (2009) 196401.

[92] P. Ganesh, J. Kim, C. Park, M. Yoon, F.A. Reboredo, P.R. Kent, Binding and diffusion of lithium in graphite: quantum Monte Carlo benchmarks and validation of van der Waals density functional methods, J. Chem. Theory Comput. 10 (12) (2014) 5318-5323.

[93] J.-C. Charlier, X. Gonze, J.-P. Michenaud, First-principles study of the stacking effect on the electronic properties of graphite (s), Carbon 32 (2) (1994) 289-299.

[94] A.N. Kolmogorov, V.H. Crespi, Registry-dependent interlayer potential for graphitic systems, Phys. Rev. B—Condensed Matter Mater. Phys. 71 (23) (2005) 235415.

[95] G. Savini, Y. Dappe, S. Öberg, J.-C. Charlier, M. Katsnelson, A. Fasolino, Bending modes, elastic constants and mechanical stability of graphitic systems, Carbon 49 (1) (2011) 62-69.

[96] P. Anees, M. Valsakumar, S. Chandra, B. Panigrahi, Ab initio study on stacking sequences, free energy, dynamical stability and potential energy surfaces of graphite structures, Modelling Simul. Mater. Sci. Eng. 22 (3) (2014) 035016.

[97] J.P. Nery, M. Calandra, F. Mauri, Ab-initio energetics of graphite and multilayer graphene: stability of bernal versus rhombohedral stacking, in: 2D Materials, vol. 8 (2021) 035006, 3.

[98] I. Kohata, K. Hisama, K. Otsuka, S. Maruyama, A machine-learning bond-order potential for exploring the configuration space of carbon, 2025, arXiv preprint arXiv:2501.11297.

[99] R. Guerrero-Avilés, M. Pelc, F.R. Geisenhof, R.T. Weitz, A. Ayuela, Rhombohedral trilayer graphene is more stable than its bernal counterpart, Nanoscale 14 (43) (2022) 16295-16302.

[100] R.R. Del Grande, M.G. Menezes, R.B. Capaz, Layer breathing and shear modes in multilayer graphene: a DFT-vdW study, J. Phys.: Condens. Matter. 31 (29) (2019) 295301.

[101] W. Ouyang, I. Azuri, D. Mandelli, A. Tkatchenko, L. Kronik, M. Urbakh, O. Hod, Mechanical and tribological properties of layered materials under high pressure: assessing the importance of many-body dispersion effects, J. Chem. Theory Comput. 16 (1) (2019) 666-676.

[102] H.S. Lipson, A. Stokes, The structure of graphite, Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 181 (984) (1942) 101-105.

[103] P. Hawtin, J. Lewis, N. Moul, R. Phillips, The heats of combustion of graphite, diamond and some non-graphitic carbons, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Sci. 261 (1116) (1966) 67-95.

[104] F. Laves, Y. Baskin, On the formation of the rhombohedral graphite modification, Z. Krist. 107 (5-6) (1956) 337-356.

[105] M.S. Dresselhaus, G. Dresselhaus, Intercalation compounds of graphite, Adv. Phys. 30 (2) (1981) 139-326.

[106] S.J. Stuart, A.B. Tutein, J.A. Harrison, A reactive potential for hydrocarbons with intermolecular interactions, J. Chem. Phys. 112 (14) (2000) 6472-6486.

[107] J. Los, A. Fasolino, Intrinsic long-range bond-order potential for carbon: Performance in Monte Carlo simulations of graphitization, Phys. Rev. B 68 (2) (2003) 024107.

[108] A.C. Van Duin, S. Dasgupta, F. Lorant, W.A. Goddard, Reaxff: a reactive force field for hydrocarbons, J. Phys. Chem. A 105 (41) (2001) 9396-9409.

[109] T. Liang, T.-R. Shan, Y.-T. Cheng, B.D. Devine, M. Noordhoek, Y. Li, Z. Lu, S.R. Phillpot, S.B. Sinnott, Classical atomistic simulations of surfaces and heterogeneous interfaces with the charge-optimized many body (COMB) potentials, Mater. Sci. Eng.: R: Rep. 74 (9) (2013) 255-279.

[110] M. Wen, E.B. Tadmor, Hybrid neural network potential for multilayer graphene, Phys. Rev. B 100 (19) (2019) 195419.

[111] I.V. Lebedeva, A.A. Knizhnik, A.M. Popov, Y.E. Lozovik, B.V. Potapkin, Interlayer interaction and relative vibrations of bilayer graphene, Phys. Chem. Chem. Phys. 13 (13) (2011) 5687-5695.

[112] M. Wen, S. Carr, S. Fang, E. Kaxiras, E.B. Tadmor, Dihedral-angle-corrected registry-dependent interlayer potential for multilayer graphene structures, Phys. Rev. B 98 (23) (2018) 235404.
[113] J.G. McHugh, P. Mouratidis, K. Jolley, Ripplocations in layered materials: Sublinear scaling and basal climb, Phys. Rev. B 103 (2021) 195436.

[114] P. Rowe, V.L. Deringer, P. Gasparotto, G. Csányi, A. Michaelides, An accurate and transferable machine learning potential for carbon, J. Chem. Phys. 153 (3) (2020).

[115] J.F. Dobson, Dispersion and induction interactions of graphene with nanostructures, Surf. Sci. 605 (17-18) (2011) 1621-1632.

[116] S.G. Srinivasan, A.C. Van Duin, P. Ganesh, Development of a ReaxFF potential for carbon condensed phases and its application to the thermal fragmentation of a large fullerene, J. Phys. Chem. A 119 (4) (2015) 571-580.

[117] C. Qian, B. McLean, D. Hedman, F. Ding, A comprehensive assessment of empirical potentials for carbon materials, APL Mater. 9 (6) (2021).

[118] R. Jones, P. Briddon, Identification of defects in semiconductors, Semicond. Semimetals 51 (1998) 287-349.

[119] C. Hartwigsen, S. Goedecker, J. Hutter, Relativistic separable dual-space Gaussian pseudopotentials from H to Rn, Phys. Rev. B 58 (1998) 3641-3662, http://dx.doi.org/10.1103/PhysRevB.58.3641.

[120] H.J. Monkhorst, J.D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B 13 (12) (1976) 5188.

[121] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (18) (1996) 3865.

[122] S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comput. Chem. 27 (15) (2006) 1787-1799.

[123] J.P. Perdew, Y. Wang, Accurate and simple analytic representation of the electron-gas correlation energy, Phys. Rev. B 45 (23) (1992) 13244.

[124] G. Kresse, J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47 (1) (1993) 558.

[125] G. Kresse, J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal-amorphous-semiconductor transition in germanium, Phys. Rev. B 49 (20) (1994) 14251.

[126] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54 (16) (1996) 11169.

[127] G. Kresse, J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6 (1) (1996) 15-50.

[128] A. Tkatchenko, M. Scheffler, Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data, Phys. Rev. Lett. 102 (7) (2009) 073005.