Molecular Dynamics Simulations of Shock Compressed Graphite

Subscriber access provided by Library, Special Collections and Museums, University of Aberdeen

Article

Molecular Dynamics Simulations of Shock Compressed Graphite Nicolas Pineau J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/jp403568m • Publication Date (Web): 28 May 2013 Downloaded from http://pubs.acs.org on June 5, 2013

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

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

Page 2 of 29

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

The Journal of Physical Chemistry

Molecular Dynamics Simulations of Shock Compressed Graphite Nicolas Pineau∗ CEA/DAM/DIF, F-91297 Arpajon, France E-mail: [email protected]

∗ To

whom correspondence should be addressed

1 ACS Paragon Plus Environment


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

Abstract We present molecular dynamic simulations of the shock compression of graphite with the LCBOPII potential. The range of shock intensities covers the full range of available experimental data, including near tera-pascal pressures. The results are in excellent agreement with the available DFT data and point to a graphite-diamond transition for shock pressures above 65 GPa, a value larger than the experimental data (20 to 50 GPa). The transition mechanism leads preferentially to hexagonal diamond through a diffusionless process but is submitted to irreversible regraphitization upon release: this result is in good agreement with the lack of highly ordered diamond observed in post-mortem experimental samples. Melting is found for shock pressures ranging from 200 to 300 GPa, close to the approximate LCBOPII diamond melting line. A good overall agreement is found between the calculated and experimental Hugoniot data up to 46 % compression rate.

Keyword: diamond, phase transition, diffusionless, LCBOPII

Introduction The behavior of bulk carbon submitted to high pressures is a field of wide interest both from the academic point of view due to the specific properties of carbon materials (versatile covalent bonding, phase changes, mechanical strength) and from the applied point of view for materials science or astrophysics (structural transitions, geophysics of planetary interior 1–3 , formation of ultrahard carbon materials in meteorites 4,5 ). The resulting number of potential applications has motivated numerous studies of the carbon phase diagram through isotropic and uniaxial compression experiments 6–18 or simulations 19–25 . For instance, the shock compression of graphite has been the subject of much experimental work yielding shock Hugoniot data for pressures up to the teraPascal regime 26–33 . Most of these studies focus on the moderate compression regime, where graphite undergoes a phase transition toward diamond (G/D transition), but very high compression experiments are also important in order to investigate the melted region of the phase diagram. The experimental data display some 2 ACS Paragon Plus Environment

Page 3 of 29

Page 4 of 29

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


variability in the observed G/D threshold pressure: the precursor work by Alder and Christian point to a transition triggered at 18 GPa and completed between 40 and 60 GPa 26 , while Gust finds a 25/35 GPa pressure range depending on the origin of the sample 27 . Later Yamada and Tanabe found the transition pressure can be lowered to below 15 GPa when using graphite samples containing voids 32 . From a mechanistic point of view, the shock-induced G/D transition is generally assumed to occur along two distinct paths: a presumably martensitic transition resulting from the 1D compression perpendicular to the basal planes seems to prevail for highly oriented graphite samples, while a diffusion-based mechanism implying a liquid intermediate appears to be more effective in highly defective or heterogeneous samples 29–33 . The variability of the measured transition pressures may relate to the amount of each mechanism involved during the experiments, the diffusion process supposingly yielding substantially lower resistance to diamond nucleation 26,27,30,33 . Overall, the experimental data suggest the Hugoniot curve and the resulting G/D transition properties depend substantially on the degree of crystallinity, purity and orientation of the shocked sample. Usually, the theoretical approaches are particularly suited to address such questions; unfortunately, simulations of the high pressure properties of graphite are scarce due to the lack of interaction models capable of describing with sufficient accuracy its subtle combination of covalent and non-bonded interactions. Only recently Mundy et al. used the multiscale shock technique (MSST) within the Density Functional Theory framework to follow the G/D transition induced by a 12 km/s shock wave, corresponding to a particle velocity and pressure of ∼ 5 km/s and ∼ 130 GPa respectively 34 . Their simulations are consistent with a diffusionless transition and lead to cubic diamond. However they were not able to yield diamond for lower shock velocities possibly due to the limited system sizes accessible to these calculations (a few hundred atoms for instance). This limitation prevents extensive ab initio studies of the shock properties of graphite over the full range of available experimental data: as a consequence, and to our knowledge, the full Hugoniot of graphite has never been obtained from atomistic simulations. Given the computational cost of the ab initio methods, reactive empirical potentials have



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

emerged as an efficient and reasonably accurate tool and are currently widely used for the simulation of shock wave processes 22–24 . In the case of carbon, several reactive potentials have been built of which the most popular are the EDIP potential 35 and the Tersoff-derived potentials REBO 36,37 , AIREBO 38 and LCBOP-like 39–42 . However most of those either lack the long-range component that is necessary to mimic the non-bonded interactions at play in graphite, or do not reproduce accurately the DFT energy barrier of the G/D transition upon isothermal compression; on the contrary, the LCBOPII potential 40–42 has been shown to accurately reproduce this barrier, making it attractive for studies involving large scale graphite or graphitic materials under standard or high pressure and temperature conditions 43 . In this study we present direct shock compression simulations of monocrystalline hexagonal (ABAB) graphite along the direction perpendicular to the basal planes, using the LCBOPII potential. We simulate a wide range of piston velocities (0.1 to 15 km/s) to match the available range of experimental shock pressures, with particular emphasis on the onset and mechanism of the G/D transition and on the stability of the diamond phase upon release.

Computational details The LCBOPII potential The LCBOPII potential is based on the approximations of Abell and Tersoff which imply that the total energy of a set of interacting particles can be split into pair contributions, each of these contributions being the sum of a repulsive and an attractive term. The repulsive term is purely pairwise while the attractive term is the product of a pairwise interaction and an environment-dependent contribution, the so-called “bond-order”: the environment of a given bond includes atoms up to third-neighbours through configurational contributions such as bond distances, bond angles, dihedrals and conjugation corrections. LCBOPII includes long-range interactions to account for non-bonded interactions and is therefore particularly well suited for graphite or graphitic materials, where π -repulsion plays a crucial role. It also contains a medium-range contribution fitted 4 ACS Paragon Plus Environment

Page 5 of 29

Page 6 of 29

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


to ab-initio calculations of the dissociation energy curves of simple, double and triple bonds, improving the reactive properties which is particularly suited for the description of liquid carbon and phase transitions. LCBOPII has been successfully applied to a variety of carbon systems including graphene, graphite, liquid carbon or nanometer-sized clusters 42–46 . An important feature of the potential for the present study is its ability to describe the barrier height of the G/D transition accurately 43 , which gives confidence in the correct description of the kinetics of this transformation through shock compression.

Methodology In order to simulate shock-compression experiments, we build a quasi 1D sample which longitudinal size is substantially larger than the transverse sizes. Periodic boundary conditions are applied in the transverse directions (y, z) only. The shock is initiated by a purely repulsive potential wall with infinite mass, the so-called piston, projected at a constant velocity u p on the sample along the x axis. All the simulations are run in the microcanonical (NVE) ensemble with a time step of 0.2 fs. The parameters of the graphite unit cell and the relative atomic positions correspond to the minimized 0 K structure, yielding an initial bulk density of 2.261 g/cm3 . The resulting pressure at 300 K (0.3 GPa) is negligible with respect to the pressures achieved through shock-compression. All the atoms are assigned initial velocities using a gaussian distribution corresponding to T = 600 K and thermal equilibration is performed before the impact of the piston, yielding a temperature close to 300 K. All the simulations have been carried out using the molecular dynamics code STAMP developed at CEA and the TERA100 supercomputer.



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

Results and discussion We ran a first set of simulations in the nonreactive shock-compression regime (no phase transition occurring), then a second one in the reactive regime (G/D transition and melting). The sample size in the transverse directions y and z was set to 3.4 × 3.4 nm: unless specified, this value was checked to eliminate any finite size effect. In order to allow the formation of a stationary shockwave in each sample while saving simulation time, two sample lengths were adopted along the shock propagation axis x: 31 nm for the non-reactive cases and 101 nm for the reactive cases to account for the resulting longer induction time. The samples contained 40320 and 134400 atoms respectively corresponding to 270 and 900 CPUs. Piston velocities u p ranging from 0.1 to 2 km/s were applied in the nonreactive cases, and from 2 to 15 km/s in the reactive case. The simulations were continued until the shockwave front reached the end of the sample, checking carefully that stationary waves had been obtained. The shock compression process yielded three distinct regimes: an elastic regime for u p ≤ 2 km/s, a G/D transition regime between 3 and 6 km/s, and a melting regime above 6 km/s. The results are synthesized in the following series of figures. In Figures 1 and 2 we present the profiles and structures of a selection of stationary shockwaves. The shock velocities Us and pressures P are plotted as a function of piston velocity u p in Figure 3, yielding the corresponding Hugoniot in the

ρ − P plane and P − T phase diagram (Figures 4 and 5). In the next subsections, we analyze these results for each regime specifically.


Page 7 of 29

Page 8 of 29

15

(a)

main plot (a) 15 km/s (b) 12 km/s (c) 8 km/s (d) 6 km/s (e) 5 km/s (f) 4 km/s (g) 2 km/s

(b) 2

10 Up (km/s)

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


1.5

(c)

inset

1 (j)

(d)

5

(h) (i)

(h) 2.00 km/s (i) 1.50 km/s (j) 1.00 km/s (k) 0.75 km/s (l) 0.50 km/s (m) 0.30 km/s

(k) (l)

0.5

(m)

(e) (f)

-150 -100

-50

0

(g)

0 -200

-150

-100 0 100 -50 50 Distance to shock front (angstrom)

150

200

Figure 1: Shock profiles for a selection of simulations with piston velocities from 0.3 to 15 km/s. Most cases display a one-wave pattern typical of quasi-elastic compression for low shock intensities, or overdriven structural transformations for high shock intensities. Only the 3 to 5 km/s cases, corresponding to the onset of the G/D transition, display a finite thickness “reaction zone”. The main plot focuses on the reactive cases (u p ≥ 2 km/s) and the inset focuses on the elastic and pseudo-elastic regimes (u p ≤ 2 km/s).

2 km/s 3 km/s 4 km/s 5 km/s 6 km/s 7 km/s 8 km/s

Figure 2: Snapshots of the shocked graphite sample for a selection of piston velocities. Up to 2 km/s, no or little structural transformation is observed. At 3 km/s the formation of diamond nuclei is observed locally but remains a slow process with respect to the achievable simulation times. At 4 and 5 km/s the formation of bulk diamond is achieved with a number of defects such as stacking faults or localized amorphous domains. At 8 km/s and above, full melting of the sample is achieved. Color code: 2-fold carbon (chain-like) in green, 3-fold (graphite-like) in blue, 4-fold (diamond-like) in red, 5-fold and beyond in yellow.



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

Figure 3: Shock velocity Us and longitudinal pressure P vs. piston velocity u p and comparison to the available experimental and DFT data. For the LCBOPII results, the color code is: orange circles for elastic, blue squares for quasi-elastic, green triangles (up) for G/D equilibrium, red triangles (down) for diamond, maroon star for diamond/liquid equilibrium, magenta crosses for liquid. The full black curve is a guide to the eye. The experimental data from Ref. 27 are given in black (dashed, dotted-dashed and dotted lines for pyrolytic, Ceylon and synthetic graphite, double-dotted-dashed line for vitrous carbon), and from Refs. 28,31 in light and dark green respectively. The black cross represents the MSST-DFT simulations from Ref. 34 . For pressure, the orange full line is calculated from the second Hugoniot relation: the good agreement between the measured (colour symbols) and calculated pressures confirms that stationarity was achieved in all simulations. The orange dashed line is an extrapolation for the purely elastic shocks calculated from the simulations with u p < 0.4 km/s.


Page 9 of 29

Page 10 of 29

Figure 4: Hugoniot curve and comparison to experimental and DFT data. The color code is identical to Figure 3. The inset focuses on the elastic and pseudo-elastic regime.

60000 900

50000 600

Temperature (K)

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


40000 300

30000 0

10

20

30

40

0 50

20000

10000

0

0

200

400

600

800

1000

Pressure (GPa)

Figure 5: Locus of the Hugoniot data in the P-T phase diagram of carbon. The color code is identical to Figure 3. The orange dashed lines are the estimated graphite/diamond, graphite/liquid and diamond/liquid coexistence lines. The inset is a focus on the shock data for graphite and G/D transition.



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

Elastic compression Elastic and pseudo-elastic compressions are found for piston velocities up to 2 km/s: they result in the preservation of the typical layered structure of graphite and agree nicely with the experimental Hugoniot data (inset in Figure 4). For piston velocities lower than 1 km/s the shock profiles display a single wave pattern characteristic of purely elastic shockwaves: the compression results mainly in the decrease of the interlayer distance between successive graphene sheets. Accordingly the Us (u p ) and P(u p ) curves yield typical linear and quadratic behaviors respectively (Figure 3). Plotting the corresponding Hugoniot data in the (P,T) phase diagram (orange symbols in Figure 5) shows that this regime matches closely the stability region of graphite, that is below the G/D coexistence line (inset in Figure 5). Above 1 km/s, the shock front tends to be less sharp due to the onset of intralayer out-of-plane oscillations which amplitude increases with shock strength. These oscillations are accompanied by non linear/quadratic behaviors of Us and P with progressive drifts toward values lower than their elastic approximation (defined as the extrapolation of the data for u p < 0.4 km/s). The location of the Hugoniot data in the LCBOPII phase diagram indicates that the shocked sample is now beyond the G/D coexistence line (blue squares), that is inside the metastable region of graphite: the onset of the out-of-plane oscillations is associated with a large increase of the Hugoniot slope dT /dP, resulting in a large increase in temperature at the expense of pressure. However these shock intensities are too weak to initiate the G/D transition, as shown by the absence of the interlayer bonds that would be necessary to nucleate sp3 clusters.


Page 11 of 29

Page 12 of 29

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


Figure 6: Influence of the section size on the local diamond structures obtained with the u p =4 km/s shocks. Top: 3.4×3.4 nm section. Bottom: 6.8×6.8 nm section. The highlighted regions (red) correspond to hexagonal diamond structures, showing increasing fractions of this phase in larger samples. The identification of the cubic and hexagonal phases are described in Figure 8.

Graphite/diamond transition G/D mechanism The graphite/diamond transition was obtained for piston velocities ranging between 3 and 6 km/s. Between 3 and 5 km/s the profiles exhibit a two-wave pattern consisting of a precursor elastic wave and a follower transformation wave. At u p = 3 km/s, we observe the formation of the first interlayer bonds as well as the first sp3 nuclei, which size is sufficient to form stable diamond seeds. Above 3 km/s, the input energy is sufficient to trigger the transition across the whole sample and the thickness of the transformation zone progressively decreases, yielding an overdriven process between 5 and 6 km/s. The detailed mechanism of the transition is presented in Figure 7(top) through the time-evolution of a three-layer slice of graphite shock-compressed at u p = 4 km/s. The transition occurs along a diffusionless mechanism 47 and leads to a cubic diamond structure through a short-lived intermediate similar to the “layered diamond” evidenced by Mundy et al. in



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

0 ps

0 ps

0.1 ps

0.2 ps

0.2 ps

0.3 ps

0.4 ps

0.5 ps

0.8 ps 1.0 ps 1.2 ps 1.4 ps

Page 13 of 29

0.9 ps

2 ps

1.1 ps

3 ps

Figure 7: Mechanism of the G/D transition seen through the time-evolution of a three-layer slice of shock-compressed graphite for u p = 4 km/s. Top: 3.4×3.4 nm section. Bottom: 6.8×6.8 nm section. The small sample leads to the thermodynamically more stable cubic diamond phase. The large sample leads to the epitaxially more favorable hexagonal diamond (lonsdaleite) resulting from direct formation of interplanar bonds through “instantaneous” 1D compression of ABAB graphite. The identification of the cubic/hexagonal phases is described in Figure 8.


Page 14 of 29

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


Cubic diamond

Hexagonal diamond

Figure 8: 2D-schematic structures of hexagonal and cubic diamond with the orientation corresponding to the one obtained in our shocked samples. This projection direction allows the identification of both phases based on their apparent local symmetry : planar symmetry perpendicular to the shock direction for the hexagonal phase (grey dashed lines), and central symmetry of the “unit cell” for the cubic phase (grey circles). All our structural analyses are based on this symmetry criterion (we checked carefully that no other orientation of one phase can yield the observed pattern of the other, which overrules the possibility of ambiguous assignment). Hexagonal diamond

ABAB graphite

intermediate

Hexagonal diamond

Cubic diamond

ABAB graphite

ABCABC graphite

Cubic diamond

Hexagonal diamond: Khaliullin's mechanism

ABAB graphite

intermediate

Hexagonal diamond

Figure 9: Proposed mechanisms for the G/D transitions observed in our simulations. The two mechanisms require lateral shifts of the graphite sheets followed by concerted atomic displacements through a diffusionless process. Note that the formation of cubic diamond involves larger lateral shifts, probably resulting in larger friction as illustrated by the thickness of the grey dashed lines located in the friction regions. We provide comparison with the mechanism used by Khaliullin et al. to generate hexagonal diamond nuclei in ABAB graphite: in this mechanism the shift of the graphite sheets is similar, then the atomic displacements differ (red arrows vs. black arrows) resulting in an alternate orientation of the hexagonal diamond phase 25 .



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

their MSST-DFT simulations 34 : this transient structure is characterised by the completion of the formation of the interlayer bonds, while the 6-membered graphite rings have not transited to their final armchair configuration (typically at t = 0.5 ps and 1.2 ps for the small and large samples in Figure 7). The timescale for the full transition process is of the order of a picosecond, one order of magnitude slower than that observed by Mundy which is consistent with the lower shock intensity in our case (4 km/s vs. 5 km/s). Careful observation of Figure 2 shows the inclusion of hexagonal diamond layers within the cubic diamond matrix, with the resulting stacking faults oriented perpendicular to the shock propagation. The predominance of the cubic phase and the homogeneity of the transition in the transverse directions suggest possible size effects resulting from the periodic boundary conditions along y and z 48 . Therefore we performed an extra simulation at 4 km/s with doubled transverse size (6.8×6.8 nm, 537600 atoms, 3600 CPUs). The shock velocity and all the thermodynamic properties computed in this simulation are within 1% of the values computed with the small sample, overruling any size-effect on macroscopic shock properties. However from a structural point of view we observe the onset of plastic deformations which are not found with the smaller sample, as well as a predominant hexagonal diamond phase (bottom of Figures 6 and 7). This result suggests that the transverse periodic boundary condition on the small sample favors a rearrangement of the initial ABAB stacking into the cubic ABCABC stacking, leading to the artificial formation of the thermodynamically more stable cubic diamond. In Figure 9, we propose schematic transition mechanisms involving possible rearrangements leading to both forms of diamond: the orientation observed in our samples (Figure 6) shows that the relevant mechanism differs slightly from the one used by Khaliullin et al. to generate hexagonal diamond nucleation in a graphite matrix 49 . The timescale is comparable for both sample sizes in spite of large variations (up to several picoseconds) in the formation time of the first interlayer bonds after compression: these discrepencies have also been observed between distinct slices of a same sample, therefore they are not significant. Beyond this induction time, the structural transition ranges between 0.5 and 1 picosecond in both cases.


Page 15 of 29

Page 16 of 29

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


Pressure threshold The G/D transformation induces a strong inflexion of the Us (u p ) and P(u p ) curves with a nearly constant shock velocity in the 3-4 km/s range (Us ∼ 10 km/s) and longitudinal pressures substantially lower than the elastic extrapolation (Figure 3). A linear Us (u p ) behavior is regained for the overdriven transition (4 km/s < u p < 6 km/s). We calculated the Us (u p ) slopes in the linear regimes corresponding to elastic graphite and G/D transition, yielding 3.72 and 1.86 respectively, somewhat larger than the experimental values (2.17 and 1.50 for Ceylan graphite 27 ): accordingly the calculated Hugoniot pressures are substantially overestimated with respect to experiments (Figure 4). Meanwhile the onset pressure for diamond nucleation is found close to 65 GPa, when experimental data point toward values between 20 GPa and 50 GPa. These overestimations are usual with MD simulations and do not necessarily result from a flaw of the LCBOPII potential, for which excellent agreement is obtained with the DFT simulations of Mundy at u p =5 km/s. They possibly result from the simulation conditions through: i) the transverse size effects due to periodic boundary conditions, ii) the non-stationarity of shockwaves due to insufficient length in the longitudinal direction, iii) the monocrystalline nature of the samples which lack the defects necessary to help nucleate transitions, or iv) the variable density of the polycrystalline experimental samples which has a strong impact on shock behavior, larger densities leading to larger shock velocities and pressures (experimental densities range from 1.54 to 2.25 g/cm3 depending on the origin or preparation conditions 27 , while the density of our initial monocrystalline sample is 2.26 g/cm3 ). Points i) and ii) have been overruled by carefully checking the convergence of the shock properties with sample dimensions. Points iii) and iv) are out of the scope of the present paper due to the potentially large variety of defects to investigate and to the sample size issues related to polycrystalline samples.

Stability of diamond upon release In order to check the stability of the diamond structure, we let the 4 km/s simulation continue after the shock has reached the end of the sample, allowing substantial release of the structure. 15 ACS Paragon Plus Environment


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

The release waves propagate through the shocked sample and lead to sample/piston separation, resulting in a sample without external interaction. In Figures 10 and 11 we present the temporal evolutions of the sample structure as well as the thickness and pressure of the diamond phase. The propagation of the release waves back and forth in the sample result in local oscillations of the measured pressure within a slowly decaying [-20 GPa;+20 GPa] range which, combined with the high shock temperatures (>3500 K), lead to regraphitization at the G/D interfaces. The regraphitization kinetics is estimated by linearly fitting the time dependence of the diamond thickness beyond the diamond peak (Figure 11, top-right, blue dotted line). We find a velocity close to 9 Å/ps, predicting a life-expectancy of ∼50 ps for the diamond cluster formed in this simulation. This result indicates that the G/D transition through diffusionless processes is probably fully reversible, leading to predominantly short-lived diamond phases due to their high structural ordering. We stress that the regraphitization velocity should not depend on the diamond thickness, which is several orders of magnitude larger in experiments (rather on the micrometer scale): we expect this velocity should rather depend on shock intensity through the temperature and the amplitude of the pressure fluctuations achieved in the shock and subsequent release waves.


Page 17 of 29

Compression

Release

up = 4 km/s

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


t

Page 18 of 29

Repulsive wall Figure 10: Structural evolution of graphite undergoing a u p = 4 km/s shock followed by release. The color code is indicative of the local structure (blue for sp2 /graphite, red for sp3 /diamond).



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

Figure 11: Time evolution of the graphite/diamond fractions (top-left), diamond thickness (topright) and pressure profile (bottom) upon release. In the latter, each colored curve corresponds to the pressure profile in the remaining diamond phase at the indicated simulation time.


Page 19 of 29

Page 20 of 29

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


Discussion These simulations give a clear picture for one of the assumed G/D mechanisms, that is the diffusionless transition upon shock compression of highly oriented graphite. The direct transition leads to highly oriented diamond clusters submitted to high temperatures, conditions which strongly favor the reverse regraphitization process upon release. This analysis is supported by the lack of experimental observations of highly oriented diamonds in post-mortem samples, when the timeresolved free surface velocity data point to a high-yield diffusionless transition (70-80%) for highly oriented graphite 5,30,50 . The substantial overestimation of the transition pressure may also have a physical explanation: the diffusional and diffusionless G/D mechanisms may yield rather distinct pressure thresholds, with a non-negligible temperature dependance for the diffusional case. Several experimental studies point to a larger pressure threshold for the diffusionless case 26,27,30,33 while others suggest otherwise 50 . Due to our simulation conditions, only the diffusionless mechanism is expected to occur thus yielding a threshold pressure corresponding to this unique process; comparatively, experiments always yield a balance between the two processes, depending on the nature of the sample, preparation conditions and experimental setup. Therefore the discrepancy between our simulations and the experiments could arise primarily from the difference in the “experimental” setups, rather than from the quality of the interaction potential: the good agreement with DFT for the G/D transition is also indicative of a reasonable LCBOPII behavior in this regime. In order to test this assumption it would be necessary to simulate the diffusional mechanism, a process that occurs most likely in highly defective regions of the material, possibly at grain boundaries (experiments point to a strong dependence of the G/D pressure and transition yield on the size-distribution of the crystallites and on temperature 29,32,33 ). Such simulations require sample sizes which are difficult to reach through conventional molecular dynamics with the current computational capabilities. Metadynamics could provide an efficient alternate approach to address this problem by allowing to estimate and compare the free energy profiles of the relevant mechanisms, although under a thermodynamic equilibrium approximation. Specific collective variables such as 19 ACS Paragon Plus Environment


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

appropriate local order parameters may be designed for that purpose 51 . However metadynamics are not available in the STAMP molecular dynamics code and such study is beyond the scope of the present paper.

Graphite/liquid transition The graphite/liquid transition is triggered around piston velocities of 7 km/s as shown by the slight inflexion in the Us (u p ) and P(u p ) curves. Accordingly the Hugoniot data follow the diamond melting line within an approximate [200 GPa;300 GPa] pressure range, confirming the prediction of Nellis et al. 31 . At these large shock intensities, a single wave pattern is observed showing that the G/L transition tends to be overdriven: the energy brought by the piston is large enough to trigger the phase transformation with negligible induction time as shown by the narrow shockcompressed graphite region. Accordingly a linear Us (u p ) relationship is regained for u p ≥ 8 km/s, as the Hugoniot curve enters the liquid region of the phase diagram. The calculated slope of 1.64 is somewhat lower than that for compressed diamond and again larger than the experimental value of 1.46 from Nellis. With bond order potentials such as LCBOPII, pressure overestimations at high compression rates are usual and arise from two sources of error. First, the insertion of second neighbour atoms in the sphere of interaction of the first neighbours induce unphysically large steric repulsion because the potential was not originally designed to address such compression rates 24 . Second, in this regime liquid carbon was shown to display a metallic behavior 16 and the role of thermal electronic excitations, not included in the LCBOPII formalism, cannot be neglected. Such effects were quantified by Romero et al. for shocked diamond, showing a substantial decrease of the shock temperature for DFT simulations with finite electronic temperature 52 : the electronic heat capacity acts as an energy sink for the system, resulting in the decrease of shock pressure and temperature at a given shock intensity. Considering these limitations, the overestimation of the shock pressure in the liquid regime seems reasonable. An appropriate treatment of the steric repulsion at high compression would certainly reduce this discrepancy, although only partially; full correction could 20 ACS Paragon Plus Environment

Page 21 of 29

Page 22 of 29

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


only be achieved by addressing the finite electronic temperature problem, a very difficult task that is usually beyond the scope of such empirical potentials.

Conclusion Our shock compression simulations of graphite yield good agreement with the available experimental and DFT data in the elastic and graphite/diamond regimes. The G/D transition is obtained for piston velocities ranging from 3 to 6 km/s via a purely diffusionless process leading to hexagonal diamond for the largest simulated samples. We find the G/D transition pressure is overestimated which we believe to result partly from the undefected/monocrystalline nature of the simulated sample, from the resulting purely diffusionless nature of the G/D transition, and possibly from an excessive stiffness of the potential. Unfortunately, the system sizes necessary to simulate the alternate diffusional process exceeds the current capabilities of our simulations. We predict the G/D transition is reversible under the high temperatures achieved, leading to the fast regraphitization of the highly-ordered diamond phase. The melting of graphite is observed above a [200 GPa ; 300 GPa] pressure range in good agreement with the Nellis estimation, a value that still needs experimental validation. The simulated shock pressures in the tera-Pascal regime are substantially overestimated, but we consider the agreement to be reasonable since LCBOPII was not designed to deal with such high compression ranges nor the resulting electronic conductivity. Only partial improvement of the potential for high-pressure properties seems achievable in that respect.

Acknowledgement N. P. acknowledges gratefully Dr. J. H. Los, Dr. L. Soulard, Dr. C. Denoual, Dr. J.-B. Maillet and Pr. Y. Bai for fruitful discussions.



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

References (1) Hubbard, W. B. Interiors of Giant Planets. Science 1981, 214, 145–149. (2) Seager, S.; Kuchner, M.; Hier-Majumder, C. A.; Militzer, B. Mass-Radius Relationships for Solid Exoplanets. The Astrophysical Journal 2007, 669, 1279–1297. (3) Madhusudhan, N.; Lee, K. K. M.; Mousis, O. A Possible Carbon-Rich Interior in Super-Earth 55 Cancri e. The Astrophysical Journal Letters 2012, 759, L40. (4) Kenkmann, T.; Hornemann, U.; Stöffler, D. Experimental Shock Synthesis of Diamonds in a Graphite Gneiss. Meteoritics and Planetary Science 2005, 40, 1299–1310. (5) Ferroir, T.; Dubrovinsky, L.; El Goresy, A.; Simionovici, A.; Nakamura, T.; Gillet, P. Carbon Polymorphism in Shocked Meteorites: Evidence for New Natural Ultrahard Phases. Earth and Planetary Science Letters 2010, 290, 150–154. (6) Bundy, F. P. Melting of Graphite at Very High Pressure. J. Chem. Phys. 1963, 38, 618–630. (7) Bundy, F. P. Direct Conversion of Graphite to Diamond in Static-Pressure Apparatus. J. Chem. Phys. 1963, 38, 631–643. (8) Pavlovskii, M. N. Shock Compression of Diamond. Sov. Phys. Solid State 1971, 13, 741. (9) Bradley, D. K.; Eggert, J. H.; Hicks, D. G.; Celliers, P. M.; Moon, S. J.; Cauble, R. C.; Collins, G. W. Shock Compressing Diamond to a Conducting Fluid. Phys. Rev. Lett. 2004, 93, 195506. (10) Nagao, H. et al. Hugoniot Measurement of Diamond Under Laser Shock Compression up to 2 TPa. Phys. Plasmas 2006, 13, 052705. (11) Brygoo, S.; Henry, E.; Loubeyre, P.; Eggert, J.; Koenig, M.; Loupias, B.; BenuzziMounaix, A.; Le Gloahec, M. R. Laser-Shock Compression of Diamond and Evidence of a Negative-Slope Melting Curve. Nature Materials 2007, 6, 274–277. 22 ACS Paragon Plus Environment

Page 23 of 29

Page 24 of 29

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


(12) Hicks, D. G.; Boehly, T. R.; Celliers, P. M.; Bradley, D. K.; Eggert, J. H.; McWilliams, R. S.; Jeanloz, R.; Collins, G. W. High-Precision Measurements of the Diamond Hugoniot in and Above the Melt Region. Phys. Rev. B 2008, 78, 174102. (13) Dewaele, A.; Datchi, F.; Loubeyre, P.; Mezouar, M. High Pressure-High Temperature Equations of State of Neon and Diamond. Phys. Rev. B 2008, 77, 094106. (14) Knudson, M. D.; Desjarlais, M. P.; Dolan, D. H. Shock-Wave Exploration of the HighPressure Phases of Carbon. Science 2008, 322, 1822–1825. (15) Bradley, D. K.; Eggert, J. H.; Smith, R. F.; Prisbrey, S. T.; Hicks, D. G.; Braun, D. G.; Biener, J.; Hamza, A. V.; Rudd, R. E.; Collins, G. W. Diamond at 800 GPa. Phys. Rev. Lett. 2009, 102, 075503. (16) Eggert, J. H.; Hicks, D. G.; Celliers, P. M.; Bradley, D. K.; McWilliams, R. S.; Jeanloz, R.; Miller, J. E.; Boehly, T. R.; Collins, G. W. Melting Temperature of Diamond at Ultrahigh Pressure. Nature Physics 2010, 6, 40–43. (17) McWilliams, R. S.; Eggert, J. H.; Hicks, D. G.; Bradley, D. K.; Celliers, P. M.; Spaulding, D. K.; Boehly, T. R.; Collins, G. W.; Jeanloz, R. Strength Effects in Diamond Under Shock Compression from 0.1 to 1 TPa. Phys. Rev. B 2010, 81, 014111. (18) Lang, J. M., Jr.; Gupta, Y. M. Strength and Elastic Deformation of Natural and Synthetic Diamond Crystals Shock Compressed along [100]. J. Appl. Phys. 2010, 107, 113538. (19) Blase, X.; Gillet, P.; San Miguel, A.; Mélinon, P. Exceptional Ideal Strength of Carbon Clathrates. Phys. Rev. Lett. 2004, 92, 215505. (20) Wen, B.; Zhao, J.; Bucknum, M. J.; Yao, P.; Li, T. First-Principles Studies of Diamond Polytypes. Diam. Rel. Mat. 2008, 17, 356–364. (21) Goldman, N.; Fried, L. E. Extending the Density Functional Tight-Binding Method to Carbon Under Extreme Conditions. J. Phys. Chem. C 2012, 116, 2198–2204. 23 ACS Paragon Plus Environment


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

(22) Zybin, S. V.; Elert, M. L.; White, C. T. Orientation Dependence of Shock-Induced Chemistry in Diamond. Phys. Rev. B 2002, 66, 220102. (23) Marks, N. A.; Lattemann, M.; McKenzie, D. R. Nonequilibrium Route to Nanodiamond with Astrophysical Implications. Phys. Rev. Lett. 2012, 108, 075503. (24) Perriot, R.; Lin, Y.; Zhakhovsky, V. V.; Pineau, N.; Los, J. H.; Maillet, J.-B.; Soulard, L.; White, C. T.; Oleynik, I. I. Shock Compression of Diamond: Molecular Dynamics Simulations using Different Interatomic Potentials. 2012. (25) Khaliullin, R. Z.; Eshet, H.; Kuhne, T. D.; Behler, J.; Parrinello, M. Nucleation Mechanism for the Direct Graphite-to-Diamond Phase Transition. Nature Materials 2011, 10, 693–697. (26) Alder, B. J.; Christian, R. H. Behavior of Strongly Shocked Carbon. Phys. Rev. Lett. 1961, 7, 367–369. (27) Gust, W. H. Phase Transition and Shock Compression Parameters to 120 GPa for 3 Types of Graphite and for Amorphous Carbon. Phys. Rev. B 1980, 22, 4744–4756. (28) Erskine, D. J.; Nellis, W. J. Shock-Induced Martensitic Phase-Transformation of Oriented Graphite to Diamond. Nature 1991, 349, 317–319. (29) Hirai, H.; Kukino, S.; Kondo, K. Predominant Parameters in the Shock Induced Transition from Graphite to Diamond. J. Appl. Phys. 1995, 78, 3052–3059. (30) Zhuk, A. Z.; Borodina, T. I.; Fortov, V. E.; Lash, A. A.; Val’Iano, G. E. Shock Metamorphism of the Graphite Quasimonocrystal. High Pressure Research 1997, 15, 245–254. (31) Nellis, W. J.; Mitchell, A. C.; McMahan, A. K. Carbon at Pressures in the Range 0.1-1 TPa (10 Mbar). J. Appl. Phys. 2001, 90, 696–698. (32) Yamada, K.; Tanabe, Y. Shock Induced Phase Transition of Oriented Pyrolytic Graphite to Diamond at Pressures up to 15 GPa. Carbon 2002, 40, 261–269. 24 ACS Paragon Plus Environment

Page 25 of 29

Page 26 of 29

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


(33) Kanel’, G. I.; Bezruchko, G. S.; Savinykh, A. S.; Razorenov, S. V.; Milyavskii, V. V.; Khishchenko, K. V. Submicrosecond Polymorphic Transformations Accompanying Shock Compression of Graphite. High Temperature 2010, 48, 806–814. (34) Mundy, C. J.; Curioni, A.; Goldman, N.; Kuo, I.-F. W.; Reed, E. J.; Fried, L. E.; Ianuzzi, M. Ultrafast Transformation of Graphite to Diamond: an Ab Initio Study of Graphite Under Shock Compression. J. Chem. Phys. 2008, 128, 184701. (35) Marks, N. A. Generalizing the Environment-Dependent Interaction Potential for Carbon. Phys. Rev. B 2000, 63, 035401. (36) Brenner, D. W. Empirical Potential for Hydrocarbons for Use in Simulating the Chemical Vapor-Deposition of Diamond Films. Phys. Rev. B 1990, 42, 9458–9471. (37) Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. A Second-Generation Reactive Empirical Bond Order (REBO) Potential Energy Expression for Hydrocarbons. J. Phys. Cond. Mat. 2002, 14, 783–802. (38) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. A Reactive Potential for Hydrocarbons with Intermolecular Interactions. J. Chem. Phys. 2000, 112, 6472–6486. (39) Los, J. H.; Fasolino, A. Intrinsic Long-Range Bond-Order Potential for Carbon: Performance in Monte Carlo Simulations of Graphitization. Phys. Rev. B 2003, 68, 024107. (40) Los, J. H.; Ghiringhelli, L. M.; Meijer, E. J.; Fasolino, A. Improved Long-Range Reactive Bond-Order Potential for Carbon. I. Construction. Phys. Rev. B 2005, 72, 214102. (41) Los, J. H.; Ghiringhelli, L. M.; Meijer, E. J.; Fasolino, A. Improved Long-Range Reactive Bond-Order Potential for Carbon. I. Construction (vol 72, pg 214102, 2005). Phys. Rev. B 2006, 73, 229901(E). (42) Ghiringhelli, L. M.; Los, J. H.; Meijer, E. J.; Fasolino, A.; Frenkel, D. Liquid Carbon: Structure Near the freezing line. J. Phys. Cond. Mat. 2005, 17, S3619–S3624. 25 ACS Paragon Plus Environment


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

(43) Los, J. H.; Pineau, N.; Chevrot, G.; Vignoles, G.; Leyssale, J.-M. Formation of Multiwall Fullerenes from Nanodiamonds Studied by Atomistic Simulations. Phys. Rev. B 2009, 80, 155420. (44) Ghiringhelli, L. M.; Los, J. H.; Meijer, E. J.; Fasolino, A.; Frenkel, D. Modeling the Phase diagram of Carbon. Phys. Rev. Lett. 2005, 94, 145701. (45) Fasolino, A.; Los, J. H.; Katsnelson, M. I. Intrinsic Ripples in Graphene. Nature Materials 2007, 6, 858–861. (46) Pineau, N.; Soulard, L.; Los, J. H.; Fasolino, A. Theoretical Study of the Nucleation/Growth Process of Carbon Clusters Under Pressure. J. Chem. Phys. 2008, 129, 024708. (47) The Transition should not be considered “martensitic” since the nature of Bonding changes substantially between two Carbon atoms belonging to neighboring Graphite sheets (nonBonded to covalent Bonding): martensitic Transitions are characterized by moderate Diffusionless Structural changes which preserve the nature of the Bonding, and the coordination numbers accordingly. (48) ABAB Graphite, the most abundent form of natural Graphite, leads to hexagonal Diamond (lonsdaleite) Under 1D-Compression, while ABCABC Graphite gives the thermodynamically more stable cubic Diamond. The Formation of cubic Diamond from ABAB Graphite more likely results from slow Processes, allowing for the relevant Structural reorganisation. (49) Khaliullin, R. Z.; Eshet, H.; Kuhne, T. D.; Behler, J.; Parrinello, M. Graphite-Diamond Phase coexistence Study Employing a Neural-Network Mapping of the Ab Initio Potential Energy Surface. Physical Review B 2010, 81, 100103(R). (50) Erskine, D. J.; Nellis, W. J. Shock-Induced Martensitic-Transformation of Highly Oriented Graphite to Diamond. J. Appl. Phys. 1992, 71, 4882–4886.


Page 27 of 29

Page 28 of 29

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


(51) Martoˇnák, R. Atomistic Simulations of Pressure-Induced Structural Transformations in Solids. Eur. Phys. J. B 2011, 79, 241–252. (52) Romero, N. A.; Mattson, W. D. Density-Functional Calculation of the Shock Hugoniot for Diamond. Phys. Rev. B 2007, 76, 214113.



initial sample velocity

release phase

graphite

diamond

shock front

Figure 12: Graphic for ToC.


time

compression phase

repulsive wall

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

Page 29 of 29

Molecular Dynamics Simulations of Shock Compressed Graphite

Recommend Documents