rsta.royalsocietypublishing.org

Research Cite this article: Glowacki DR, Rodgers WJ, Shannon R, Robertson SH, Harvey JN. 2017 Reaction and relaxation at surface hotspots: using molecular dynamics and the energy-grained master equation to describe diamond etching. Phil. Trans. R. Soc. A 375: 20160206. http://dx.doi.org/10.1098/rsta.2016.0206

Reaction and relaxation at surface hotspots: using molecular dynamics and the energy-grained master equation to describe diamond etching David R. Glowacki1,2,3 , W. J. Rodgers1 , Robin Shannon1,3 , Struan H. Robertson4 and Jeremy N. Harvey5 1 School of Chemistry, University of Bristol, Bristol BS8 1TS, UK 2 Department of Computer Science, University of Bristol, Bristol

Accepted: 11 January 2017 One contribution of 11 to a theme issue ‘Theoretical and computational studies of non-equilibrium and non-statistical dynamics in the gas phase, in the condensed phase and at interfaces’. Subject Areas: physical chemistry Keywords: heterogeneous chemistry, reaction dynamics, molecular dynamics, statistical mechanics, chemical vapour deposition, diamond Author for correspondence: David R. Glowacki e-mail: [email protected]

Electronic supplementary material is available at https://doi.org/10.6084/m9.figshare.c. 3685630.

BS8 1UB, UK 3 Department of Mechanical Engineering, Stanford University, 452 Escondido Mall, Stanford, CA 94305, USA 4 Dassault Systémes BIOVIA, 334 Cambridge Science Park, Cambridge CB4 0WN, UK 5 Department of Chemistry, KU Leuven, Celestijnenlaan 200F, 3001 Heverlee, Belgium DRG, 0000-0002-9608-3845 The extent to which vibrational energy transfer dynamics can impact reaction outcomes beyond the gas phase remains an active research question. Molecular dynamics (MD) simulations are the method of choice for investigating such questions; however, they can be extremely expensive, and therefore it is worth developing cheaper models that are capable of furnishing reasonable results. This paper has two primary aims. First, we investigate the competition between energy relaxation and reaction at ‘hotspots’ that form on the surface of diamond during the chemical vapour deposition process. To explore this, we developed an efficient reactive potential energy surface by fitting an empirical valence bond model to higher-level ab initio electronic 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

It has been known for some time that the dynamics of vibrational energy flow can play a pivotal role in determining the rates and pathways of unimolecular and bimolecular reactions in the gas phase [1–4]. The extent to which vibrational energy transfer dynamics can impact reaction outcomes in condensed phases remains an active research question. Over the past few years, a combination of molecular dynamics and master equation studies have provided examples of reaction outcomes in liquids which depend on energy transfer efficiencies [5–8]. Chemistry and reaction dynamics at the gas–surface interface [9–11] provide an interesting application domain for investigating the role that energy transfer plays in governing reaction outcomes because it offers an intermediate regime between the gas and liquid phases. Developing predictive chemical models which are capable of providing detailed microscopic insight into gas–surface kinetics and dynamics is an important area within chemical dynamics research [11]. Not only do a wide range of important industrial, catalytic and environmental processes occur at the gas–surface interface [12,13], but reaction dynamics studies at the gas– surface interface open up a wide range of fundamental research questions [11,14–16], including the following. To what extent is it possible for gas–surface kinetics and dynamics to be described using statistical models, or are they inherently dynamical processes? What role does energy transfer at the gas–surface interface play in determining chemical reaction outcomes [17]? And how does charge transfer and breakdown of the Born–Oppenheimer approximation impact surface reactivity? This paper focuses on the first two questions listed above—i.e. how well can statistical models describe reactivity at the gas–surface interface, and to what extent can we develop quantitative models of energy transfer? The energy transfer question is particularly interesting given that recent advances in both experimental and computational methods have revealed a range of gas–surface reaction systems which appear to exhibit some degree of mode selectivity [16]. For example, in dissociative chemisorption experiments of CH2 D2 on a metal surface, Beck et al. showed that dissociation was more likely if two quanta of energy were placed into a single CH stretching motion than if slightly more energy was distributed across two CH stretching motions—i.e. the more highly excited bond breaks more readily [18]. From a fundamental perspective, these studies provide insight into the time scales for energy redistribution and dissipation at the gas–surface interface, and also help us to understand the appropriateness of the ergodicity assumption in gas–surface dynamics. This in turn enables us to determine whether the

........................................................

1. Introduction

2

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

structure theory. We then ran 160 000 NVE trajectories on a large slab of diamond, and the results are in reasonable agreement with experiment: they suggest that energy dissipation from surface hotspots is complete within a few hundred femtoseconds, but that a small fraction of CH3 does in fact undergo dissociation prior to the onset of thermal equilibrium. Second, we developed and tested a general procedure to formulate and solve the energy-grained master equation (EGME) for surface chemistry problems. The procedure we outline splits the diamond slab into system and bath components, and then evaluates microcanonical transitionstate theory rate coefficients in the configuration space of the system atoms. Energy transfer from the system to the bath is estimated using linear response theory from a single long MD trajectory, and used to parametrize an energy transfer function which can be input into the EGME. Despite the number of approximations involved, the surface EGME results are in reasonable agreement with the NVE MD simulations, but considerably cheaper. The results are encouraging, because they offer a computationally tractable strategy for investigating non-equilibrium reaction dynamics at surfaces for a broader range of systems. This article is part of the themed issue ‘Theoretical and computational studies of nonequilibrium and non-statistical dynamics in the gas phase, in the condensed phase and at interfaces’.

(a)

(b)

3

SCH3 ∗ → S + CH3 ∗

SCH3 → SCH3

(R2a) (R2b)

The source of the internal energy indicated in (R2a) is H atoms, which combine with dangling CH2 radicals at the diamond surface, resulting in a local ‘hotspot’ containing approximately 100 kcal mol–1 excess energy. Within the nascent ‘hotspot’, it may be possible for CH3 dissociation (R2a) to occur before energy diffusion to the diamond bulk is complete (R2b), as illustrated in scheme 1. In what follows, we have examined the mechanistic proposal shown in scheme 1 using both MD simulations and non-equilibrium statistical mechanics calculations. For the former, we

........................................................

kinetics at surfaces is amenable to modelling using statistical rate theories (e.g. the energy-grained master equation (EGME)), or whether more detailed approaches based on explicit molecular dynamics are indeed required. In this paper, we investigate these questions by carrying out a detailed study of dissociation at the surface of a diamond slab shown in figure 1. Through comparisons between molecular dynamics (MD) simulations and EGME calculations, we show that dissociation at the diamond surface is in fact amenable to statistical treatments. There are two aspects which are key to the accuracy of the statistical approach outlined herein: (i) the system must be separated into ‘system’ and ‘bath’ components; and (ii) it must be possible to construct reasonably accurate models for energy transfer between the system and the bath. Our strategy for constructing an energy transfer model for surface reactions is straightforward: we have developed energy transfer functions which can be used within the EGME and which reasonably reproduce the time-dependent energy transfer profiles obtained in non-equilibrium MD simulations. Chemical vapour deposition (CVD) is a common way of making diamond in the laboratory; however, there remains significant uncertainty regarding the microscopic chemical mechanisms that drive deposition and etching kinetics at diamond surfaces. The specific motivation of the method described in this paper is to provide further insight into the etching behaviour observed in CVD experiments. In particular, it has been shown that etching occurs in atomic H atmospheres. For example, net growth of single-crystalline diamond is observed under CVD conditions using an H2 microwave plasma in the presence of CH4 , whereas in the absence of CH4 , etching is observed [19]. Spontaneous thermal dissociation of C–C bonds has been suggested as a possible etching mechanism; however, computational studies suggest that the barriers for C–C dissociation are too high to explain the observed etching rates [20]. An alternative explanation invokes non-equilibrium kinetics at the diamond surface according to the following mechanism (where ‘S’ indicates the diamond surface):

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

Figure 1. The C1498 H592 diamond [100] slab that was modelled in this work. (a) A ‘top-down’ view, and (b) a ‘side-on’ view. The surface is hydrogen-capped except for a single –CH3 group, which is located near the centre of the slab, and easier to see in the ‘side-on’ view. (Online version in colour.)

H

H

H

H H

* H

H

H

Scheme 1. Mechanistic proposal for etching at the diamond surface under CVD conditions. Hydrogen atoms combine with dangling CH2 radicals at the diamond surface, creating a local ‘hotspot’ with approximately 100 kcal mol−1 excess energy, which is energetically sufficient to dissociate a CH3 radical. carried out NVE MD simulations of a diamond slab. (In the microcanonical (NVE) ensemble, amount of substance (N), volume (V) and energy (E) are conserved.) C–C bond dissociation within these simulations was possible using an empirical valence bond (EVB) potential energy surface fitted to CASPT2 and CCSD(T) electronic structure theory. For the latter, we formulated and solved a stochastic EGME model that treats the competition between (i) CH3 dissociation (R2a) and (ii) energy diffusion out of the nascent hotspot into the bulk diamond (R2b). Energy transfer within the EGME model was treated using a function whose parameters were fitted to results obtained from the non-equilibrium MD simulations. As shown in what follows, the EGME (when compared to the results of MD simulations) is able to capture reasonably well the CH3 dissociation kinetics so long as the energy transfer time scales are appropriately modelled. Because diamond is an efficient heat conductor, energy transfer from the system to the bulk is very fast. Nevertheless, both the MD simulations and the EGME suggest a small probability (approx. 0.1%) that CH3 dissociation occurs prior to the completion of energy diffusion. This dissociation probability yields a phenomenological etching rate in line with experiment. The extent to which EGME models can be applied beyond the gas phase remains an open question, with a critical issue being how to treat energy transfer beyond isolated binary collision conditions. The approach that we have used in this work—i.e. fitting energy transfer parameters based on non-equilibrium MD simulations—is one which has previously been used with success to describe energy transfer both in liquids [21] and also in the gas phase [22]. To the best of our knowledge, this study represents one of the first attempts to apply this methodology to reaction dynamics at the gas–surface interface. The remainder of this paper is organized as follows. First, we describe the ab initio calculations we carried out in order to develop an accurate analytical –CH3 dissociative potential energy curve at the diamond surface. Second, we describe three different sets of MD simulations which we carried out in order to examine CH3 dissociation from the diamond surface. These included: (i) 160 000 non-equilibrium NVE trajectories in which the CH stretch of the surface methyl group was ‘plucked’ with the quantity of energy that would be available immediately following the H association step shown in scheme 1; (ii) thermal sampling along the –CH3 dissociation coordinate using the boxed molecular dynamics (BXD) method [23,24] in order to ascertain the free energy to dissociation and the corresponding thermal dissociation rate coefficient; and (iii) long NVE trajectories from which we backed out energy–energy correlation functions which we could use to fit the energy transfer parameters within our EGME model. Finally, we describe the EGME we formulated to model CH3 dissociation from the diamond surface, and show comparisons with the MD results.

2. Ab initio calculations for methyl dissociation A key component of the present study is the calculation of an accurate potential energy curve for the bond dissociation step (R2a). To avoid the expense of carrying out electronic structure

........................................................

H

H

4

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

H

+ H

A

B

C

D

5

relaxed geometriesa

rigid geometriesb

B3LYP/6-311g*

BDE0 (C–C)c 84.5

BDE(C–H)d 98.0

b

BDE(C–C)d 76.1

BDE0 (C–C)c 102.4

B3LYP-D3/6-311G*

88.4

98.7

80.1

n.a.

MP2-F12/aug-cc-pVTZf

99.0

100.5

90.6

116.7

90.7

99.6

82.4

108.0

93.7

100.5

85.4

110.7

.......................................................................................................................................................................................................... .......................................................................................................................................................................................................... ..........................................................................................................................................................................................................

CCSD-F12/aug-cc-pVTZf

..........................................................................................................................................................................................................

f

CCSD(T)-F12/aug-cc-pVTZ

..........................................................................................................................................................................................................

c

e

CASSCF/cc-pVTZ

75.7

n.a

67.3

CASPT2/cc-pVTZ

90.3c

n.a.

82.0e

91.0e

..........................................................................................................................................................................................................

108.1e

..........................................................................................................................................................................................................

a Energies obtained using the indicated method with the structures optimized using B3LYP. b Energies obtained using the indicated method, and with fragment structures frozen to correspond with the B3LYP geometries at the neopentane minimum. c Energies do not include zero-point energy corrections. d Energies include B3LYP/6-311g* zero-point energy corrections. e Dissociation limiting energies calculated for the supermolecule with r = 10 Å. CC f The cc-pVTZ basis set (not augmented) was used on hydrogen.

calculations on a large diamond slab like that shown in figure 1, we instead carried out electronic structure calculations on a sequence of increasingly larger diamond proxy model compounds, the sequence of which is shown in figure 2. Using a range of methods, we initiated our model calculations by examining the dissociation of neopentane (model compound A in figure 2) to CH3 + C(CH3 )3 . All density functional theory calculations reported in this work were obtained using the Gaussian 09 suite of programs [25], while all post-SCF (self-consistent field) and multi-reference calculations were obtained using the MOLPRO suite of programs [26]. According to the NIST webbook [27], the respective enthalpies of formation at 298 K of neopentane, methyl and t-butyl are –40.14 ± 0.15, 34.821 and 11.0 ± 0.7 kcal mol–1 , suggesting a bond dissociation enthalpy of 86.0 ± 0.8 kcal mol–1 . Calculations at the B3LYP/6-311G(d) level of theory predict that the dissociation energy at 0 K is 1.77 kcal mol–1 smaller than the dissociation enthalpy at 298 K, so we take the experimental value of the dissociation energy at 0 K to be 84.2 kcal mol–1 . Table 1 shows the results of calculated bond energies using a variety of methods. It can be seen that the best calculated value at the CCSD(T)-F12 level of theory is 85.4 kcal mol–1 , in very good agreement with the experimental value. The (H3 C)3 CCH2 –H bond energy was also calculated at the B3LYP and CCSD(T) levels of theory, and the corresponding energies are shown in table 1. In addition to the single point energy calculations given in table 1, we also carried out a set of scans along the C–C dissociation coordinate, plots of which are shown in figure 3. Because

........................................................

Table 1. Calculated neopentane C–C (and C–H) bond dissociation energies (BDEs) in units of kcal mol–1 .

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

Figure 2. Model compounds A–D on which we carried out electronic structure calculations, chosen to mimic a diamond slab with a dangling –CH3 group.

(a)

(b)

6

relative energy (kcal mol–1)

60

40 B3LYP CASSCF CASPT2

20

0 1.5

2.5 3.0 2.0 C–C distance (Å)

3.5

1.5

2.0 2.5 3.0 C–C distance (Å)

3.5

Figure 3. Scans along the –CH3 dissociation coordinate for neopentane, calculated at various levels of theory. Results obtained from (a) relaxed and (b) rigid scans.

neopentane dissociation goes to a singlet diradical at large separations, an accurate treatment of the wave function along the dissociation coordinate requires the use of multi-reference methods. Our approach was to carry out both CASSCF(2,2)/cc-pVTZ and CASPT2(2,2)/cc-pVTZ energy calculations in which a two-electron, two-orbital active space was chosen to correspond to the singly occupied p orbital localized on each carbon during dissociation. Figure 3a shows CASPT2 and CASSCF results obtained at structures derived from relaxed B3LYP/6-311G(d) scans along the C–C reaction coordinate constrained to a C3V symmetry. Table 1 shows the CASSCF and CASPT2 results at large separations (i.e. with rCC = 10 Å), using the B3LYP optimized geometries. Note that the B3LYP calculations involved careful testing for a broken-symmetry unrestricted solution at each bond length, with the closed-shell solution found to be more stable until rCC = 2.4 Å. The expectation value for the S2 operator for the Kohn–Sham orbitals then increases rapidly from 0.19 at rCC = 2.6 Å to 0.92 at rCC = 3.4 Å. It should also be noted that the CASSCF and CASPT2 curves show no barrier to dissociation in excess of the endothermicity. The zero-point-corrected CASPT2 bond dissociation energy (BDE) of 82.0 kcal mol–1 is in good agreement with both experiment and the CCSD(T)-F12 BDE. Figure 3b shows the results of CASSCF and CASPT2 rigid scans along the C–C dissociation coordinate of neopentane. During these scans, C3V symmetry was enforced, and the starting structure was the optimized B3LYP structure for neopentane. At large separations, the energies shown in figure 3b are larger than those shown in figure 3a, as the forming methyl and t-butyl radicals are not allowed to relax to their optimum structure. For example, the CASPT2 rigid-scan BDE (taken as the potential energy for rCC = 10 Å) is 108.1 kcal mol–1 , which is 17.8 kcal mol–1 larger than the energy obtained using B3LYP optimized geometries (i.e. a relaxed scan). Table 1 shows that the BDEs computed without fragment relaxation are all approximately 18 kcal mol–1 greater than the relaxed bond energies. The neopentane dissociation data in table 1 and figure 3 were used to fit an analytical potential energy surface model to accurately represent (R2a), which is discussed in further detail below. In order to investigate the effect that the total system size had on the calculated –CH3 dissociation energy, we performed additional calculations on the incrementally larger diamond proxy models (structures B, C and D) in figure 2. Table 2 reports C–CH3 and CH2 –H BDEs calculated for each of these structures using the B3LYP-D3 level of theory. The results in table 2

........................................................

80

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

100

BDE(C–H)

B3LYP-D3 BDE(C–C)

model B

87.1

BDE(C–H)

96.4

90.8

97.1

model C

77.1

90.5

82.4

91.1

model D

73.6

90.8

81.3

91.8

.......................................................................................................................................................................................................... .......................................................................................................................................................................................................... ..........................................................................................................................................................................................................

show that the smaller model B leads to an overestimate of the C–CH3 and –CH2 –H BDEs at the diamond [100] surface, due to the lack of steric hindrance at the surface. In this respect, models C and D should be closer to the energies that we would expect at the surface of the diamond slab. The data in table 2 suggest that this is in fact the case, as the calculated bond energies appear to be converging to a limit as the system size increases. Note that the B3LYP bond energies are smaller than the B3LYP-D3 energies (especially for BDE(C–C)), owing to neglect of the attractive dispersive interactions between the methyl radical and the surface. A similar effect can also be seen for the C–C bond energy of neopentane in table 1. It is also worth noting that B3LYP-D3 underestimates the neopentane BDEs compared to both the experimental and the CCSD(T)-F12 values (both of which are in good agreement). Assuming that B3LYP-D3 makes a similar underestimate in the bond energy at the diamond [100] surface, our best estimates of the surface C−CH3 BDEs and CH2 −H BDEs (based on the model D B3LYP-D3 value and the difference between the B3LYP-D3 and CCSD(T)-F12 neopentane values) are 86.6 and 93.6 kcal mol−1 , respectively.

3. Analytical potential energy surface for CH3 dissociation at the diamond surface Non-equilibrium etching dynamics leading to dissociation of –CH3 from the diamond slab shown in figure 1 and scheme 1 turns out to be an extremely rare event. Accumulating the statistics required in order to assess the probability of this event required us to run 160 000 trajectories. This required an efficient dissociative potential energy surface (PES) with accurate energetics. Ab initio direct dynamics proved too computationally demanding for our purposes. Hybrid methods such as combined quantum mechanical and molecular mechanical (QM/MM) approaches might be tractable for optimizations and energy calculations [20,28–30], but are likewise too computationally demanding to carry out a large number of trajectory simulations. Instead, we chose to build on our previous work to date [8,31] using a parallel EVB [32] algorithm which we have implemented within the CHARMM molecular dynamics software [33]. Our preferred way to use the EVB method relies on fitting to higher-level electronic structure theory calculations. In this case, we chose to fit to the CASPT2 calculations shown in table 1 and figure 3, with a few caveats: it turns out that the CASPT2 method returns a relative energy at large C–C distances that is slightly smaller than our best estimate of the bond energy, discussed above. Furthermore, the CASPT2 energies along the dissociation curve cannot readily be corrected for zero-point energy. Accordingly, we scaled the CASPT2 BDEs along the figure 3 dissociation curve by D0 (CCSD(T))/De (CASPT2), where D0 (CCSD(T)) is 85.4 kcal mol–1 (the zero-point-energycorrected bond energy computed using CCSD(T)-F12) and De (CASPT2) is 90.3 kcal mol–1 (the bond energy calculated using CASPT2 without zero-point correction). This yields a set of relative energies Eab initio (ri ), which are shown in figure 4a, and which were used to carry out the PES fitting described below.

........................................................

B3LYP BDE(C–C)

7

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

Table 2. Zero-point-corrected BDEs (kcal mol–1 ) obtained for the diamond proxy models B–D shown in figure 2, at the B3LYP/6-311G(d) and B3LYP-D3/6-311G(d) levels of theory. (Geometry optimization and zero-point energy corrections were calculated at the B3LYP/6-31G(d) level of theory.)

(a)

(b)

8

90 relative energy (kcal mol–1)

60 50 40 30 20 10 0

ab initio energy EVB fit 1.5

2.0 2.5 3.0 reaction coordinate (Å)

diamond scan 3.5

1.5

2.0 2.5 3.0 reaction coordinate (Å)

3.5

Figure 4. (a) The results of relaxed scans along the neopentane C–C bond, comparing the energies obtained from E ab initio and λ0 of the fitted EVB pseudo-Hamiltonian in equation (3.1). (b) The results of a relaxed scan along the C–C bond in the diamond slab shown in figure 1, using Hamiltonian parameters derived from the neopentane fit depicted in panel (a). (Online version in colour.)

Within the EVB approach, the potential energy for a given set of nuclear coordinates q is obtained as the lowest eigenvalue λ0 of a pseudo-Hamiltonian matrix,   V1 (q) + ε1 H12 (q) H(q) = . (3.1) V2 (q) + ε2 H12 (q) The diagonal elements Vj (q) are molecular mechanics energies for the product and reactant structures. In this work, the diagonal elements of the matrix, V 1 (q) and V 2 (q), were obtained from the Merck molecular mechanics force field (MMFF94) [34] using appropriate bonding terms for the bonded and unbonded states, and modified to treat sp2 -hybridized carbon radical centres like CH3 .1 The role of the energy offsets εj (which are constant for all structures) is to correctly describe the electronic energy difference between reactant and product states. In the case that the respective coupling elements, H12 (q = P) and H12 (q = R), are close to zero near the product and reactant geometries, then the reaction energy is (V 2 (q = P) + ε2 ) − (V 1 (q = R) + ε1 ). We represented offdiagonal elements H12 (q) using a linear combination of N Gaussian functions (each of which depended on the C–C distance r) for the dissociating bond. The Gaussians had the form     N  1 r − Bk 2 Ak exp − H12 (r) = , (3.2) 2 Ck k=1

where Ak , Bk and Ck are the respective amplitude, centre and width parameters for a particular Gaussian function. In the present case, the potential was fitted to the electronic structure data mentioned earlier using a combination of two Gaussian functions for the off-diagonal term (i.e. N = 2 in equation (3.2)). The fit was chosen to minimize a least-squares metric, χ 2 (Ak , Bk , Ck ; 1 ≤ k ≤ N) =



[λ0 (ri ) − Eab initio (ri )]2 .

(3.3)

i 1

We modified the MMFF protocol so that it did not throw errors if an sp2 -hybridized carbon centre was bonded to three hydrogens. The CH bonding parameters in CH3 were taken to be identical to those used for the CH bonds in =CH2 . Similarly, the improper torsion parameters for CH3 (which ensure a planar equilibrium geometry) were taken to be identical to those for =CH2 .

........................................................

70

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

80

Table 3. Optimized parameters for the Gaussian functions making up H12 in equation (3.2). B (centre)

C (width)

1

38.35

2.32

0.27

2

133.34

2.95

0.44

.......................................................................................................................................................................................................... ..........................................................................................................................................................................................................

In this expression, λ0 (ri ) is the optimized value of the lowest eigenvalue of the matrix of equation (3.1), with the C–C distance r frozen to the same value ri as that in the electronic structure calculations, but with all other structural parameters optimized. The value of the offset ε1 was chosen such that V 1 + ε1 is equal to zero at the minimum of the potential energy curve, and such that V 2 + ε2 is equal to the ab initio dissociation energy at large r. The best-fitting parameters for each Gaussian are shown in table 3. The value of H12 is close to zero at distances r that are less than 1.54 Å and greater than 4.5 Å—i.e. the potential is described almost entirely by V 1 in the reactant region and by V 2 in the product region. A comparison between the fitted λ0 and Eab initio as a function of r is shown in figure 4a. Owing to the difficulty of calculating ab initio values for the potential energy curve of CH3 dissociation from anything larger than structure D in figure 2, the EVB potential for CH3 dissociation was constructed as follows: The functional form and parameter values required to specify H12 were obtained from fits to ab initio energies along the neopentane dissociation curve (shown in figure 4a); the offset parameters ε1 and ε2 were specified so that V 1 + ε1 is equal to zero and V 2 + ε2 equals the –CH3 dissociation energy of structure D (see figure 2) at large r. A relaxed scan along the C–C reaction coordinate for the final parametrized diamond EVB model is shown in figure 4b. The fact that the neopentane fits in figure 4a can produce reasonable dissociation curves at the diamond surface (figure 4b) indicates that the fitted EVB parameters have some degree of transferability.

4. Molecular dynamics simulations of –CH3 dissociation at the diamond surface (a) NVE simulations Having obtained an accurate and efficient –CH3 dissociation potential at the diamond surface, we investigated the relative rates of processes (R2a) and (R2b) using non-equilibrium NVE MD simulations. Here, ‘non-equilibrium’ refers to the fact that the initial conditions of each NVE trajectory were modified from coordinates and velocities selected from an NVT ensemble (details below), in order to emulate the nascent vibrationally excited SCH3 * species that we expect would be formed upon hydrogen atom addition to SCH2 . (In the canonical (NVT) ensemble, amount of substance (N), volume (V) and temperature (T) are conserved.) MD trajectories were propagated in order to observe the fraction of trajectories in which prompt methyl formation was observed prior to the establishment of thermal equilibrium. The procedure we used to carry out these MD simulations was as follows: A 1.6 µs (1 fs time step) NVT simulation of SCH3 was carried out using a Langevin thermostat with a friction coefficient of 10 ps–1 and heat bath of 1300 K (a typical temperature at which carbon etching is observed experimentally under [8] CVD growth conditions). Temperature equilibration was achieved within a few picoseconds, after which time coordinates and velocities from this trajectory were sampled every 10 ps (our tests showed this time to be a sufficiently long interval to ensure that the initial conditions were not significantly correlated to one another). This allowed us to generate initial conditions for 160 000 NVE simulations. Prior to running the NVE simulations, a non-equilibrium kinetic energy impulse was applied to one of the hydrogen atoms on the methyl group by modifying its NVT velocity in the direction of the CH bond, a local mode excitation strategy that we have used in previous work [8,35]. The magnitude of this velocity perturbation was chosen to correspond to a particular quantity of kinetic energy (KE), where

........................................................

A (amplitude)

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

Gaussian

9

Table 4. The fraction of prompt CH3 dissociation from the surface as a function of three different energetic √impulses applied to the CH stretch. A Poisson distribution of rare events (N) was assumed, where the error is proportional to N. CH3 dissociation events (N) 11

% dissociation 0.018 ± 0.006

..........................................................................................................................................................................................................

128

60 000

65

0.065 ± 0.010

148

40 000

92

0.092 ± 0.016

.......................................................................................................................................................................................................... ..........................................................................................................................................................................................................

KE = EC–H + EH (T) + Eexcess . The C–H bond energy, EC–H , was specified to be 93.6 kcal mol–1 (our best guess for the diamond surface C–H BDE, discussed in relation to table 2). Energy EH (T) is the typical thermal translation energy for a hydrogen atom at 1300 K. And Eexcess was an additional quantity of energy added to some of the runs in order to (i) increase the likelihood of a relatively rare –CH3 dissociation event, (ii) explore the sensitivity of the prompt –CH3 dissociation to initial energy and (iii) provide further data points for comparison with the non-equilibrium statistical mechanics approach described below. In the end, we ran three sets of NVE trajectories, each with a kinetic energy impulse KE = 108, 128 and 148 kcal mol–1 . The non-equilibrium NVE trajectories were run for a total of 10 ps (1 fs time step). All of the dynamics were carried out using the leapfrog Verlet integration scheme, and the 858 atoms at the boundaries of the C1498 H592 diamond slab model (figure 1) were frozen, leaving 1232 atoms that were free to move (for 3696 total degrees of freedom (d.f.)). As discussed below, the maximum excess energy added during our MD simulations was 148 kcal mol−1 . In the limit that all of this energy is equipartitioned into the diamond slab on the time scale of a single MD simulation, this corresponds to an average of 0.04 kcal mol−1 per degree of freedom (i.e. 148.32 kcal mol−1 /3696 d.f.), approximately 50 times smaller than kT (1.99 kcal mol−1 ) at 1000 K. This gave us confidence that our simulated diamond slab was large enough so as not to give rise to MD artefacts due to the system heating up as a result of adding non-equilibrium excess energy. The trajectories were subsequently analysed in order to identify prompt –CH3 dissociation events. Any trajectory in which the dissociating C–C bond length exceeded 3 Å was counted as a dissociation event (inspection of dissociative trajectories showed that this was an effective point of no return). The results of the simulations are shown in table 4. Dissociation events were indeed observed over this energy range, which shows that (for a small number of trajectories) surface dissociation can occur prior to energy dissipation. In general, the number of dissociation events increases with the internal energy introduced via the non-equilibrium impulse applied to the C–H stretch. Despite the fact that these dissociation events are rare, there were enough events to count, which allowed us to derive approximate Poisson error estimates. The low probability of –CH3 dissociation arises from the fact that the excess energy in the C–H stretch is quickly dissipated into the bulk of the diamond slab (a point which is discussed in further detail below, with a typical time-dependent plot shown as electronic supplementary material, figure S1). When dissociation did take place, it necessarily occurred prior to dissipation of the initial C–H stretch energy to the bulk, and was therefore an extremely rapid process—i.e. within 100–200 fs.

(b) Free energy sampling and thermal dissociation rate coefficient The results in table 4 give the fraction of prompt dissociation events observed when the system has an energy which is initially in excess of the –CH3 BDE, and therefore provide insight into the relative rate of C–H energy dissipation versus CH3 dissociation. The energy dissipation rates in this system are extremely high. All of the observed dissociation events take place before the excess energy is dissipated and the system achieves thermal equilibrium, and are therefore very fast. In order to get a handle on how the time scales observed for the prompt dissociation events compare to the time scales observed for thermal dissociation, we calculated the absolute rate of the thermal CH3 dissociation (i.e. the dissociation rate coefficient at typical thermal energy distributions). We note that the time scales of the NVE simulations described above are too short to observe any

........................................................

trajectories 60 000

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

energy (kcal mol−1 ) 108

10

30

11

energy (RT)

15 700 K 800 K 900 K 1000 K 1100 K 1200 K 1300 K

10

5

0 1.5

2.0

2.5 3.0 reaction coordinate (Å)

3.5

4.0

Figure 5. Temperature-dependent free energy surfaces of –CH3 dissociation at the diamond surface obtained using BXD at a range of temperatures. (Online version in colour.) thermal dissociation events. In order to back out the kinetics, we therefore used an accelerated free energy sampling technique called boxed molecular dynamics (BXD) [23,24,36]. The idea in BXD is to define a collective variable (or set of collective variables) which describes reaction progress, and then splice it into a set of ‘boxes’, which essentially correspond to hard-wall potential boundaries. A box is defined as a region of configuration space that lies between two boundaries; within any given box, a trajectory runs on the unbiased PES. If the trajectory crosses a particular box boundary, a velocity inversion operation [36] is performed to keep it within the specified box. BXD simulations are run by locking the system within a set of adjacent ‘boxes’, and then performing statistical analysis of the time spent in each box, and the relative number of hits at the boundaries that define the box. These quantities define box-to-box rate coefficients, which can then be used to calculate a potential of mean force, which is independent of the boundary locations. Choosing BXD boundaries is analogous to the process of specifying umbrellas (in umbrella sampling). A key difference is the fact that umbrella sampling requires two parameters per umbrella (location and force constant), whereas BXD requires only one parameter (location). While the latest BXD algorithms are able to adaptively determine the most computationally optimal box boundary locations [36] without requiring user specification, this capability was not available at the time the work described in this paper was carried out, and BXD boundaries were specified the ‘oldfashioned’ way (via user trial and error). In order to calculate thermal dissociation rate coefficients, we performed BXD simulations from 700 K to 1300 K in 100 K intervals. The collective variable which we used to define progress along the –CH3 dissociation coordinate was the C–C bond distance. This coordinate was split into user-specified boxes of varying sizes, depending on the C–C bond distance: between 1.4 and 2 Å the boxes were 0.1 Å wide; between 2 and 3 Å the boxes were 0.2 Å wide; and between 3.3 and 4.0 Å the system was defined by a single box. Simulations were then run in each of these boxes, using an NVT ensemble with a Langevin heat bath (T = 1300 K, friction = 15 ps–1 ). Each simulation used a time step of 0.25 fs, and lasted for 187.5 ps. These settings ensured convergence of the BXD sampling within each box. The free energy profiles calculated from BXD at each temperature are shown in figure 5. All of the free energy surfaces show the expected minimum at the C–C equilibrium bond length, then a

........................................................

20

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

25

12

12

10

ln(k)

4 2 0 –2 –4 –6 0.7

0.9

1.1 103/T (K–1)

1.3

1.5

Figure 6. Arrhenius plot of table 5 rate constants. The points indicate the rate coefficients calculated from the Eyring equation using the BXD free energies of activation.

Table 5. Free energy data obtained from BXD free energy profiles. The energies obtained are in units of RT and kcal mol−1 ; the reported rate constants were obtained by plugging the reported activation free energies into the Eyring equation. temperature (K)

activation free energy (RT)

activation free energy (kcal mol−1 )

rate constant (s–1 )

700

27.84

38.72

1.2 × 101

800

24.42

38.82

4.1 × 102

900

21.28

38.06

1.1 × 104

1000

18.72

37.20

1.5 × 105

1100

16.86

36.85

1.1 × 106

1200

14.88

35.48

8.6 × 106

1300

13.55

35.00

3.5 × 107

.......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... ..........................................................................................................................................................................................................

steep increase, followed by a maximum at rC–C ≈ 3.4 Å, followed by a decrease corresponding to the gain in entropy as the methyl radical dissociates. The maximum corresponds to the variational transition state for methyl loss. The origins of the variational transition-state barrier are primarily entropic (since there is no barrier on the underlying PES, see figure 4), and easiest to understand by thinking about the reverse process (i.e. CH3 association on the diamond surface): the barrier reflects the entropy penalty that CH3 pays for forming covalent bonds at the diamond surface. Table 5 gives Ga (T), the free energy of activation for CH3 dissociation at the diamond surface, calculated as the difference between the free energy at the variational transition state and the free energy minimum near the C–C equilibrium bond distance. Table 5 also shows the calculated thermal dissociation rate coefficients calculated using the Eyring equation, where k(T) = (kB T/h) exp(−Ga (T)/RT). Figure 6 shows an Arrhenius plot of k(T) plotted against 1/T. In a recent kinetic Monte Carlo modelling study of diamond growth under CVD conditions [37], it was found that good agreement with observed rates of growth could be found when assuming a rate constant for etching provided with an Eyring expression together with a temperature-independent free energy of activation of 200 kJ mol−1 (48 kcal mol−1 ). While this is

........................................................

6

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

8

The results discussed above provide significant microscopic insight into the relative time scales of prompt (non-equilibrium) dissociation of CH3 at the diamond surface, and thermal (equilibrium) dissociation. However, the computational cost of the simulations outlined in the previous section is considerable: in addition to the thermal free energy sampling which we carried out using BXD, we also ran 160 000 non-equilibrium NVE trajectories. A key aim of the work described in the final part of this paper is to investigate whether it is possible to use a more computationally tractable approach for investigating non-equilibrium chemistry at surfaces, which avoids the need for large numbers of NVE trajectories, but which is able to provide results that are in line with the NVE trajectory predictions. To this end, we formulated a non-equilibrium EGME model [38] to describe the competition between cooling and dissociation of the vibrationally excited SCH3 * at the diamond surface. All of the master equation and statistical mechanics calculations described below were carried out using MESMER [38], which is an open-source, cross-platform master equation code that we have been actively developing over the past few years. The master equation modelling approach outlined below builds on previous work in which we have carried out studies of non-equilibrium reaction dynamics in liquids [6]. The EGME is a computationally efficient framework typically used to describe the competition between relaxation and reaction dynamics for small molecules in the gas phase. It is essentially an energy-resolved system–bath model which describes the competing reaction and relaxation dynamics on a set of interlinked chemical intermediates. For this particular system, where we aim to model the competition between (R2a) and (R2b), the ansatz is as follows:  E=∞  E=∞ dp(E) = F(E ← E )p(E ) dE − F(E ← E)p(E) dE − kd (E)p(E), (5.1) dt E=0 E=0 where p(E) is the population of SCH3 at a particular energy E, F(E ← E) is a function that describes the probability that the SCH3 population at energy E is transferred (via bath interactions) to the SCH3 population at energy E, and kd (E) is the energy-resolved rate coefficient for dissociation of SCH3 at a particular energy to S + CH3 , typically calculated as  E=E E=E0

kd (E) =

ρTS (E ) dE

hρSCH3 (E)

,

(5.2)

where ρ SCH3 (E) is the density of states of SCH3 and ρ TS (E) is the density of states of the CH3 dissociation transition state. Equation (5.2) works well when the reaction path in question has a well-defined PES barrier; however, it was not a good choice for the present case, where figure 4 shows that the PES does not have a well-defined barrier. As a result, we used an alternative approach, which allows a set of microcanonical rate coefficients kd (E) to be calculated from a modified Arrhenius rate expression for the corresponding canonical rate coefficients k(β), where β = (κT)−1 . So long as k(β) can be represented by the modified Arrhenius expression,  k(β) = A0

β0 β

n exp(−βEa ),

(5.3)

........................................................

5. Master equation model of non-equilibrium dissociation

13

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

larger than the values shown in table 5, it should be noted that this value in fact corresponds to the activation free energy for loss of an adsorbed ‘carbon’, which is typically present as a −CH2 group inserted into a surface C–C bond. In the presence of excess H atoms in the diamond growth plasma, this can convert to a surface −CH3 group, but this process is endothermic by approximately 80 kJ mol−1 (19 kcal mol−1 ) [20]. This suggests a free energy barrier for breaking the surface −CH3 bond of about 29 kcal mol−1 , which is in satisfactory agreement with the BXD results in table 5, considering the uncertainty in the experimental analysis of the etching rate, and the inaccuracies of our EVB potential.

Diagonalizing the matrix M provides a solution for p, the population vector p(E), i.e. p(t) = U exp(Λt)U−1 p(0),

(5.6)

where U is the matrix of eigenvectors of M, and Λ a diagonal matrix of the corresponding eigenvalues. The EGME has advanced to the point that it is routinely able to provide nearquantitative predictions of gas-phase kinetics for reactive systems comprising multiple wells and transition states. Part of the reason for its success in the gas phase arises from the fact that the system/bath boundary is relatively well defined: the system comprises the reactive molecule of interest (which has available to it a range of reactive pathways), and the bath comprises third bodies which undergo collisions with the reactive system of interest. For example, the atmospheric bath is made up of N2 and O2 , both of which can undergo elastic and inelastic collisions with the reacting species. Third-body collisions of this sort can be treated (often to a relatively high degree of accuracy) using standard functional forms that have their origins in isolated binary collision models, quantum mechanical scattering calculations and classical trajectory simulations. Two important questions that arise in trying to apply system–bath models like the EGME to systems other than the gas phase are as follows: (i) How do we define the system/bath boundary? (ii) How do we treat energy transfer probabilities between the system and the bath? The approach that we have taken in this work is to use a simple distance-based criterion to define the system, which is shown in figure 7. We simply defined the reactive ‘system’ to comprise those atoms which are no more than four covalent bonds away from the C–H stretch to which the excess energy was initially added in the NVE simulations. Atoms which are more than four atoms away from the reactive system are treated as part of the bath. This definition led to a reactive system with 16 atoms. The microcanonical rate coefficients kd (E) required for solving equation (5.1) were obtained via equation (5.4) and require specification of the SCH3 density of states ρ(E). These were calculated in a locally modified version of CHARMM as follows: (i) a geometry optimization was carried out on the full diamond slab (system plus bath) shown in figure 4; (ii) the bath atoms in figure 7 were then frozen; (iii) a Hessian calculation was performed on the figure 7 system atoms; (iv) the system Hessian was then diagonalized, yielding a set of 3n vibrational frequencies (once the translational and rotational degrees of freedom were removed); and (v) energy-dependent densities of states for SCH3 ρ(E) were then calculated using the Beyer–Swinehart algorithm [38] within the harmonic oscillator approximation for the vibrational degrees of freedom only (i.e. rotational densities of states were set to zero). The kd (E) calculated using this procedure are shown in figure 8. The energy transfer function F(E ← E ) required to solve equation (5.1) was obtained by fitting to results obtained from MD simulations, together with linear response theory. Specifically, a single 1 ns NVE simulation was performed (time step = 1 fs), using the full reactive EVB potential energy surface, and starting from an equilibrated structure generated during an NVT simulation at 1300 K. Atomic velocities were saved at each time step during the simulation, and the total

........................................................

The A0 , B0 , n and Ea parameters in equations (5.3) and (5.4) were determined by fitting to the data in figure 6, and are discussed in further detail in the electronic supplementary material; in equation (5.4) Q(β) corresponds to the SCH3 partition function. Equation (5.1) represents a set of coupled differential equations which is typically discretized into a set of contiguous intervals or ‘grains’ with energy spacing E and written as a matrix equation: d p = Mp. (5.5) dt

14

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

then it is possible to calculate kd (E) using an inverse Laplace transform (ILT) (where L−1 [ . . . ] indicates an inverse Laplace transform of some argument [ . . . ]), i.e. [38,39]   1 (5.4) kd (E)ρ(E) = A0 β0 n L−1 Q(β) n exp(−βEa ) . β

15 ........................................................

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

Figure 7. The 16 ‘system’ atoms at the surface of the diamond slab (including the pendent CH3 ) are shown as spheres, which are embedded in bath atoms (shown as a wireframe structure). (Online version in colour.)

1012

k(E) (s–1)

1010 108 106 104 102 1 0

40

80

120

160

200

240

energy (kcal mol–1)

Figure 8. Value of k(E) for SCH3 dissociation to S + CH3 at the diamond surface, calculated via ILT as described in the text. kinetic energy of the atoms within the reactive subsystem was then computed from these velocities, in line with the strategy we have used successfully in previous work [8,21,31,40]. Specifically, the atomic kinetic energies were summed to construct a kinetic energy time series KE(t) for the ‘system’ atoms in figure 7. The average total energy within the system as a function of time E(t) was calculated by the virial theorem as 2 KE(t)τ , where the angled brackets indicate averages, and τ is a user-specified time window over which the averaging is carried out. This is a strategy that we have used successfully in a number of previous studies [5,8,31,40], which has given quantitative agreement with experiment, so long as τ is chosen to span several of the slowest vibrational periods within the set of ‘system’ atoms. For this system, we found that our results more or less converged so long as τ was greater than 100 fs. The E(t) time series was

16

MD correlation function

700

energy (kJ mol–1)

500

400

300 10–17

10–16

10–15

10–14 10–13 time (s)

10–12

10–11

Figure 9. Comparison of energy correlation function obtained from MD (line) with that obtained through fitting to the EGME as described in the text (points). (Online version in colour.)

then used to determine C(t), the time correlation function describing the total energy within the ‘system’ atoms, C(t) = δE(0)δE(t) = E(0)E(t) − E2 .

(5.7)

We note that, for this non-equilibrium, non-harmonic system, application of the virial theorem is not strictly appropriate. Ideally, we would instead compute the potential energy contribution from the subsystem, to yield V(t). However, the energy expression used to propagate our MD simulations is not separable, so this is not straightforward. Nevertheless, except at very short time, the energy will be spread between many modes, so that roughly harmonic behaviour should be expected, as well as a reasonable balance between kinetic and potential energy. For this reason, we invoke the virial theorem, because it provides an approximate way to estimate the overall energy of the subsystem as a function of time. We note that an alternative approximation would have been to consider only the kinetic energy of the subsystem—this obviously decays on exactly the same time scale that we obtain here. The resultant correlation function, shown in figure 9, provides information on the rate of energy transfer between the ‘system’ atoms in figure 7 and the ‘bath’ atoms. It shows that energy dissipation is complete within approximately 100 fs, and thus helps to explain the very short time scales observed for the prompt dissociation events in the previously described NVE simulations. It should be noted that the shape of the correlation function will be influenced to a certain extent, at short times, by the averaging over a period τ . With the MD correlation functions in hand, we then fit energy transfer functions F(E ← E ) for use in the EGME. These fits were carried out by varying the parameters in a Gaussian energy transfer function of the form   (E − Eavg )2 , (5.8) F(E ← E ) = ω exp − 2α 2 using a nonlinear least-squares algorithm to minimize the difference between (a) the MD correlation function in figure 9 and (b) the average time-dependent internal energy of SCH3 calculated in MESMER. In equation (5.8), the parameter definitions are as follows: E = E − E;

........................................................

600

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

MESMER energy data

10–1

17

fraction dissociating

10–4

10–5

450

500 550 energy (kJ mol–1)

600

Figure 10. Comparison of the fraction of promptly dissociating CH3 calculated from the NVE MD simulations (crosses) and the EGME (solid line). (Online version in colour.)

Eavg is the most probable quantity of energy transferred; α is the width of the energy transfer distribution; and ω specifies the average frequency (s−1 ) at which energy transfer events occur. In the current simulations Eavg , α and ω were determined to have values of 5.0 × 105 cm−1 , 4.5 × 105 cm−1 and 5.9 × 1014 s−1 , respectively. The approximate nature of the EGME model means that it is not possible to model the finer structure of the MD correlation functions, which arise from the detailed dynamics of coupled system–bath vibrations, and also from the ballistic nature of energy transfer at short times. Having calculated both kd (E) and F(E ← E ), we were then able to solve equation (5.1), and evaluate the fraction of prompt dissociation events predicted by the EGME approach compared to the full NVE simulations. The results, shown in figure 10, were generated by changing the EGME initial conditions vector p(0) in equation (5.6) to correspond to delta-functions centred at the energies in table 4. The results in figure 10 indicate that the statistical EGME model does a reasonable job predicting the competition between reaction and relaxation within ‘system’ atoms. The agreement is particularly good at the lower energies, and slightly worse at higher energies. This is perhaps not an entirely unexpected result, since it is probably the case that non-statistical ballistic effects (which have a more accurate physical treatment in the MD simulations) are more important at higher energies. The agreement between the MD and EGME results, while not perfect, is certainly surprising considering the number of approximations involved in formulating the EGME model, including: (i) the statistical approximation that is the basis for microcanonical transition-state theory, which neglects ballistic dynamical effects, and assumes that energy within the system is distributed ergodically; (ii) the fact that the boundary chosen to divide the system and the bath cuts across covalent bonds within the diamond slab; (iii) the fact that the energy transfer function F(E ← E ) was fitted to a correlation function obtained from linear response theory using MD results from a single equilibrium NVE simulation; (iv) the approximate nature of the energy transfer function F(E ← E ), which is unable to capture the finer dynamical structure of the correlation function obtained from MD simulations; and (v) the assumption that the rate at which the bath removes energy from the system is independent of where precisely that energy is located—i.e. the system is homogeneously coupled to the bath phonon modes.

........................................................

10–3

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

10–2

6. Conclusion

grant no. EP/G00224X/1. Funding for W.J.R. was provided by the EPSRC (grant no. EP/H043292/1), and funding for R.S. was provided by the US Air Force Office of Scientific Research (AFOSR) under contract no. FA9550-16-1-0051.

References 1. Miller JA, Klippenstein SJ. 2006 Master equation methods in gas phase chemical kinetics. J. Phys. Chem. A 110, 10 528–10 544. (doi:10.1021/jp062693x) 2. Vereecken L, Glowacki DR, Pilling MJ. 2015 Theoretical chemical kinetics in tropospheric chemistry: methodologies and applications. Chem. Rev. 115, 4063–4114. (doi:10.1021/ cr500488p)

........................................................

Competing interests. The authors declare that there are no competing interests. Funding. Funding for D.R.G. is from the Royal Society. Funding for J.N.H. was provided by EPSRC Programme

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

In this work, we have carried out a detailed study of etching kinetics at the diamond surface under typical CVD conditions using MD simulations, and also a statistical EGME model. The aim of the work was twofold. First, we wanted to investigate the possibility that surface ‘hotspots’ might lead to prompt CH3 dissociation—i.e. dissociation that occurs prior to the onset of thermal equilibrium. To this end, we developed an efficient reactive PES by fitting an EVB model to higherlevel ab initio electronic structure theory. Using this PES, we carried out 160 000 NVE trajectories, and observed that energy dissipation from surface hotspots on diamond is indeed very fast, but that a small fraction of CH3 does nevertheless undergo dissociation prior to the onset of thermal equilibrium. The results of the MD simulations are in reasonable agreement with experimentally observed etching rates. Second, we wished to investigate the extent to which a significantly cheaper computational model—namely, a statistical EGME model—could capture the NVE trajectory results. To facilitate this aspect of the work, we outlined a general procedure for formulating and solving the EGME for surface chemistry. The idea is to split the system into system and bath components, and then carry out microcanonical transition-state theory in the configuration space of the system atoms. Energy transfer from the system to the bath is estimated using linear response theory from a single long MD trajectory, and fitted to an energy transfer function which can be input into the EGME. Despite the number of approximations involved in formulating this model, the results obtained are in reasonable agreement with the results obtained from 160 000 NVE MD simulations, potentially offering a computationally tractable strategy for investigating nonequilibrium reaction dynamics at surfaces for a broader range of systems beyond the diamond system investigated here. The results obtained in this paper are encouraging insofar as they suggest that it is possible to treat complicated surface dynamics using models which are computationally much less expensive than full MD simulations; however, there are also a range of issues that need to be investigated in further detail moving forward. For example, it would be good to investigate the extent to which the EGME results depend on different methods for drawing a separation between the system and the bath degrees of freedom: How do the results change with differing system/bath boundaries? How sensitive are the – CH3 dissociation probabilities to system size—i.e. how large or small does the system need to be for the EGME results to deviate significantly from the MD results? The form of the energy transfer function also needs further study—e.g. do energy transfer functions that capture the finer structure observed in MD simulations improve the results? And how would such a model perform if it used energy transfer models which explicitly describe energy transfer probabilities in terms of the spectral coupling between the system vibrations and phonon vibrations within the bath? Finally, it would be good to obtain further insight into the energetic regimes over which linear response theory provides an accurate description of energy dissipation.

18

19 ........................................................

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

3. Nesbitt DJ, Field RW. 1996 Vibrational energy flow in highly excited molecules: role of intramolecular vibrational redistribution. J. Phys. Chem. 100, 12 735–12 756. (doi:10.1021/ jp960698w) 4. Glowacki DR, Lockhart J, Blitz MA, Klippenstein SJ, Pilling MJ, Robertson SH, Seakins PW. 2012 Interception of excited vibrational quantum states by O2 in atmospheric association reactions. Science 337, 1066–1069. (doi:10.1126/science.1224106) 5. Dunning G, Glowacki D, Preston T, Greaves S, Greetham G, Clark I, Towrie M, Harvey J, OrrEwing A. 2015 Vibrational relaxation and microsolvation of DF after F-atom reactions in polar solvents. Science 347, 530–533. (doi:10.1126/science.aaa0103) 6. Glowacki DR, Liang C, Marsden SP, Harvey JN, Pilling MJ. 2010 Alkene hydroboration: hot intermediates that react while they are cooling. J. Am. Chem. Soc. 132, 13 621–13 623. (doi:10.1021/ja105100f) 7. Orr-Ewing AJ, Glowacki DR, Greaves SJ, Rose RA. 2011 Chemical reaction dynamics in liquid solutions. J. Phys. Chem. Lett. 2, 1139–1144. (doi:10.1021/jz2002716) 8. Glowacki DR, Orr-Ewing AJ, Harvey JN. 2015 Non-equilibrium reaction and relaxation dynamics in a strongly interacting explicit solvent: F + CD3 CN treated with a parallel multi-state EVB model. J. Chem. Phys. 143, 044120. (doi:10.1063/1.4926996) 9. Bagot PAJ, Waring C, Costen ML, McKendrick KG. 2008 Dynamics of inelastic scattering of OH radicals from reactive and inert liquid surfaces. J. Phys. Chem. C 112, 10 868–10 877. (doi:10.1021/jp8024683) 10. Perkins BG, Häber T, Nesbitt DJ. 2005 Quantum state-resolved energy transfer dynamics at gas−liquid interfaces: IR laser studies of CO2 scattering from perfluorinated liquids. J. Phys. Chem. B 109, 16 396––16 405. (doi:10.1021/jp0511404) 11. Wodtke AM, Matsiev D, Auerbach DJ. 2008 Energy transfer and chemical dynamics at solid surfaces: the special role of charge transfer. Prog. Surf. Sci. 83, 167–214. (doi:10.1016/ j.progsurf.2008.02.001) 12. Sander R. 1999 Modeling atmospheric chemistry: interactions between gas-phase species and liquid cloud/aerosol particles. Surv. Geophys. 20, 1–31. (doi:10.1023/A:1006501706704) 13. Nathanson GM, Davidovits P, Worsnop DR, Kolb CE. 1996 Dynamics and kinetics at the gas−liquid interface. J. Phys. Chem. 100, 13 007–13 020. (doi:10.1021/jp953548e) 14. Donald S, Navin J, Harrison I. 2013 Methane dissociative chemisorption and detailed balance on Pt(111): dynamical constraints and the modest influence of tunneling. J. Chem. Phys. 139, 214707. (doi:10.1063/1.4837697) 15. Navin JK, Donald SB, Harrison I. 2014 Angle-resolved thermal dissociative sticking of light alkanes on Pt(111): transitioning from dynamical to statistical behavior. J. Phys. Chem. C 118, 22 003–22 011. (doi:10.1021/jp505660p) 16. Utz AL. 2009 Mode selective chemistry at surfaces. Curr. Opin. Solid State Mater. Sci. 13, 4–12. (doi:10.1016/j.cossms.2009.01.004) 17. Cavanagh RR, King DS, Stephenson JC, Heinz TF. 1993 Dynamics of nonthermal reactions: femtosecond surface chemistry. J. Phys. Chem. 97, 786–798. (doi:10.1021/j100106a002) 18. Beck RD, Maroni P, Papageorgopoulos DC, Dang TT, Schmid MP, Rizzo TR. 2003 Vibrational mode-specific reaction of methane on a nickel surface. Science 302, 98–100. (doi:10.1126/science.1088996) 19. Villalpando I, John P, Porro S, Wilson JIB. 2011 Hydrogen plasma etching of diamond films deposited on graphite. Diamond Relat. Mater. 20, 711–716. (doi:10.1016/j.diamond.2011. 03.007) 20. Cheesman A, Harvey JN, Ashfold MN. 2008 Studies of carbon incorporation on the diamond {100} surface during chemical vapor deposition using density functional theory. J. Phys. Chem. A 112, 11 436–11 448. (doi:10.1021/jp8034538) 21. Glowacki DR, Lightfoot R, Harvey JN. 2013 Non-equilibrium phenomena and molecular reaction dynamics: mode space, energy space and conformer space. Mol. Phys. 111, 631–640. (doi:10.1080/00268976.2013.780100) 22. Jasper AW, Pelzer KM, Miller JA, Kamarchik E, Harding LB, Klippenstein SJ. 2014 Predictive a priori pressure-dependent kinetics. Science 346, 1212–1215. (doi:10.1126/science. 1260856) 23. Glowacki DR, Paci E, Shalashilin DV. 2009 Boxed molecular dynamics: a simple and general technique for accelerating rare event kinetics and mapping free energy in large molecular systems. J. Phys. Chem. B 113, 16 603–16 611. (doi:10.1021/jp9074898)

20 ........................................................

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 375: 20160206

24. Glowacki DR, Paci E, Shalashilin DV. 2011 Boxed molecular dynamics: decorrelation time scales and the kinetic master equation. J. Chem. Theory Comput. 7, 1244–1252. (doi:10.1021/ ct200011e) 25. Frisch MJ et al. 2009 Gaussian 09, Revision A.1. Wallingford, CT: Gaussian, Inc. 26. Werner H-J, Knowles PJ, Knizia G, Manby FR, Schütz M. 2012 Molpro: a general-purpose quantum chemistry program package. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2, 242–253. (doi:10.1002/wcms.82) 27. Linstrom PJ, Mallard WG (eds). 2015 NIST Chemistry WebBook, NIST Standard Reference Database Number 69. Gaithersburg, MD: National Institute of Standards and Technology. (doi:10.18434/T4D303) 28. Cheesman A, Harvey JN, Ashfold MN. 2005 Computational studies of elementary steps relating to boron doping during diamond chemical vapour deposition. Phys. Chem. Chem. Phys. 7, 1121. (doi:10.1039/b418664h) 29. Richley JC, Harvey JN, Ashfold MN. 2009 On the role of carbon radical insertion reactions in the growth of diamond by chemical vapor deposition methods. J. Phys. Chem. A 113, 11 416– 11 422. (doi:10.1021/jp906065v) 30. Richley JC, Harvey JN, Ashfold MN. 2012 Boron incorporation at a diamond surface: a QM/MM study of insertion and migration pathways during chemical vapor deposition. J. Phys. Chem. C 116, 18 300–18 307. (doi:10.1021/jp305773d) 31. Glowacki DR, Orr-Ewing AJ, Harvey JN. 2011 Product energy deposition of CN + alkane H abstraction reactions in gas and solution phases. J. Chem. Phys. 134, 214 508. (doi:10.1063/ 1.3595259) 32. Warshel A, Weiss RM. 1980 An empirical valence bond approach for comparing reactions in solutions and in enzymes. J. Am. Chem. Soc. 102, 6218–6226. (doi:10.1021/ja00540a008) 33. Brooks BR et al. 2009 CHARMM: the biomolecular simulation program. J. Comput. Chem. 30, 1545–1614. (doi:10.1002/jcc.21287) 34. Halgren TA. 1996 Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94. J. Comput. Chem. 17, 490–519. (doi:10.1002/(SICI)1096987X(199604)17:5/63.0.CO;2-P) 35. Glowacki DR, Reed SK, Pilling MJ, Shalashilin DV, Martínez-Núñez E. 2009 Classical, quantum and statistical simulations of vibrationally excited HOSO2 : IVR, dissociation, and implications for OH + SO2 kinetics at high pressures. Phys. Chem. Chem. Phys. 11, 963–974. (doi:10.1039/B816108A) 36. O’Connor M, Paci E, McIntosh-Smith S, Glowacki DR. 2016 Adaptive free energy sampling in multidimensional collective variable space using boxed molecular dynamics. Faraday Discuss. 195, 395–419. (doi:10.1039/C6FD00138F) 37. Rodgers WJ, May PW, Allan NL, Harvey JN. 2015 Three-dimensional kinetic Monte Carlo simulations of diamond chemical vapor deposition. J. Chem. Phys. 142, 214707. (doi:10.1063/1.4921540) 38. Glowacki DR, Liang C-H, Morley C, Pilling MJ, Robertson SH. 2012 MESMER: an opensource master equation solver for multi-energy well reactions. J. Phys. Chem. A 116, 9545–9560. (doi:10.1021/jp3051033) 39. Davies JW, Green NJB, Pilling MJ. 1986 The testing of models for unimolecular decomposition via inverse Laplace transformation of experimental recombination rate data. Chem. Phys. Lett. 126, 373–379. (doi:10.1016/S0009-2614(86)80101-4) 40. Glowacki DR, Rose RA, Greaves SJ, Orr-Ewing AJ, Harvey JN. 2011 Ultrafast energy flow in the wake of solution-phase bimolecular reactions. Nat. Chem. 3, 850–855. (doi:10.1038/ nchem.1154)