Simulation of the Elastic and Ultimate Tensile Properties of Diamond

Article pubs.acs.org/JPCA

Simulation of the Elastic and Ultimate Tensile Properties of Diamond, Graphene, Carbon Nanotubes, and Amorphous Carbon Using a Revised ReaxFF Parametrization Benjamin D. Jensen,†,‡ Kristopher E. Wise,† and Gregory M. Odegard*,‡ †

Advanced Materials and Processing Branch, NASA Langley Research Center, Hampton, Virginia 23681, United States Department of Mechanical Engineering - Engineering Mechanics, Michigan Technological University, Houghton, Michigan 49931, United States

Downloaded by UNIV OF CALIFORNIA SANTA BARBARA on September 5, 2015 | http://pubs.acs.org Publication Date (Web): September 2, 2015 | doi: 10.1021/acs.jpca.5b05889

‡

S Supporting Information *

ABSTRACT: In light of the enduring interest in using nanostructured carbon materials as reinforcing elements in composite materials, there is a significant need for a reliable computational tool capable to predict the mechanical properties, both elastic properties and properties at the point of fracture, in large-scale atomistic simulations. A revised version of the ReaxFF reactive force field parametrization for carbon, ReaxFFC‑2013, was recently published and is notable because of the inclusion of density functional theory (DFT)-derived mechanical data for diamond and graphite in the fitting set. The purpose of the present work is to assess the accuracy of this new force field for predicting the mechanical properties for several allotropes of carbon, both in the elastic regime and during fracture. The initial discussion focuses on the performance of ReaxFFC‑2013 for diamond and graphene, the two carbon forms for which mechanical properties were included in the parametrization data set. After it is established that simulations conducted with the new force field yield results that agree well with DFT and experimental data for most properties of interest, its transferability to amorphous carbon and carbon nanotubes is explored. ReaxFFC‑2013 is found to produce results that, for the most part, compare favorably with available experimental data for single and multiwalled nanotubes and for amorphous carbon models prepared over a range of densities. Although there is opportunity for improvement in some predicted properties, the ReaxFFC‑2013 parametrization is shown to generally perform well for each form of carbon and to compare favorably with DFT and experimental data. bond scission and formation, in MD simulations.10 This capability enables the simulation of complex mechanical phenomena at a small multiple of the cost of a traditional MD simulation, rather than the orders of magnitude greater expense associated with QM calculations. ReaxFF simulations have been shown to accurately describe a number of reactive processes, including the catalytic formation of carbon nanostructures,11−13 the oxidation of hydrocarbons,14 the formation of silicon carbide,15 the thermal decomposition of high-energy materials,16,17 and the thermal decomposition of polydimethylsiloxane.18 All of the potential energy terms involving bonded interactions in ReaxFF depend on the bond order function, a continuous function of the interatomic distance that goes to zero beyond covalent bond distances. Dispersion and electrostatic nonbonded interactions complete the ReaxFF. As with other MD force fields, the accuracy of the results produced with ReaxFF depends on the quality of its parametrization,

1. INTRODUCTION Because of the challenges associated with making direct measurements of the mechanical properties of carbon nanostructures, modeling and simulation have been extremely useful in advancing the understanding of the behavior of these materials under various loading conditions. Classical molecular dynamics (MD) simulations using fixed valence force fields, for example, have played an invaluable role in exploring the small deformation mechanical properties of carbon nanotubes and graphene.1−5 Studies of large deformation processes like fracture, on the contrary, have generally relied on the use of quantum mechanical (QM) methods because of their ability to model the breaking of chemical bonds.6−9 However, the generality and accuracy of QM calculations come at a much greater computational cost, which limits their use to relatively small model systems. This is a significant problem because, for a number of technologically important problems, it is important to be able to simulate large model systems over large deformations with near QM-quality results. The ReaxFF reactive force field is an empirical potential developed to directly simulate chemical changes, such as covalent © XXXX American Chemical Society

Received: June 19, 2015 Revised: August 28, 2015

A

DOI: 10.1021/acs.jpca.5b05889 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article


The Journal of Physical Chemistry A

AC predicted by the two parameter sets. All structural images were rendered using Ovito.27 Two 3D periodic MD models of diamond were built for use in this work. The first contained 8000 atoms and was oriented such that the primary crystallographic directions (⟨100⟩, ⟨010⟩, and ⟨001⟩) aligned with the x, y, and z axes of the model, respectively. For this structure, shown in Figure 1a, the simulation cell

particularly when it is applied to systems not considered during the fitting process. For example, the ReaxFFCHO parametrization was originally developed by Chenoweth et al.14 to model oxidation of hydrocarbons. Though it performs well in those simulations, recent work has identified a number of shortcomings of the ReaxFFCHO parametrization when it is used to simulate the mechanical response of carbon materials.19 In particular, it overestimates the Young’s moduli of graphene and carbon nanotubes and yields Poisson’s ratios for these materials and diamond that are in particularly poor agreement with both experimental and QM-derived values. The unreliable performance of the ReaxFFCHO parametrization for these systems is not particularly surprising because no mechanical response data for the condensed phases of carbon were used in its development. Srinivasan et al.20 recently published a new ReaxFF parametrization, referred to here as ReaxFFC‑2013, to address the mechanical property prediction problem and to make other improvements. The developers of ReaxFFC‑2013 sought to improve the mechanical predictions by constructing a new fitting set that included data on the bulk expansion and contraction of diamond and graphite, and uniaxial tension and compression data for diamond strained in the ⟨100⟩ direction. The objective of this study is to assess the accuracy of the ReaxFFC‑2013 parameter set in predicting the mechanical properties of a range of condensed carbon materials. To evaluate the performance of ReaxFFC‑2013 for materials with mechanical properties included in the fitting set, results are presented first for the elastic and fracture properties of diamond and graphene. For these systems, the ReaxFFC‑2013 results are benchmarked against both experimental data and DFT predictions, and compared with simulation results from the older ReaxFFCHO parametrization. After this evaluation, the transferability of ReaxFFC‑2013 to chemically related materials whose mechanical properties were not included in the fitting set is tested. Specifically, the elastic properties and fracture behavior of amorphous carbon (AC) and carbon nanotubes (CNTs) are simulated using ReaxFFC‑2013 and compared with experimental results.

Figure 1. (a) Unrotated and (b) rotated equilibrated MD diamond models.

dimensions were 3.57 nm in each direction. The second model was composed of 8064 atoms and was rotated so that the ⟨112⟩, ⟨110⟩, and ⟨111⟩ directions aligned with the x, y, and z axes of the simulation, respectively. For this system, shown in Figure 1b, the cell dimensions were 3.50 × 3.54 × 3.71 nm in the x, y, and z directions, respectively. The graphene model, shown in Figure 2, is composed of 6696 atoms. The initial simulation cell dimensions were

2. COMPUTATIONAL DETAILS The DFT simulations were performed using the VASP software package.21 The PAW−PBE plane-wave functional with Grimme’s DFT-D2 van der Waals correction term, referred to below as PBE-D2, was used.22−24 Two-atom primitive cells were used for both graphene and diamond calculations. A 15 × 15 k-point grid with an 800 eV energy cutoff was used for graphene whereas a 16 × 16 × 16 grid and 700 eV cutoff was used for diamond. In both cases, the k-point grids and energy cutoffs were carefully checked for convergence. All computed stresses were based on the virial stress measured at the deformed geometries and are, therefore, true stresses. All molecular dynamics simulations described in this work were done using ReaxFF, as implemented in the open source program LAMMPS.25,26 The ReaxFFC‑2013 parametrization was used for all systems considered in this study, and diamond and graphene were also simulated with the ReaxFFCHO parametrization to evaluate the relative performance of the ReaxFFC‑2013 parametrization. The ReaxFFC‑2013 parameter set used in this work was kindly provided by the authors of ref 20 in advance of publication. Slight adjustments were made to some of the parameters prior to final publication, but testing done with the final version revealed that the modifications had an insignificant impact on the properties of interest in this work. Figures S1−S4 in the Supporting Information demonstrate the close agreement in the stress−strain behaviors of diamond, graphene, CNTs, and

Figure 2. Equilibrated MD graphene model.

13.25 × 13.18 × 25.00 nm in the armchair, zigzag, and normal directions, respectively. The graphene model used 3D periodic boundary conditions and the out-of-plane box dimension was set to be sufficiently large to avoid interactions across the periodic boundary. Stresses were computed using an assumed graphene thickness of 0.34 nm. All MD simulations of graphene and diamond used the same equilibration procedure. Each system was first heated from 0 to 300 K over 5 ps using a Langevin thermostat. A Berendsen thermostat and barostat were then used to maintain the system at 300 K and zero pressure for 16 ps. The barostat was then removed and the system maintained at 300 K for 2 ps. A time step of 0.2 fs was used in all simulations. Because carbon nanotubes are found in both single and multiwall configurations, it was of interest to compare the mechanical properties of different configurations. The first model, composed of 14 single wall nanotubes (SWNTs) in a 2 × 7 hexagonal array and containing 26 880 atoms, is shown in Figure 3a. The SWNTs were of (20,0) chirality and were 1.5 nm in diameter. The second model, containing 4 multiwall nanotubes (MWNTs) in a 2 × 2 hexagonal arrangement and a total of 26 496 atoms, is shown in Figure 3b. The nanotubes in the MWNT array were built with three concentric CNTs with chiralities of (14,0), (23,0), B


Article



Figure 3. Equilibrated (a) SWNT array and (b) MWNT array.

and (32,0) and had an outer diameter of 2.5 nm. In both cases, the periodic cells were 10.3 nm long in the axial direction. Equilibration of the CNT systems was achieved by first heating the system from 0 to 300 K over 100 ps, followed by a temperature hold of 300 K using a Langevin thermostat and zero pressure using a Berendsen barostat. The thermostat was then switched to a Berendsen, the barostat removed, and the system maintained at 300 K for an additional 50 ps. The densities of the equilibrated systems were 1.26 and 1.99 g/cm3 for the SWNT and MWNT systems, respectively. The experimentally determined mechanical properties of AC vary widely, but they are known to be highly dependent on density. To determine if the simulation approach used in this work is capable of reproducing this behavior, three AC simulation models were built with densities spanning the range found in the experimental work. Specifically, the model systems, each containing 6720 atoms, were built at densities of 2.35, 2.82, and 3.37 g/cm3. This was done by first placing the carbon atoms at random positions within a simulation cell of the proper volume to yield initial densities of 2.50, 3.0, and 3.5 g/cm3, which were reduced to the values given above by the equilibration procedure. Each system was heated to 3000 K under constant volume conditions at a rate of 200 K/ps using a Berendsen thermostat, followed by cooling to 300 K at a rate of 5 K/ps. At 300 K, the volume constraint was removed and the system was equilibrated for 50 ps using a Berendsen barostat. Finally, the volume was fixed at the equilibrated value and the system was allowed to relax for an additional 50 ps at 300 K. All simulations conducted at temperatures above 300 K used a time step of 0.1 fs, which was increased to 0.2 fs for simulations at 300 K. Panels a−c of Figure 4 show representative AC structures equilibrated at densities of 2.37, 2.82, and 3.35 g/cm3, respectively, with the atoms colored according to the valence of the bonds around the carbon atoms. All atom pairs connected with a bond order greater than 0.3 were considered to be covalently bound. It is clear from Figure 4 that the ratio of sp3 to sp2 bonding increased significantly with increasing density. This observation was quantified by calculating the radial distribution function and reduced radial distribution function for each of the three AC models; these are shown in Figure 5a,b, respectively. Also shown in Figure 5 are corresponding experimental distribution functions derived from neutron diffraction data by Gilkes et al.28 Because of the extreme difference in system size between the simulations and experiment, the simulated peaks have been broadened to match the width of the first peak in the experimental radial distribution function by convolution with Gaussian functions with σ = 0.1125, 0.089, and 0.081 Å for the 2.37, 2.82, and 3.35 g/cm3 systems, respectively. Inspection of these plots reveals that the first peak shifts to the right as the density increases, indicating growth of the average carbon− carbon bond distance. This shift provides further evidence that the sp3 content of the system increases with increasing density.

Figure 4. Equilibrated AC systems with final densities of (a) 2.37 g/cm3, (b) 2.82 g/cm3, and (c) 3.35 g/cm3. In all images, sp3 bonded atoms are shown in red, sp2 in green, and sp in blue.

It can also be seen that the locations of the second, third, and fourth peaks are in good agreement with those in the experimental reduced radial distribution function. Finally, note that the experimental peak, measured for a system with a density of 3.00 g/cm3, falls between the simulations peaks for the 2.82 and 3.35 g/cm3 model systems. Readers seeking further discussion of the radial distribution functions for AC are referred to the recent work of Caro et al.29 C


Article



The stiffness coefficients for each material were determined using MD through a set of six simulations. For each simulation, the box and atomic coordinates were scaled in one of the six strain directions, three axial and three shear, whereas the remaining five strain components were set to zero. Each simulation therefore produced one column of the stiffness tensor. Strained systems were run for 10 ps at constant volume and a temperature of 300 K using a Berendsen thermostat. The potential energy and stress in the systems were found to plateau well within the 10 ps relaxation time. The pressure differences between the strained and unstrained systems were calculated using the final 5.0 ps of these simulations for all systems except for the CNT systems. Because of the low stiffness in the transverse and shear directions relative to the thermal fluctuations, it was necessary to run CNT systems for a total of 50 ps before and after straining. The pressure difference between the strained and unstrained systems was calculated using the final 40 ps of these simulations to obtain converged values. For all materials, strains of 0.10%, 0.25%, 0.50%, and 1.00% were investigated, and the AC systems were also strained to 2.00%. Using strains of 0.10% and 0.25% resulted in significant scatter and, in some cases, unphysical results, due to the magnitude of the thermal fluctuations relative to the small changes in stress. Strains of 0.50% and larger were found to produce reliable results and those results were averaged together here, as were the symmetric components of the stiffness tensor. The DFT results were also computed by straining the lattice in the three normal and three shear directions and computing the resulting stresses. The PBE-D2 stress−strain response simulations were computed through a series of static minimizations at fixed strains. Specifically, at each axial strain increment of 1%, the system was strained at five subincrements in the transverse direction and a parabolic fit of the five transverse strain energies was made. The energy corresponding to the minimum of this curve, which represents the zero transverse stress configuration, was used for subsequent mechanical property calculations. The ultimate tensile strength and fracture behavior of the models were determined in MD simulations by deforming each material in uniaxial tension until the measured stress dropped sharply, indicating fracture. During the tensile simulations, contractions transverse to the tensile axis were allowed by maintaining the pressures at zero in the transverse directions using a Berendsen barostat. The models were deformed at a true strain rate of 1.5 ns−1 while the temperature was maintained at 300 K using a Berendsen thermostat. This strain rate was selected on the basis of previous simulations of graphene, CNT, and diamond using the ReaxFFCHO parametrization.19 Additional simulations were performed to ensure that this strain rate was also appropriate for the AC system, as described in the Supporting Information. The Supporting Information also discusses the analysis used to select a time step 0.2 fs for the fracture simulations.

Figure 5. (a) Radial distribution function g(r) and (b) reduced radial distribution function G(r) of simulated AC structures compared to experimental results of Gilkes et al.28

The elastic mechanical response of each system was determined by computing the full stiffness tensor using an equivalentcontinuum mechanics approach30 where stresses and strains are related via the stiffness tensor as shown in eq 1. ⎡C C12 C13 0 0 0 ⎤⎡ ε11 ⎤ ⎡ σ11 ⎤ ⎢ 11 ⎥⎢ ⎥ ⎢σ ⎥ ⎢ C22 C23 0 0 0 ⎥⎢ ε22 ⎥ ⎢ 22 ⎥ ⎢ ⎥⎢ ε ⎥ ⎢ σ33 ⎥ ⎢ C33 0 0 0 ⎥ 33 ⎢ ⎥ ⎢ ⎥=⎢ σ C44 0 0 ⎥⎥⎢ 2ε23 ⎥ ⎢ 23 ⎥ ⎢ ⎢ ⎥ ⎢ σ13 ⎥ ⎢ symm C55 0 ⎥⎢ 2ε13 ⎥ ⎢ ⎥ ⎢ ⎥⎢ ⎥ σ ⎣ 12 ⎦ ⎢ C66 ⎥⎦⎣ 2ε12 ⎦ ⎣

(1)

The σij are components of the stress tensor, Cij are components of the stiffness tensor, and εij are components of the strain tensor. The stiffness tensor and engineering constants (Young’s modulus, Poisson’s ratio, and shear modulus) are inversely related as shown in eq 2, where Ei is the Young’s modulus in the i direction, vij is the Poisson’s ratio associated with the i and j directions, and Gij is the shear modulus in the ij plane. ⎡C11 C12 C13 0 0 0 ⎤−1 ⎢ ⎥ C22 C 23 0 0 0 ⎥ ⎢ ⎢ ⎥ C33 0 0 0 ⎥ ⎢ ⎢ ⎥ C44 0 0 ⎥ ⎢ ⎢ symm C55 0 ⎥ ⎢ ⎥ ⎢⎣ C66 ⎥⎦ ⎡ 1 − ν12 − ν13 ⎤ 0 0 0 ⎥ ⎢ E1 E1 ⎢ E1 ⎥ ⎢ ⎥ − ν23 1 ⎢ 0 0 0 ⎥ E2 E2 ⎢ ⎥ ⎢ ⎥ 1 ⎢ 0 0 0 ⎥ ⎢ ⎥ E3 ⎥ =⎢ 1 ⎢ ⎥ 0 0 ⎢ ⎥ G23 ⎢ ⎥ ⎢ ⎥ 1 symm 0 ⎢ ⎥ G13 ⎢ ⎥ ⎢ 1 ⎥ ⎢ ⎥ G12 ⎦ ⎣

3. RESULTS AND DISCUSSION The presentation of results is divided into two sections: the first describes the elastic and fracture response of diamond and graphene, whereas the second describes the elastic and fracture response of amorphous carbon and carbon nanotubes. The first section focuses on graphene and diamond systems, which represent complete sp2 and sp3 hybridized bonding, respectively, because they are relatively simple and well-represented in the ReaxFFC‑2013 training set. The ReaxFFC‑2013 parametrization results are compared to those from previous ReaxFFCHO parametrization simulations, PBE-D2 predictions, and experimental results. The second section focuses on the transferability of the ReaxFFC‑2013

(2) D


Article


components C11 and C12, via eq 3, which may be reduced to v = C12/C11. The Young’s modulus is reduced from 1024 GPa, with ReaxFFCHO, to 715 GPa, with ReaxFFC‑2013. For this property, the fortuitous cancelation of error in the ReaxFFCHO C11 and C12 stiffness tensor components results in a closer match to the experimental and PBE-D2 results. These results point to a need for further optimization of the ReaxFFC‑2013 parametrization with a focus on improving the C12 component. Finally, the lattice constant predicted by ReaxFFC‑2013 is also improved relative to ReaxFFCHO and is in good agreement with the PBE-D2 and experimental results. Experimental and calculated elastic properties of graphene are presented in Table 2. Because of the difficulties of measuring


parametrization to CNTs and to AC structures, which are characterized by bonding patterns that are mixtures of sp2 and sp3. 3.1. Diamond and Graphene. The experimental values for the stiffness tensor components of diamond assume a conventional crystallographic orientation with the ⟨100⟩, ⟨010⟩, and ⟨001⟩ crystallographic axes aligned with the x, y, and z directions. Therefore, only the stiffness tensor of the unrotated diamond simulation cell, shown in Figure 1a, was computed using ReaxFF. Diamond has cubic crystal symmetry, which reduces the number of independent stiffness tensor components to three: C11, C12, and C44. In this case, eq 2 reduces to eq 3. ⎡C11 C12 C12 0 0 0 ⎤−1 ⎢ ⎥ C11 C12 0 0 0 ⎥ ⎢ ⎢ ⎥ C11 0 0 0 ⎥ ⎢ ⎢ C44 0 0 ⎥ ⎢ ⎥ ⎢ symm C44 0 ⎥ ⎢ ⎥ ⎢⎣ C44 ⎥⎦

Table 2. Elastic Properties of Graphene C11 (GPa) C12 (GPa) C66 (GPa) E (GPa) υ lattice constant (Å)

⎡ 1 −v −v ⎤ 0 0 0⎥ ⎢ E E ⎢E ⎥ ⎢ ⎥ −v 1 0 0 0⎥ ⎢ E E ⎢ ⎥ 1 ⎢ ⎥ 0 0 0⎥ ⎢ E ⎥ =⎢ 1 ⎢ ⎥ 0 0 ⎢ ⎥ G ⎢ ⎥ 1 ⎢ ⎥ symm 0 ⎢ ⎥ G ⎢ ⎥ 1⎥ ⎢ ⎢⎣ G ⎥⎦

a

(3)

Table 1. Elastic Properties of Diamond

a

expt

PBE-D2

ReaxFFC‑2013

ReaxFFCHO

107931,32 12531,32 57931,32 1065a 0.116a 3.56733

1076 127 574 1061 0.118 3.563

1027 566 252 715 0.551 3.578

6431 5897 457 1024 0.917 3.619

PBE-D2

ReaxFFC‑2013

ReaxFFCHO

1067 148 459 1046 0.139 2.467

1005 505 186 751 0.502 2.462

5320 4662 311 1235 0.876 2.462

Computed from E and υ using eq 4.

the mechanical properties of graphene, a Poisson’s ratio and lattice constant measured on graphite are listed in Table 2 as reference values. An experimental Young’s modulus of 1000 ± 100 GPa has been reported for graphene,34 which is in good agreement with previous ab initio results of 1050 GPa.9 The Poisson’s ratio of graphene is often estimated to be the same as the in-plane Poisson’s ratio of graphite, which is 0.165.35 This is in good agreement with ab initio results, which predict a Poisson’s ratio for graphene of 0.164−0.186.9,36,37 Finally, an experimental shear modulus of 280 ± 36 GPa has been reported for graphene.38 The shear modulus is equivalent to the C66 stiffness tensor component and is listed as such in Table 2. If graphene is taken to be a planar 2D material, the stress and strain components transverse to the plane are eliminated in eq 1, and eq 2 reduces to eq 4. Like diamond, graphene is isotropic, which makes C11 = C22. Experimental stiffness tensor components in Table 2 are computed using eq 4.

Experimental and calculated elastic properties of diamond are shown in Table 1.

C11 (GPa) C12 (GPa) C44 (GPa) E (GPa) υ lattice constant (Å)

expt 1028a 170a 28038 100034 0.16535 2.45939

⎡ C11 C12 ⎢ ⎢C12 C11 ⎢ 0 ⎣ 0

Computed from C11 and C12 using eq 2.

⎡ 1 −v ⎤−1 0⎥ ⎢ E 0 ⎤ ⎢ E ⎥ ⎥ ⎢ −v 1 ⎥ 0 ⎥ =⎢ 0⎥ E E ⎥ ⎥ C66 ⎦ ⎢ 1⎥ ⎢ 0 ⎢⎣ 0 G ⎥⎦

(4)

The PBE-D2 results in Table 2 agree well with the experimental data, particularly the C11 and C12 tensor components and the lattice constant. The PBE-D2 and experimental results differ more substantially for the shear tensor component C66, although it is not clear if this arises from experimental uncertainties or from problems in the DFT calculations. If, for example, planewrinkling modes contribute to the lower experimental result, the DFT calculations would overestimate the shear modulus because a two-atom unit cell is not large enough to capture that deformation mode. These PBE-D2 results are, however, in good agreement with previous periodic DFT results, which report a shear modulus of 408 GPa using the PBE functional,40 as well as an experimental graphite in-plane shear modulus of 435 GPa.35

The experimental reference values for the Young’s moduli (E) and Poisson’s ratios (υ) shown in Table 1 were computed from the stiffness tensor components reported in refs 31 and 32, using eq 3. It can be seen that excellent agreement is found between experimental results and PBE-D2 predictions for every property listed in Table 1. Using these values as the basis for comparison, the ReaxFFC‑2013 parametrization is found to yield much more accurate C11 and C12 terms than those from ReaxFFCHO simulations. The ReaxFFC‑2013 Poisson’s ratio of 0.551, though still too large, is substantially improved in comparison to the ReaxFFCHO result of 0.917. This improvement in Poisson contraction response relates to the improvement in the stiffness tensor E


Article


included in the DFT-derived fitting set, more accurately than the older ReaxFFCHO parametrization. The greatest improvements were found in the C11 and C12 terms and, in the case of diamond, the lattice constant. There was a substantial improvement for the Poisson’s ratio in both materials, but additional optimization of the parameter set may be required to further improve agreement with experiment. The remainder of this section is dedicated to the fracture behavior of diamond and graphite. Because the fracture of diamond is highly sensitive to defects,41 there are very few experimental data available for use in validating the simulations, which are defect free. Because of this, the MD results will have to be compared with DFT predictions. However, because the simulation cell sizes that can practically be used in DFT calculations are too small to capture a physically realistic crack initiation and propagation process, the point of maximum stress on the stress− strain curve for diamond will, by necessity, have to serve as a proxy for the ultimate strength. Although there is clearly a need for more sophisticated DFT-based investigations of the diamond fracture process, they are beyond the scope of the present work. Although similar problems exist for DFT predictions of graphene fracture, experimental data are available for validating the MD simulations. The stress−strain responses for deformations along the ⟨111⟩, ⟨112⟩, ⟨110⟩, and ⟨100⟩ crystallographic directions of diamond are shown in Figure 7a−d, respectively. Each panel compares the results from simulations using the two ReaxFF parametrizations with the PBE-D2 results. Except for the ⟨100⟩ results shown in Figure 7d, both ReaxFF parametrizations yield similar results, although the ReaxFFCHO set shows a better agreement for the elastic modulus, which was expected on the basis of the results discussed above. For the deformation in the ⟨100⟩ direction, the ReaxFFC‑2013 parameter set better reproduces the large strain behavior predicted by the PBE-D2 calculations. This improvement stems from the inclusion of stress−strain data for this plane in the ReaxFFC‑2013 fitting set. The data shown graphically in Figure 7 are summarized as ultimate strengths and strains in Table 3. In addition to the results from both ReaxFF parameter sets and the current PBE-D2 predictions, Table 3 also includes the results from previous DFT calculations reported in the literature. Although


Comparing the ReaxFF results to the PBE-D2 and experimental results, the ReaxFFC‑2013 parametrization is again found to give dramatic improvements over the ReaxFFCHO results for both the C11 and C12 terms. As with diamond, ReaxFFC‑2013 yields a substantially improved treatment of the Poisson contraction response. The C66 component from both ReaxFF parametrizations is much closer to the experimental result than the PBE-D2 predictions are. The better agreement found for the MD results appears to be a result of the use of a much larger simulation cell, which is capable of capturing plane-wrinkling modes, as shown in Figure 6. Finally, the predicted lattice constant, which is the same

Figure 6. Wrinkle deformation modes for graphene sheared 1% using ReaxFFC‑2013. Shear strained box dimensions exaggerated for visibility.

in both of the ReaxFF parametrizations, is in good agreement with the PBE-D2 and experimental results. Considering the results presented to this point, it is clear that the new ReaxFFC‑2013 parametrization predicts many of the mechanical properties of diamond and graphene, the two materials

Figure 7. Comparison of ReaxFF and PBE-D2 stress−strain response strained along (a) ⟨111⟩, (b) ⟨110⟩, (c) ⟨112⟩, and (d) ⟨100⟩ directions of diamond. F


Article


graphene strained in the armchair direction and results in the flattened portion of the stress−strain response that starts around 0.16 strain. The magnitude of the potential energy bifurcation is substantially reduced with ReaxFFC‑2013 compared to ReaxFFCHO. When strain is applied in the zigzag direction, however, with ReaxFFC‑2013, no energy bifurcation is observed. In either strain direction, the ReaxFFC‑2013 parameters better reproduce the PBE-D2 predictions of ultimate stress than do the previous ReaxFFCHO parameters. Table 4 summarizes the ultimate stress and strain properties of graphene in the armchair and zigzag directions. Because

Table 3. Fracture Properties of Defect Free Diamond DFT (lit.) ⟨111⟩ ultimate strength (GPa) ⟨112⟩ ultimate strength (GPa) ⟨110⟩ ultimate strength (GPa) ⟨100⟩ ultimate strength (GPa) ⟨111⟩ ultimate strain ⟨112⟩ ultimate strain ⟨110⟩ ultimate strain ⟨100⟩ ultimate strain

PBE-D2

ReaxFFC‑2013 ReaxFFCHO

90,42 95,43 826 936

86

145

149

96

130

138

130,42 1146

119

153

179

225,42 2306

210

184

182

0.12 0.14 0.22 0.31

0.18 0.18 0.23 0.27

0.17 0.17 0.21 0.19

Table 4. Fracture Properties of Graphene


expt armchair ultimate stress (GPa) zigzag ultimate stress (GPa) armchair ultimate strain zigzag ultimate strain

there is some scatter in the previously published DFT results, the results obtained in the present work are consistent with them. The closest agreement between PBE-D2 and ReaxFF predictions is found for diamond strained in the ⟨100⟩ direction. For that system, both ReaxFF parametrizations yield ultimate strength values within 13% of the PBE-D2 result, although the ReaxFFCHO parametrization predicts failure at a much lower strain. When diamond is strained in the remaining three crystallographic directions, both ReaxFF parametrizations substantially overpredict the ultimate stress, exceeding the PBE-D2 predictions by 30% to 70%. This is not totally unexpected because only ⟨100⟩ strained diamond was included in the training set for ReaxFFC‑2013. The predicted stress−strain response of graphene is depicted in Figure 8a,b, which shows the results from the PBE-D2

130

34

0.2534

PBE-D2

ReaxFFC‑2013 ReaxFFCHO

114

125

134

139

138

96

0.24

0.24

0.21

0.27

0.28

0.20

the direction of applied strain could not be determined in the experimental work of Lee et al.,34 only one experimental value is given for stress and strain. For both ultimate stress and strain, the ReaxFFC‑2013 values agree very well with the PBE-D2 results and bracket the experimental number. The ReaxFFCHO set yields poorer predictions for both strain directions relative to both experiment and PBE-D2 results. To summarize, the results indicate that the ReaxFFC‑2013 predictions are in better agreement with PBE-D2 and experiment than are the predictions that used the previous ReaxFFCHO parameters. Both ReaxFF parameter sets overpredict the ultimate strength of diamond in several planes. Including these planes in the training set of future ReaxFF parametrizations may improve the fit. The most notable remaining deficiency of the ReaxFFC‑2013 parametrization is the systematic underprediction of the Young’s modulus. Incorporation of additional elastic data in the fitting set would likely alleviate this problem. 3.2. Carbon Nanotubes and Amorphous Carbon. The results in the previous section indicate that the ReaxFFC‑2013 parametrization is an improvement over ReaxFFCHO for most mechanical properties of diamond and graphene. It also makes property predictions that are in reasonably good agreement with experimental data and PBE-D2 results. Having developed a familiarity with the performance of the ReaxFFC‑2013 parameter set for diamond and graphene, in this section its ability to predict mechanical properties for CNTs and AC is evaluated. The bonding encountered in these materials is more complex than in either graphene or diamond, which are essentially pure sp2 and sp3, respectively. This provides a challenging test for the ReaxFFC‑2013 parametrization because neither of these materials’ mechanical responses were represented in the fitting data used to parametrize the force field. As in the previous section, the elastic properties of the two materials are presented first, followed by their fracture behavior. The predicted stiffness tensor components and corresponding engineering properties for the SWNT and MWNT systems are shown in Table 5. All nanotubes were aligned such that their axial direction corresponds to the subscript “1” in eq 5 and Table 5; the two equivalent transverse directions correspond to the subscripts “2” and “3”. Because the systems are transversely isotropic, eq 2 reduces to eq 5. The shear stiffness tensor component C66,

Figure 8. Predicted stress−strain response of graphene strained in the (a) armchair and (b) zigzag direction using MD and PBE-D2 techniques.

calculations and both ReaxFF parametrizations for strains along the armchair and zigzag directions, respectively. Consistent with experimental observations,34 all of the ReaxFF simulations result in brittle fracture. The stress−strain response predicted with ReaxFFCHO plateaus at around 0.1 strain in both strain directions. A previous work determined that this is associated with bifurcation of the potential energy into regions of high and low energy.19 As the system is strained, the high potential energy regions grow in size, resulting in a relatively flat stress−strain response. For the ReaxFFC‑2013 parametrization, the same behavior is observed for G


Article


this variability, three different model systems, spanning much of the experimental range of densities, have been constructed and simulated. The predicted stiffness tensor components and corresponding engineering properties of these models are shown in Table 6. The AC systems are isotropic, so the relationship

Table 5. Elastic Properties of CNT Systems


C11 (GPa) C22 (GPa) C12 (GPa) C23 (GPa) C44 (MPa) E1 (GPa) E2 (GPa) υ12 υ23

SWNT array

MWNT array

463 33 30 32 582 436 2 0.002 0.184

699 25 18 22 1193 685 5 0.003 0.044

Table 6. Elastic Properties of Amorphous Carbon Systems density (g/cm3) C11 (GPa) C12 (GPa) E (GPa) G (GPa) υ sp content sp2 content sp3 content

which corresponds to shear along the axis of the CNTs, was too low relative to the thermal fluctuations to calculate accurately (

Simulation of the Elastic and Ultimate Tensile Properties of Diamond

Recommend Documents