Phase transitions between carbon allotropes are calculated using the generalized solid-state nudged elastic band method. We find a new reaction mechanism between graphite and diamond with nucleation characteristics that has a lower activation energy than the concerted mechanism. The calculated barrier from graphite to hexagonal diamond is lower than to cubic diamond, resolving a conflict between theory and experiment. Transitions are calculated to three structures of cold compressed graphite: bct C4, M, and Z-carbon, which are accessible at the experimentally relevant pressures near 17 GPa.
Carbon has many crystalline phases, including graphite, cubic diamond (CD), hexagonal diamond (HD), and cold compressed graphite (CCG). The transition mechanisms between these phases are not thoroughly understood. There is no consensus, for example, of the structure of CCG.1–3 The transitions from graphite to CD and to HD have, on the other hand, been well documented in experiments. It has been found that both transitions can occur at 15 GPa static pressure and the temperature required to form HD is lower than CD,4,5 indicating that the barrier to form HD is lower at this pressure.
Challenges still exist in the theoretical explanation of these two-phase transitions. Molecular dynamics studies are limited by the vibrational time scale of atoms, and so are often done at pressures several times higher than in experiments so that a transition, which is a rare event, can be observed in the accessible simulation time scale. The problem with this approach is that the transition mechanism can change with increasing pressure.6 Another approach is to use transition state theory to calculate the barrier and mechanism for the activated process. Previous saddle point calculations have focused on the concerted mechanism, where the barrier to HD is higher than to CD, in disagreement with experiments.7,8 The concerted mechanism can be active at pressures higher than 15 GPa or under shock wave compression,9,10 but at 15 GPa static pressure there must be a lower energy pathway. Employing a neural network potential and including 145000 atoms, Khaliullin et al. found the structure of a stable nucleus, but not a reaction pathway of its formation.11 In this simulation, the lowest pressure required to stabilize a diamond nucleus in graphite was 30 GPa, twice the experimental pressure. CCG can also be synthesized at low pressure (17 GPa) and at lower temperatures than is required to form diamond, indicating that a different mechanism with a lower barrier is present.3 To resolve these issues, we use a computational method that can reveal the mechanism of solid-solid phase transitions between carbon phases in the low-pressure regime.
In this Communication, we present a new transition mechanism at low pressure that was found using the generalized solid-state nudged elastic band (G-SSNEB) method. Forces and energies were calculated with density functional theory (DFT) using a simulation cell containing only 40 atoms. The rate-limiting step in this mechanism corresponds to the nucleation of several diamond layers in the graphite structure. The diamond phase then grows, layer by layer, overcoming a series of smaller barriers. Our result shows that the highest barrier to HD is lower than to CD at 15 GPa. The barriers to form previously proposed structures of CCG are compared, and found to be lower than to CD though still higher than to HD, pointing to a selection criterion for identifying the structure of CCG.
Minimum energy paths (MEPs) were calculated using the climbing-image G-SSNEB method.12 This nascent method is appropriate for solid-solid phase transitions that are described by changes in both atomic coordinates and lattice vectors. Neglecting to include either of these degrees of freedom in the reaction coordinate (e.g., by minimizing them along a reaction path) leads to a bias in the calculation and a risk of missing the correct MEP.12 An adaptive nudged elastic band approach was used to increase the resolution near saddle points.13
Hexagonal graphite (HG) was chosen as the initial state in the formation of HD, and rhombohedral graphite (RG) for CD, based upon symmetry considerations.11 Supercells of 1 × 1 × 4 unit cells were used for the HG-HD transition and 2 × 1 × 3 for RG-CD, so that both transitions were represented by four carbon atoms per layer and a similar number of layers (8 for HG and 9 for RG).
Energies and forces between the carbon atoms were calculated with DFT using a plane wave basis set for valence electrons and the projector augmented wave method for core electrons, as implemented in the Vienna ab initio simulation package.14 The electron exchange-correlation energy was calculated in the generalized gradient approximation (GGA) with the Perdew-Wang (PW91) functional.15 Two levels of convergence were used, a lower level with the standard carbon pseudopotential and a plane wave energy cutoff of 520 eV, and then a higher level to converge the enthalpy barriers using hard pseudopotentials and a cutoff energy of 910 eV. A semi-empirical dispersion term was added to describe van der Waals interactions.16
Figure 1 shows the enthalpy landscapes of transitions from RG to CD and from HG to HD at 15 GPa. The nucleation mechanisms (green lines) have lower barriers than the concerted mechanisms (red lines). Furthermore, in the nucleation mechanism, the overall barrier to form CD is higher than HD. The common feature of the nucleation mechanisms in Fig. 1 is that there are a series of barriers and minima along the path, and the initial barrier is higher than the following ones. This landscape is characteristic of a nucleation process followed by growth of the nucleus. Unlike the concerted mechanism, where all the atoms in the supercell transform to the new phase simultaneously, the nucleation mechanism involves the transformation of several layers into a two-dimensional diamond nucleus that subsequently grows layer-by-layer into the graphite. We note, however, that the nucleation and growth occurs along one dimension; the transition is still concerted within each layer. The lateral interfaces, which must also contribute to the energy of a truly local nucleus, are not considered in our small simulation cells. The nucleation mechanism found here may be observable in the case of a small graphite particle where the nucleus is able to extend across the particle so that the lateral interfaces are eliminated. Despite not providing true nucleation barriers, our calculations can be used to compare barriers to different phases when the lateral sizes of the nuclei are comparable. We note that the barriers of the nucleation mechanisms are largely determined by the enthalpy to form the nuclei, C1 and H1.
Enthalpy landscapes and nucleation structures of (a) RG-CD and (b) HG-HD phase transitions at 15 GPa. Arrows on the initial structures indicate the relative movement of atoms. Two unit cells are shown for each structure. The undulations in the minimum energy paths for the nucleation mechanisms correspond to layer-by-layer growth of the diamond phase.
Enthalpy landscapes and nucleation structures of (a) RG-CD and (b) HG-HD phase transitions at 15 GPa. Arrows on the initial structures indicate the relative movement of atoms. Two unit cells are shown for each structure. The undulations in the minimum energy paths for the nucleation mechanisms correspond to layer-by-layer growth of the diamond phase.
For the RG to CD transition, the nucleus is a mixed CD and M-carbon phase (or W-carbon, as one cannot distinguish from one layer). In each basin along the minimum enthalpy landscape, both ends of the nucleus are capped with M-carbon; on the plateaus, one side is M-carbon and the other is a CD/RG interface. The difference in termination means that the M/RG interface has lower energy than the CD/RG interface, although the CD/RG interface is unavoidable as the CD phase grows. There is no such mixed interface for the HG to HD transition so the HD/HG interface is more stable. While the enthalpy of HD is slightly higher than CD, the energy of the HD/HG interface is much lower than that of CD/RG and the overall nucleation barrier is also lower.
Due to our limited supercell size, the nuclei have periodic structure in the graphite planes and the interface between phases is only present between layers. As the nucleus grows, the interface size and structure does not change. The enthalpy change along the reaction path is only due to differences in enthalpy of the graphite and diamond phases. When the pressure is sufficiently high so that the diamond phase is more stable than graphite, the enthalpy profile is exothermic after the nucleus forms.
Since the highest barrier along each MEP is due to the formation of the smallest nucleus, only this barrier is calculated when considering transformations from graphite to CCG. Initial structures of the CCG nuclei were generated based upon those found for CD and HD. Each was taken as a final state of a G-SSNEB calculation. If an intermediate state was found by the G-SSNEB, it was set as the final state and the calculation was restarted. This procedure was repeated until no intermediate state was found, so that the stable nucleus was directly connected to the initial state by a single barrier.
Bct C4, M, and Z-carbon are calculated as candidates of the kinetically accessible CCG phase.17–20 The configurations of the CCG nuclei at 15 GPa are shown in Fig. 2. The structure of the Z-carbon nucleus is only slightly distorted from bulk. Two C atoms above and below the octagon in the center are displaced away from the 4-carbon ring positions, so the 6 + 4 + 6 rings merge into large 12-carbon rings. Interestingly, the distortion does not occur in bct C4, where no 6-carbon rings exist. A Brønsted-Evans-Polanyi (BEP) relationship, showing a linear trend between the barrier and enthalpy of forming the smallest nucleus, is given in Fig. 3. To further save computational effort, we use the enthalpies of the critical nuclei and the BEP relationship to determine the overall phase transition barriers.
The smallest nuclei of bct C4, M, and Z-carbon. Two unit cells are shown for bct C4 and M-carbon; one unit cell is shown for Z-carbon.
The smallest nuclei of bct C4, M, and Z-carbon. Two unit cells are shown for bct C4 and M-carbon; one unit cell is shown for Z-carbon.
The phase transition barriers are shown as a function of pressure in Fig. 4. The barriers are only plotted at pressures where the new phases are more stable than graphite. The three forms of CCG become stable above 15–25 GPa, consistent with experiment and recent calculations also using GGA functionals with long-range dispersion.21 Earlier local density approximation (LDA) calculations underestimate the enthalpy differences between graphite and CCG.8,20
Pressure dependence of the maximum barrier for each phase transition. Dashed lines indicate regions where the new phases are less stable than graphite.
Pressure dependence of the maximum barrier for each phase transition. Dashed lines indicate regions where the new phases are less stable than graphite.
The barriers from RG and HG to the M-phase of CCG are the same because RG and HG differ only in their long-range stacking of layers, which does not effect the formation energy of the critical nucleus. The barrier to Z and M-carbon are similar, although Z-carbon has a somewhat more negative activation volume (the slope in Fig. 4), so the transition to Z becomes favorable as compared to M above 20 GPa. An interesting and important result is that enthalpy changes to the final phases are not proportional to the barriers to reach them. It is the formation energy of the critical nucleus that determines the activation energy. For example, bct C4 has higher enthalpy than M-carbon,21 but the barrier is lower due to a lower energy interface with graphite. In the formation of metastable forms of carbon, the kinetics can be more important than the thermodynamics.
Figure 4 shows that the formation of CCG is favorable as compared to CD below pressures of 27 GPa. This is only true for the nucleation mechanism found here; CD is favored via the concerted mechanism.8,22 Our ordering of barrier heights is consistent with the experimental observation that CCG can be synthesized at a lower temperature than CD near 17 GPa. On the other hand, our calculations also show that HD should form more readily than CCG at this pressure, which is not observed. One possible resolution is that the kinetically accessible phase of CCG has not been found yet; new forms are still being discovered.23 Our model also has some limitations: the lateral interfaces of the critical nuclei are not included, nor defects, which may play a role in the nucleation process.
While not all issues are resolved, our calculations show that the kinetics of the nucleation mechanisms are qualitatively different from the concerted mechanisms. The nucleation mechanisms are understood in terms of the structure and enthalpy of a series of stable nuclei (the intermediate minima on the MEP in Fig. 1). The enthalpy of these nuclei are described by two components; one from the interface energy between the phases, and the other from the enthalpy of the new phase. The former is positive with respect to graphite and the latter is negative. The ratio of these two components depends on the nucleus size; it reaches a maximum at the critical nucleus. Therefore, the enthalpy of the critical nucleus determines the barrier of the phase transition, for example C1 and H1 in Fig. 1. In order to compare competing phase transition barriers, it is important to find the critical nuclei. In Ref. 11, for example, the enthalpy of diamond nuclei in graphite are calculated, but the size of the nuclei is chosen arbitrarily and there is no guarantee that these will be close to the critical size.
Under increasingly high pressure conditions, the concerted mechanism will be faster than nucleation. To compare these two mechanisms, we can write the enthalpy barriers as ΔE + PΔV, where ΔE is the energy barrier, ΔV is the activation volume, and P is the pressure. For any phase transition from graphite, ΔV is negative in both mechanisms, as can be seen from the slopes in Fig. 4. In a system with a large number of atoms, the volume at the saddle point for the concerted mechanism is smaller than a nucleation mechanism. At increasing pressure, the PΔV term will dominate the enthalpy barrier and the concerted mechanism, with a more negative ΔV, will eventually be spontaneous.
It will be interesting to investigate the size at which the reactions favor a truly local nucleation mechanism, as was done with the G-SSNEB method for CdSe.12 In a large graphite simulation at pressures between 15 and 27 GPa, the critical nucleus must be limited in lateral extent, and it may involve more graphite layers. Our current calculations correspond to the limit of layer-by-layer propagation of the new phase in graphite.
In conclusion, we are able to calculate new nucleation mechanisms with a modest number of atoms over a wide pressure range using the G-SSNEB method. The barrier from graphite to three candidate CCG structures, bct C4, M, and Z-carbon, are calculated to be lower than to CD in the pressure range at which they are synthesized in experiment. The barriers to HD and to CD are in the same order as seen in experiments, demonstrating the usefulness of the G-SSNEB method in exploring solid-state potential energy surfaces.
The method development was supported by the National Science Foundation under CHE-1152342 and the calculations were supported as part of the program “Understanding Charge Separation and Transfer at Interfaces in Energy Materials (EFRC:CST)” an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0001091. We acknowledge the Texas Advanced Computing Center for computational resources, Qiang Zhu and Maximilian Amsler for insightful discussions, and Carol King for help with the article.