Reaction Profiles and Energy Surfaces of Compressed Species - The

Article pubs.acs.org/JPCA

Reaction Profiles and Energy Surfaces of Compressed Species Jacob Spooner,† Brandon Yanciw,† Brandon Wiebe,‡ and Noham Weinberg*,†,‡ †

Department of Chemistry, Simon Fraser University, Burnaby, BC V5A 1S6, Canada Department of Chemistry, University of the Fraser Valley, Abbotsford, BC V2S 7M8, Canada

‡

S Supporting Information *

ABSTRACT: Both experiment and first principles calculations unequivocally indicate that properties of elements and their compounds undergo a tremendous transformation at ultrahigh pressures exceeding 1 Mbar due to the fact that the difference between intra- and intermolecular interactions disappears under such conditions. Yet, even at much milder pressures of 50−300 kbar, when molecules still retain their individual identity, their chemical properties and reactivity change dramatically. Since kinetics and mechanisms of condensed-phase reactions are described in terms of their potential energy (PES) or Gibbs energy (GES) surfaces, chemical effects of high pressure can be assessed through analysis of pressure-induced deformations of GES of solvated reaction systems. We use quantum mechanical and molecular dynamics simulations to construct GES and reaction profiles of compressed species, and analyze how topography of GES changes in response to compression. We also discuss the important role of volume profiles in assessing pressure-induced deformations and show that the high-pressure GES are well described in terms of these volume profiles and the reference zero-pressure GES.

1. INTRODUCTION It is well-known1 and well-documented2 that high pressures have a profound effect on reaction rates over even a relatively modest range of 1−10 kbar. By 100−300 kbar the strength of intermolecular interactions reaches the level of intramolecular forces,3 which unveils completely new patterns of chemical reactivity.4 Reaction profiles and potential energy surfaces are standard theoretical tools in discussing reactivity and reaction mechanisms at ambient conditions.5 They are often augmented by the consideration of volume profiles6 when high pressure reaction mechanisms are analyzed. In most cases, such volume profiles are treated as three-point schematics describing the volumetric properties of the reactant, transition, and product states.2 Their continuous-curve versions have also been explored in the past7 in the context of studies of high pressure effects on the structure and energy of transition states (TS), as it has been argued that high pressure reactions are more appropriately represented by their enthalpy profiles H(x) = U(x) + PV(x) along reaction coordinate x than by the energy profiles U(x), with the expansion term PV(x) obtained by multiplication of pressure P by volume profile V(x). The major limitation of that approach stemmed from the lack of an adequate computational technique to produce the necessary volume profiles. In the present work we address this problem by proposing a method for the calculation of continuous volume profiles based on molecular dynamics (MD) simulations and use the volume profiles thus obtained to construct high pressure free energy profiles for the model system. We then extend this approach to a multidimensional case and discuss free energy surfaces for model and realistic linear triatomic systems X−H−Y (X, Y = H, F, Cl, Br) in compressed media. © 2013 American Chemical Society

2. COMPUTATIONAL DETAILS MD calculations were performed using the GROMACS package.8 The MD trajectories were obtained using leapfrog integration with 1 fs time step in a cubic box with periodic boundary conditions and an interaction cutoff radius of 0.9−1.0 nm. The systems were maintained at a constant pressure and constant temperature conditions using Berendsen temperature and pressure coupling9 with τT = 0.1 ps and τP = 0.5 ps. Pressure effects of a compressed medium are exerted on an embedded reaction system through solute−solvent interactions and are therefore dependent on the nature of these interactions. To simplify the discussion and to retain the focus of the work on the effects of external pressure, we limited ourselves to consideration of nonpolar solvent−solute systems, where the intermolecular interactions are the least complicated. Monoatomic Lennard-Jones (LJ) particles with parameters m = 40 amu, σ = 0.35 nm, ε = 2.0 kJ mol−1 were used as a model solvent. The same particles were used as constituent elements for the model diatomic and triatomic solutes discussed in the subsequent sections. Organic molecules were described by the OPLS force field.10 A 1000-molecule solvent was used for volume profile calculations in section 3, and a 256-molecule solvent was used for the more computationally demanding volume surfaces in section 4. Quantum-mechanical (QM) calculations were performed with the Gaussian 09 suite.11 All optimized minima and TSs Received: October 23, 2013 Revised: December 10, 2013 Published: December 29, 2013 765

dx.doi.org/10.1021/jp410496t | J. Phys. Chem. A 2014, 118, 765−777

The Journal of Physical Chemistry A

Article

Figure 1. Volume profile ΔV(x) of a model diatomic at the reference pressure of 1 kbar (red diamonds). The linear profile predicted by Stearn−Eyring model1b is shown in blue for comparison. The insets, labeled by the values of x, illustrate the evolution of the sizes and shapes of the solvent cavity.

Figure 3. Volume profile ΔV(x) (red diamonds) and length profile L(x) (blue line) for a model collinear triatomic exchange reaction at the reference pressure of 1 kbar.

Figure 4. Volume profile ΔV(L) of a model triatomic at the reference pressure of 1 kbar (red diamonds). The linear profile predicted by Stearn−Eyring model1b is shown in blue for comparison. Its equation and correlation coefficient are shown in the upper right part of the figure.

Figure 2. PES for model homonuclear linear triatomic system A−B−C (A = B = C). The white line is the BEBO reaction path. The contour spacing is 1.5 kJ/mol.

were obtained at B3LYP/6-31++G(d,p) level and verified by frequency analysis. High pressure potential energy surfaces were generated using HF/6-31G, which offers a reasonable compromise between computational cost and quality, adequate to the exploratory nature of the calculations. Graphic renderings of MD trajectories shown in Figures 1, 5, and 7 and similar images included in the Supporting Information were produced using the VMD package.12

3. VOLUME PROFILES We have recently proposed an MD-based approach toward calculation of molecular volumes, within which the volume V(Z) of a solute molecule Z is found as a difference between the volume of pure solvent consisting of N identical molecules S,

Figure 5. Solvent cavities around the reactant state (a) and TS (b) of abstraction of hydrogen atom (white) by bromine (red) from toluene (blue and white wire frame). 766



Article

Table 1. B3LYP/6-31G++G(d,p)-Optimized RCH and RHBr Distances in TSs of Hydrogen Transfer Reactions of Series A and B and the Experimental Values18 of Activation Volume Differences ΔΔV‡ reference series A

series B

substrate

RCH (nm)

RHBr (nm)

ΔΔV (cm3/mol)

toluene m-xylene p-xylene m-bromotoluene p-chlorotoluene ethylbenzene cumene diphenylmethane

0.1616 0.1619 0.1628 0.1589 0.1601 0.1579 0.1515 0.1423

0.1534 0.1535 0.1533 0.1540 0.1540 0.1546 0.1570 0.1624

0 1.5 3.0 −3.0 −1.2 −4.8 −5.5 −2.4

isomerization process between two states of the diatomic represented by two stable minima of potential (eq 4) located at x = 0.15 nm (“short state”) and x = 0.35 nm (“long state”). The minima are separated by a barrier of 9 kJ/mol with a TS at x = 0.25 nm. MD simulations were performed at T = 300 K and P = 1 kbar, under which conditions the solvent was liquid.14 As can be seen from the insets in Figure 1, the solvent cavities retain nearly constant cross-sections σ for x ranging between 0.1 and 0.3 nm. As a result, the volume profile for this system is well described by the linear equation ΔV(x) = σ(x − x0) of the Stearn−Eyring model.1b The radius of the circular cross-section estimated from this linear dependence is 0.31 nm. At x > 0.3 nm, a narrowing develops in the middle of the cavity, resulting in a marked nonlinearity of the volume profile. Linear Triatomic System. The volume profile for the bistable diatomic is predictably a monotonic function of x since this variable describes the size of the system. The same can be expected for other simple processes, such as dissociation/ association,15 where a single distance parameter serves as a reaction coordinate for the duration of the entire process. The situation is quite different in most other cases, where the reaction coordinate involves more than one internal coordinate of a system. Substitution reactions constitute an important case; collinear exchange reactions A + BC → AB + C being a classic example of this type. A linear triatomic system ABC that undergoes this transformation is described by two geometrical parameters: A−B distance RAB and B−C distance RBC. A typical potential energy surface16 (PES) for a homonuclear case A = B = C is shown in Figure 2 (the numerical values of the potential energy function U(RAB,RBC) determining this PES can be found in the Supporting Information). The reaction path can be approximated reasonably well by the Bond Energy Bond Order17 (BEBO) parametric equation

V(NS), and the volume V(Z + NS) of the solvent−solute system:13 V (Z) = V (Z + N S) − V (N S)

(1) ‡

The method was used to calculate volumes of activation ΔV as a difference ΔV ‡ = V (Z‡) − V (Z0)

(2) ‡

between the volumes of reaction system Z in TS Z and reactant state Z0.This approach can easily be extended to calculation of continuous volume profiles ΔV(x) if an arbitrary transient structure Z(x) along reaction coordinate x is used instead of the TS structure Z‡ in eq 2: ΔV (x) = V (Z(x)) − V (Z0)

(3)

Bistable Diatomic. An example of a volume profile obtained using the above procedure13c is shown in Figure 1 for a model bistable diatomic described by a quartic double-well potential 1 U (x) = λx 2(1 − γx 2) (4) 2 −1 −2 −2 with λ = −3600 kJ mol nm and γ = 50 nm . The variable bond length x serves as a reaction coordinate that describes the

⎧ RAB = α − β ln(n) ⎪ ⎨ ⎪ ⎩ RBC = α − β ln(1 − n)

(5)

Figure 6. TSs for reactions of series A and B plotted in coordinates (RCH, RHBr) together with BEBO reaction path (eq 8). 767



Article

Figure 7. Plot of experimental activation volume differences,17 ΔΔV‡, vs the differences in C−Br distances, ΔL‡, in TSs of hydrogen abstraction by bromine radical from substituted toluenes: blue circles, series A; red squares, series B. The insets show solvent cavity cross-sections orthogonal to C−H−Br line and passing through H atom.

with α = 0.093 nm and β = 0.028 nm. The bond order parameter n can serve as a nonuniformly scaled reaction coordinate with values of n = 0, 0.5, and 1 for reactant, TS, and product, respectively. Alternatively, a uniformly scaled reaction coordinate x (with the reference x = 0 for TS) can be defined as a signed length of the reaction path n

x(n) = β

∫0.5

η2 + (1 − η)2 η(1 − η)

dη

cross-section in this case is 0.30 nm, which is consistent with the value found for the bistable diatomic. Stearn−Eyring Model. It appears that the Stearn−Eyring model ΔV (x) = σ[L(x) − L(x0)]

(7)

works reasonably well for different types of volume profiles as long as the cavity cross-section σ remains constant. Successful application of the model to various types of chemical reactions1 seems to support that. It must be emphasized that the requirement of σ = const needs to be satisfied only for the part of the cavity adjacent to the reaction center. The rest of the system does not change in the course of reaction and its volume remains constant, thus giving zero contribution to the volume profile. As an illustration, we will consider the case of two reaction series of hydrogen atom abstraction by bromine radical from the aliphatic hydrogen of substituted toluenes:18

(6)

The volume profile along coordinate x calculated for the above triatomic system at T = 300 K and P = 1 kbar is shown in Figure 3. As could be expected, the graph represents a curve with a minimum at x = 0 corresponding to the TS. Remarkably, the volume profile matches closely the linear dimension of the triatomic L(x) = RAB + RBC. The graph becomes linear if, instead of reaction coordinate x, the volume profile is plotted against parameter L (Figure 4), as it was originally suggested by Stearn and Eyring.1b The estimated radius of the circular

X−Ph−CH 2H + Br • → X−Ph−CH 2• + HBr 768

(series A)



Article

Figure 8. Exact and approximate energy profiles for a bistable diatomic (left) and a model collinear triatomic exchange reaction (right): magenta, gas-phase energy profile U(x), eq 4; green, reference Gibbs energy profile ΔG0(x) at P0 = 1 kbar; red, ΔG(x;P) at P = 1.5 kbar; blue, ΔG(x;P) at P = 3 kbar; solid lines, exact ΔG(x;P); open circles, approximate ΔG(x;P), eq 12.

is indeed observed for series A where cavity cross-sections are roughly equal due to the fact that substituents X are removed sufficiently far from the reaction center. In the case of series B, the substituents are introduced directly into the reaction center, which results in a substantial variation of σ and, consequently, in the loss of linearity between ΔΔV‡ and ΔL‡.

and Ph−CY′Y″H + Br • → Ph−CY′Y″• + HBr

(series B)

with X = H, p-CH3, m-CH3, p-Cl, m-Br, and Y′Y″ = HH, CH3H, CH3CH3, PhH. The images of the solvent cavities around the reactant state and TS of the unsubstituted toluene (X = H, Y′Y″ = HH) are shown in Figure 5. The major change in the cavity occurs in the course of reaction in the region adjacent to the reaction center C−H−Br (the “neck” of the rubber-duck-shaped cavity). The situation is similar in other reactions of series A and B (see Supporting Information). In all cases, the C−Br distance serves as the linear dimension parameter, L = RCBr = RCH + RHBr. The values of RCH and RHBr for the B3LYP/6-31++G(d,p)-optimized19 TSs of reactions of both series are listed in Table 1 along with the experimental values18 of activation volume differences ΔΔV‡ = ΔV‡ − ΔV‡ref determined with reference to the reaction of unsubstituted toluene. As can be seen from Figure 6, these TSs cluster fairly close to the common BEBO reaction path ⎧ RHBr = 0.142 − 0.039 ln(n) ⎪ ⎨ ⎪ ⎩ R CH = 0.108 − 0.039 ln(1 − n)

4. REACTION PROFILES AND ENERGY SURFACES FOR MODEL SYSTEMS Reaction Profiles. Volume is an isothermal pressure derivative of Gibbs energy. Therefore, the Gibbs energy profile ΔG(x;P) at a given pressure P can be obtained from the Gibbs energy profile ΔG0(x) = ΔG(x;P0) at a reference pressure P0 by integration of the pressure-dependent volume profile ΔV(x;P): ΔG(x ; P) = ΔG0(x) +

ΔV (x ; p)dp

(11)

If the pressure dependence of ΔV(x) can be neglected, eq 11 reduces to a much simpler approximate relationship ΔG(x ; P) ≈ ΔG0(x) + (P − P0)ΔV (x)

(8)

(12)

We used these equations to evaluate the effect of pressure on the shapes of reaction profiles for the short-to-long isomerization of a bistable diatomic and model collinear triatomic exchange reaction described in the previous section. The necessary reference Gibbs energy profiles ΔG0(x) were obtained from the gas-phase energy profiles U(x) by adding to them MD solvation Gibbs energies acquired by thermodynamic integration21 at a reference pressure P0 = 1 kbar. The results are shown in Figure 8. The approximate pressure-dependent profiles ΔG(x;P) calculated using eq 12 compare favorably with the exact profiles obtained directly from U(x) using thermodynamic integration at the target pressure P. Despite neglecting pressure dependence of the volume profile, approximate eq 12 works quite well. This should probably be attributed to the fact that the major pressureinduced restructuring of the solute happens along its softest coordinate x, which is considered explicitly in the equation. This restructuring is substantial: both the barrier height and

ΔΔV ‡ = ΔV (x‡) − ΔV (x‡ ref ) (9)

where x‡ is the reaction coordinate of a TS, and x0 is assumed to be the same for all reactions. If the cavity cross-section σ for a given reaction is close to that of the reference reaction, σref, then eq 9 simplifies to ΔΔV ‡ = σ[L(x‡) − L(x‡ ref )] = σ ΔL‡

P

0

This allows one to use the same reaction coordinate x for all reactions of these series. If the Stearn−Eyring model applies, the value of ΔΔV‡ can be expressed in terms of eq 7 as = σ[L(x‡) − L(x0)] − σref [L(x‡ ref ) − L(x0)]

∫P

(10)

That is, if σ is constant in a reaction series, the relative values of activation volumes, ΔΔV‡, can be expected to be proportionate to the relative values of C−Br distances ΔL‡ in the TSs. As can be seen in Figure 7, such linear correlation 769



Article

Figure 10. Linear X−H−X system in a 1D Ne matrix: X is green, H is white, and Ne is light blue.

Figure 11. Linear X−H−X system in a 3D Ne matrix: X is green, H is white, and Ne is light blue.

⎧⎛ dΔG(x ; P) ⎞ ⎛ d2ΔG(x ; P) ⎞⎫ aP − a0 ≈ −⎨⎜ ⎟/⎜ ⎟⎬ ⎠ ⎝ dx dx 2 ⎠⎭ ⎩⎝ x=a ⎪

⎪

⎪

⎪

(13)

0

By virtue of eq 12 and since the first derivative of ΔG0(x) at its stationary point a0 is zero, eq 13 reduces to ⎧⎛ dΔV (x) ⎞ ⎛ d2ΔG(x ; P) ⎞⎫ aP − a0 ≈ −(P − P0)⎨⎜ ⎟/⎜ ⎟⎬ dx 2 ⎠⎭ ⎩⎝ dx ⎠ ⎝ x=a ⎪

⎪

⎪

⎪

= −(P − P0)ΔV ′(a0)/ΔG″(a0 ; P)

0

(14)

For reactant and product minima, the second derivative ΔG″(a0;P) is positive. Therefore, for the bistable diatomic, where the volume profile ΔV(x) is an increasing function of x, both short- and long-state minima shift to the left with pressure. On the contrary, since ΔG″(a0;P) is positive for the TS it shifts to the right. In the case of the triatomic exchange reaction, ΔV′(x) is negative for the reactant, positive for the product, and zero for the TS. As a result, the reactant minimum moves to the right, the product shifts to the left, and the TS does not change its position with pressure. Energy Surfaces. The bistable diatomic is a one-dimensional (1D) system and its energy profile U(x) constitutes its full potential energy surface (PES). The situation is different for multidimensional systems, such as the linear triatomic system where the energy profile U(x) is just a cross-section of its 2D PES. Solvation modifies this PES and, consequently, the reaction path, which redefines reaction coordinate x and energy profile U(x). Compression is expected to further magnify this effect. Therefore, above a certain threshold, the effect of pressure on such systems is more appropriately described by its effect on their PESs rather than reaction profiles. Equations 11 and 12 can easily be upgraded to a multidimensional case by simple replacement of reaction coordinate x with a multidimensional coordinate vector x:

Figure 9. Reference 1 kbar GES ΔG0(x) (a) and the exact (b) and approximate (c) pressure-dependent GES ΔG(x;P) of the model linear triatomic system at P = 3 kbar. The contour spacing is 2 kJ/mol.

positions of stationary points (minima and TS) are affected by pressure. The direction of these shifts can be assessed using eq 12.7b If x = a0 is a stationary point of ΔG0(x) and x = aP is the corresponding stationary point of ΔG(x;P), their relative shift aP − a0 can be estimated in the first iteration of the Newton−Raphson optimization22 as

Table 2. Geometrical Parameters of Stationary Points and Reaction Barriers of Three High Pressure GESs for the Collinear Triatomic Exchange Reaction reactant

TS

pressure (kbar)

energy surface

RAB (nm)

RBC (nm)

RAB (nm)

RBC (nm)

reaction barrier (kJ/mol)

1.0 3.0

reference GES (Figure 9a) exact GES (Figure 9b) approximate GES (Figure 9c)

0.0916 0.0912 0.0914

0.1897 0.1737 0.1761

0.1094 0.1086 0.1085

0.1094 0.1086 0.1085

14.9 13.4 13.1

770



Article

Table 3. Parameters of the HF/6-31G PES of Symmetric X−H−X Systems (X = F, Cl, Br) Embedded in Ne Matrices of Varying Number of Shells, with the Outer Shell Frozen; Parameters of the Gas-Phase PES Are Also Included for Comparison reactant system F−H−F

compression zero low

medium

high

Cl−H−Cl

zero low

medium

high

Br−H−Br

zero low medium high

TS

matrix size

RXH (nm)

RHX (nm)

RXH (nm)

gas phase 1-shell 2-shell 3-shell 1-shell 2-shell 3-shell 1-shell 2-shell 3-shell gas phase 1-shell 2-shell 3-shell 1-shell 2-shell 3-shell 1-shell 2-shell 3-shell gas phase 1-shell 2-shell 1-shell 2-shell 1-shell 2-shell 3-shell

0.092 0.091 0.091 0.091 0.089 0.090 0.090 0.087 0.087 0.088 0.130 0.129 0.129 0.129 0.128 0.128 0.128 0.133 0.133 0.135 0.141 0.141 0.141 0.143 0.143 0.145 0.144 0.144

0.214 0.189 0.189 0.189 0.140 0.142 0.142 0.111 0.113 0.118 0.301 0.228 0.229 0.229 0.169 0.169 0.169 0.133 0.133 0.135 0.280 0.223 0.230 0.170 0.173 0.145 0.144 0.144

0.113 0.110 0.112 0.113 0.107 0.106 0.105 0.098 0.098 0.100 0.154 0.153 0.153 0.150 0.145 0.144 0.145 single minimum single minimum single minimum 0.167 0.165 0.165 0.155 0.155 single minimum single minimum single minimum

reaction barrier (kJ/mol)

no no no

no no no

152 130 124 124 65 65 65 14 17 23 82 58 56 54 9 8 8 reaction reaction reaction 56 38 27 2 2 reaction reaction reaction

Table 4. Parameters of the HF/6-31G PES of Asymmetric X−H-Y Systems (X, Y = F, Cl, Br) Embedded in the 2-Shell Ne Matrix with the Outer Shell Frozen; Parameters of the Gas-Phase PES Are Also Included for Comparison XH + Y

X + HY

system

compression

RXH (nm)

RHY (nm)

RXH (nm)

F−H−Cl

gas phase low medium high gas phase low medium high gas phase low medium high

0.092 0.092 0.090 0.089 0.092 0.092 0.089 0.089 0.130 0.129 0.128 0.130

0.275 0.221 0.177 0.148 0.261 0.224 0.158 0.152 0.272 0.241 0.180 0.147

0.248 0.190 0.149

F−H−Br

Cl−H−Br

ΔG(x ; P) = ΔG0(x) +

∫P

0.246 0.193

0.295 0.231

0.130 0.128 0.125 single 0.142 0.140 single single 0.142 0.141 single single

RXH (nm) 0.120 0.118 0.116 minimum 0.124 0.120 minimum minimum 0.159 0.155 minimum minimum

reaction barrier (kJ/mol) RHY (nm) 0.148 0.147 0.138 0.159 0.158

0.163 0.162

XH + Y 164 137 77 no 165 140 no no 80 65 no no

X + HY 95 72 19 reaction 78 53 reaction reaction 60 40 reaction reaction

at P = 3 kbar. The reference-pressure GES ΔG0(x) and the exact GES ΔG(x;P) were obtained from the gas-phase PES U(x) of Figure 2 by augmenting them with MD solvation Gibbs energies generated by thermodynamic integration at P0 = 1 kbar and P = 3 kbar, respectively. The results are shown in Figure 9. As can be seen from the figure and Table 2 listing parameters of the GESs, the exact and approximate surfaces are fairly similar. As in the case of Gibbs energy profiles, reaction barriers and stationary points of the GESs are strongly affected by pressure. The pressure-induced displacement of the stationary points is

P

ΔV (x ; p)dp

0

≈ ΔG0(x) + (P − P0)ΔV (x)

TS

RHY (nm)

(15)

ΔV(x;P) and ΔV(x) thus become multivariable functions and should therefore more appropriately be referred to as volume surfaces rather than volume profiles. We calculated the ΔV(x) surface for the model linear triatomic for P0 = 1 kbar and used it, in conjunction with eq 15, to construct an approximate pressure-dependent Gibbs energy surface (GES) for this system 771



Article

Figure 12. PES for FHF system at ambient conditions (a), low compression (b), medium compression (c), and high compression (d). The corresponding TS positions (marked by symbols) shift in the diagonal direction with increasing pressure. The arrows project TS positions between lower and higher compression PESs. The contour spacing is 10 kJ/mol.

roughly described by eq 16, generalizing eq 14 to a multidimensional case: ⎧⎛ ∂ΔV (x) ⎞ ⎛ ∂ 2ΔG(x ; P) ⎞⎫ ⎟⎬ λiP − λi0 ≈ −(P − P0)⎨⎜ ⎟/⎜ ∂λi 2 ⎠⎭ ⎩⎝ ∂λi ⎠ ⎝ x=a ⎪

⎪

⎪

⎪

0

(16)

where λiP and λi0 are normal coordinates at the stationary points a and a0 of ΔG(x;P) and ΔG0(x). The normal modes for the TS are λ1 = L and λ2 = x. Since ∂ΔV/∂x = 0, x‡ = x‡0, i.e., the TS does not shift along the reaction coordinate in the multidimensional case either. At the same time, ∂ΔV/∂L > 0 and ∂2ΔG/∂L2 > 0, and as a result, the TS is somewhat shifted with pressure in the diagonal direction toward the lower left corner. For the reactant and product minima, the normal modes are λ1 ≈ RAB and λ2 ≈ RBC, and the minima shift in both directions due to the fact that ∂ΔV/∂λi > 0 and ∂2ΔG/∂λi2 > 0. Since the second derivative of ΔG(x;P) along the reaction coordinate x is smaller than that across, both reactant and product minima are shifted toward the TS to a much greater degree than in the lateral direction. This lateral deformation of GES can likely be ignored at lower pressures, and the reaction path remains practically unchanged, which makes volume and reaction profiles a meaningful tool for studying high pressure effects over the 0−10 kbar range. However, as we show below, the situation becomes dramatically different as pressure is further increased.

Figure 13. High compression PES for the Cl−H−Cl system. The contour spacing is 10 kJ/mol.

5. ENERGY SURFACES FOR HYDROGEN TRANSFER REACTIONS AT EXTREME PRESSURES A close realistic analog of the model triatomic system discussed in the previous sections is a hydrogen transfer reaction X + HY → XH + Y between halogen atoms (X and Y = F, Cl, Br). We modeled the effect of extreme pressures on their PESs by 772



Article

Figure 14. Plots of the minima and TS positions for X−H−Y systems at ambient conditions (purple), and low (blue), medium (green), and high (black) compression. The dashed lines are BEBO reaction paths.

embedding the reaction systems into compressed 1D and 3D neon matrices as shown in Figures 10 and 11. The PES were calculated by scanning RXH and RHY on a 25 × 25 grid with a step of 0.005 nm. One-Dimensional Matrix. Linear triatomic systems X−H−Y were inserted coaxially within a chain of Ne atoms (one to three Ne atoms on each side), and positions of the terminal Ne atoms were adjusted to attain the desired degree of compression, while the geometry of the rest of the system was optimized. The PESs were then constructed by scanning the RXH and RHY parameters of X−H−Y with positions of the terminal Ne atoms held fixed and the rest of parameters of the

system optimized. The force exerted onto the system was calculated as a negative derivative F = −dU/dL of the optimized energy U(L) with respect to distance L between terminal Ne atoms. We considered three regimes, referred to as low, medium, and high compression, and corresponding to compression forces of 0.15, 0.75, and 2.0 hartree/nm, respectively. These forces were estimated to be roughly equivalent to pressures of 150, 900, and 3000 kbar in the case of the F−H−F system in 3D neon matrix. Varying the number of Ne atoms in the solvent matrix did not seem to have a significant effect on the energy surfaces (see Table 3). The parameters of the stationary points obtained for symmetric and asymmetric X−H−Y systems are summarized in 773



Article

Tables 3 and 4. In all cases, pressure has a tremendous effect on the PES topography. With increasing pressure, the reactant and product minima shift toward the TS, as predicted by eq 16, and the reaction barrier gets lower. The trend is illustrated in Figure 12 comparing PESs for the F−H−F system. At sufficiently high compression the systems collapse into a single minimum, like in the case of the high-pressure PES for Cl−H− Cl shown in Figure 13. The full set of PESs for all X−H−Y systems at various compressions is included in the Supporting Information. Figure 14 displays the pressure-dependent behavior of the stationary points and reaction pathways for these systems. At lower compressions, the reaction paths are fairly similar to the zero-pressure paths despite the notable shift of reactant and product minima along the reaction coordinate. However, as pressure increases and deformation of the PES becomes more significant, the reaction paths shift in the diagonal direction toward the lower left corner. This behavior is common for both symmetric and asymmetric systems, with one notable difference; whereas the single minimum of a highly compressed PES for a symmetric system appears as the limit of the TS sequence, for the asymmetric system it appears as the limit for the sequence of the exothermic minima. The pressure dependence of the reaction barriers was roughly exponential (see Figure 15).

Figure 16. Correlation between reaction barriers EXHY for X−H−Y systems obtained by direct QM calculations and from the Marcus relation: red, zero compression; purple, low compression; green, medium compression.

Figure 17. Correlation between activation energies EXHY and the squares of the H−X bond elongation (R‡HX−RHX,e)2 in the TSs of F−H−F (blue circles), Cl−H−Cl (green squares), and Br−H−Br (red triangles).

reduces to the enthalpy function ΔH(x;P). In addition, for a 1D model, the expansion term PΔV(x) becomes FΔL(x), where L(x) is the linear size of the system and F is the compression force, playing the roles of 1D volume and 1D pressure, respectively. By virtue of the linearity of the X−Y−H system, L(RXH,RHY) = RXH + RHY, and with T = 0 and P0 = 0, eq 15 turns into

Figure 15. Dependence of the reaction barriers EXHY for X−H−Y systems on the strength of compression force F.

The calculated activation energies EXHY(P) obeyed the Marcus relation23 EXHY = (EXHX + EYHY)/(2)(1 + 1/2(EXHY − EYHX)/ (EXHX + EYHY))2 (see Figure 16), despite the fact that of all symmetric X−H−X cases, only EFHF correlated well with the square (R‡HX − RHX,e)2 of the H−X bond elongation in the TS relative to the equilibrium bond length (see Figure 17), as it could be expected on the basis of a simple two-parabola model. Approximate Energy Surfaces. Since approximate eq 15 offers a considerable computational saving in comparison to the full quantum-mechanical calculations, it is important to know if the description of high-pressure PESs provided by this equation is sufficiently accurate. Since all calculations in this section refer to zero temperature, the Gibbs energy function ΔG(x;P)

H(RXH , RHY ; F ) ≈ U0(RXH , RHY ) + FL(RXH , RHY ) = U0 + F(RXH + RHY )

(17)

The approximate high-pressure PESs H(RXH,RHY;F) were derived from the zero-pressure PES U0(RXH,RHY) using eq 17 and showed a remarkable similarity to the exact PES described in the previous subsection. Parameters of the approximate PESs for symmetric X−H−X systems are listed in Table 5 in comparison with the corresponding parameters of the exact PESs. Three-Dimensional Matrix. The 1D model of Figure 10 was adequate to describe the effects of the longitudinal compression 774



Article

Table 5. Parameters of the Approximate PES of Symmetric X-H-X Systems (X = F, Cl, Br) Obtained Using eq 17; the Figures in Parentheses Represent Their Relative Deviation from the Corresponding Parameters of the Exact PES Listed in Table 3 reactant system F−H−F

compression low medium high

Cl−H−Cl

low medium high

Br−H−Br

low medium high

average deviation

TS

RXH (nm)

RHX (nm)

RXH (nm)

0.091 (0%) 0.089 (1%) 0.088 (1%) 0.129 (0%) 0.128 (0%) 0.133 (0%) 0.141 (0%) 0.142 (1%) 0.143 (1%) 0.5%

0.185 (2%) 0.144 (1%) 0.117 (4%) 0.231 (1%) 0.170 (1%) 0.133 (0%) 0.231 (0%) 0.173 (0%) 0.143 (1%) 1%

0.110 (2%) 0.106 (0%) 0.110 (12%) 0.151 (1%) 0.143 (1%) single min (0%) 0.164 (1%) 0.154 (1%) single min (1%) 2%

reaction barrier, (kJ/mol) 115 (7%) 63 (3%) 20 (18%) 50 (11%) 8 (0%) no reaction (0%) 33 (22%) 2 (0%) no reaction (0%) 7%

Figure 18. Pressure dependence of the lengths of long (red squares) and short (green circles) HH bonds in the H−H−H− system.

Softer systems with lower activation barriers undergo a transition to a single-minimum state at lower pressures as evidenced by the comparison of hydrogen transfer reactions involving fluorine, chlorine, and bromine atoms. Unsaturated systems represent another example of softer species. It has been argued earlier7b that the reaction barrier for Diels−Alder dimerization of cyclopentadiene gets reduced from about 70 kJ/mol at ambient conditions to 20 kJ/mol at 50 kbar. If a high-pressure single minimum replaces the barrier of a thermoneutral process, as in the case of symmetric X−H−X systems, the high pressure species associated with that minimum are unstable at low pressures when the barrier is restored. The transition between high- and low-pressure forms is therefore expected to be reversible in such cases, as it was with the H3− anion in polyhydrides.24 On the contrary, for exothermic processes, where high pressures simultaneously eliminate the barrier and the endothermic minimum, the high-pressure state is associated with the exothermic minimum and therefore remains stable at lower pressures. An example of that kind is high pressure polymerization of aromatic compounds,26 where static pressure of 100−350 kbar eliminates the barrier separating the endothermic minimum of the monomer state and the exothermic minimum of the polymer state. Approximate expressions 12 and 15, together with the compressed matrix approach used in this work, promise remarkable savings and computational efficiency in producing highpressure reaction profiles and energy surfaces, the knowledge of which is crucial to understanding reactivity and spectral properties of compressed species.

but left the lateral compression of a reaction system outside of its scope. It was also unsuitable for calculation of the pressure equivalent of a given compression force. To address these issues, we considered a 3D matrix built by encasing the linear model in a rectangular 3 × 3 × 9 prism of uniformly spaced Ne atoms as shown in Figure 11. As in the case of the 1D matrix, the compressing forces were calculated as negative partial derivatives of energy with respect to the linear dimensions of the prism: Fi = −dU/dLi (i = x,y,z). The pressure equivalents were found by dividing the forces by the appropriate cross-section areas to which the forces were applied. Both isotropic (Px = Py = Pz) and axially anisotropic (Px = Py ≠ Pz) compressions were considered for the symmetric F−H−F system. The resulting PESs (see Supporting Information) were very similar to those shown in Figure 12 for the 1D case. The longitudinal compression had a profound effect on the PES topography, whereas the effect of the lateral compression was rather marginal.

6. CONCLUDING REMARKS The reported results indicate that reaction profiles and energy surfaces of compressed systems can undergo a tremendous change with pressure, sometimes resulting in the complete annihilation of the reaction barrier separating reactants and products. This is in accord with the findings of the recent studies of potassium, rubidium, cesium, and barium polyhidrides,24 which observed formation of a symmetric linear H3− anion at elevated pressures. At lower pressures this ion disproportionated into a linear complex H−···H2, stable at ambient conditions.25 The exact pressure of this transition depended on the nature of the metal cation and crystal stoichiometry, and for RbH5, it was about 300 kbar.24b To check the accuracy of the predictions based on our system-in-matrix calculations, we optimized the geometry of the linear H−H−H− system embedded in a compressed 3D Ne matrix, as shown in Figure 11 for F−H−F. The calculations were performed with PBE/6-311+ +G(d,p) to match the level of theory used in refs 24 and 25b. The results presented in Figure 18 predict a transition to the symmetric H3− anion above 300 kbar.

■

ASSOCIATED CONTENT

S Supporting Information *

Numerical values of the potential energy surface for the model triatomic; images of the solvent cavities around the reactant state and TS of the substituted toluene systems; compressed potential energy surfaces for all X−H−Y systems; potential energy surfaces for the F−H−F system embedded in a 3D Ne matrix; full citation for ref 11. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*(N.W.) E-mail: [email protected]. 775



Article

Notes

(8) (a) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS-A Message Passing Parallel Molecular-Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43−56. (b) Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: A Package for Molecular Simulation and Trajectory Analysis. J. Mol. Model. 2001, 7, 306−317. (c) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (d) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Rfficient, LoadBalanced, and Scalable Molecular Simulation. J. Chem. Theory Comut. 2008, 4, 435−447. (9) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics With Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (10) (a) Jorgensen, W. L.; Tirado-Rives, J. The OPLS Potential Functions for Proteins: Energy Minimizations for Crystals of CyclicPeptides and Crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (b) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225−11236. (11) Frisch, M. J.; et al. Gaussian 09, revision A.1; Gaussian, Inc.: Wallingford, CT, 2009 (see full reference in Supporting Information). (12) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (13) (a) Deglint, E.; Martens, H.; Edwards, E.; Boon, N.; Dance, P.; Weinberg, N. Molecular Dynamics Calculations of Activation Volumes. Phys. Chem. Chem. Phys. 2011, 13, 438−440. (b) Wiebe, H.; Spooner, J.; Boon, N.; Deglint, E.; Edwards, E.; Dance, P.; Weinberg, N. Calculation of Molecular Volumes and Volumes of Activation Using Molecular Dynamics Simulations. J. Phys. Chem. C 2012, 116, 2240−2245. (c) Spooner, J.; Wiebe, H.; Boon, N.; Deglint, E.; Edwards, E.; Yanciw, B.; Patton, B.; Thiele, L.; Dance, P.; Weinberg, N. Molecular Dynamics Calculation of Molecular Volumes and Volumes of Activation. Phys. Chem. Chem. Phys. 2012, 14, 2264− 2277. (14) Nicolas, J.; Gubbins, K.; Streett, W.; Tildesley, D. Equation of State for the Lennard-Jones Fluid. Mol. Phys. 1979, 37, 1429−1454. (15) (a) Crick, I. S.; Tregloan, P. A. Temperature and Pressure Dependent Equilibrium Studies of 1,4,8,11-Tetramethyl-1,4,8,11tetraazacyclotetradecane Nickel(II) in Coordinating Solvents. Inorg. Chim. Acta 1988, 142, 291−299. (b) Zhang, M.; van Eldik, R.; Espenson, J. H.; Bakac, A. Volume Profiles for the Reversible Binding of Dioxygen to Cobalt(II) Complexes: Evidence for a SubstitutionControlled Process. Inorg. Chem. 1994, 33, 130−133. (16) Atkins, P. W.; De Paula, J. Physical Chemistry, 9th ed.; Oxford University Press: New York, 2010. (17) Johnston, H. S.; Gas Phase Reaction Rate Theory; Ronald Press Co.: New York, 1966. (18) (a) Zhulin, V. M.; Botnikov, M. Y.; Milyavskaya, I. K. Linear Correlation Between Constants of Substituents and Volume Effects of Activation. Russ. Chem. Bull. 1975, 24, 424. (b) Zhulin, V. M.; Botnikov, M. Y.; Milyavskaya, I. K. Bromination of Hydrocarbons by N-Bromosuccinimide at High Pressures. Russ. Chem. Bull. 1975, 24, 1038−1041. (c) Zhulin, V. M.; Botnikov, M. Y.; Milyavskaya, I. K. Position of the Enol Proton in the Chelate Forms of Unsymmetrical βdiketones. Russ. Chem. Bull. 1977, 26, 517−521. (19) Although somewhat inadequate for open-shell systems, B3LYP/ 6-31G++G(d,p) calculations produce good geometry and decent energy parameters for the TS of hydrogen transfer from the methyl group.20 (20) (a) Zheng, J.; Zhao, Y.; Truhlar, D. G. Thermochemical Kinetics of Hydrogen-Atom Transfers Between Methyl, Methane, Ethynyl, Ethyne, and Hydrogen. J. Phys. Chem. A 2007, 111, 4632−4642. (b) Vandeputte, A. G.; Sabbe, M. K.; Reyniers, M.-F.; van Speybroeck, V.; Waroquier, M.; Marin, G. B. Theoretical Study of the Thermodynamics and Kinetics of Hydrogen Abstractions from Hydrocarbons. J. Phys. Chem. A 2007, 111, 11771−11786.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support by the Natural Sciences and Engineering Research Council of Canada (NSERC) and the University of the Fraser Valley is gratefully acknowledged. This research has been enabled by the use of WestGrid computing resources, which are funded in part by the Canada Foundation for Innovation, Alberta Innovation and Science, BC Advanced Education, and the participating research institutions. WestGrid equipment is provided by IBM, Hewlett-Packard, and SGI.

■

REFERENCES

(1) (a) Evans, M. J.; Polanyi, M. Some Applications of the Transition State Method to the Calculation of Reaction Velocities, Especially in Solution. Trans. Faraday Soc. 1935, 31, 875−893. (b) Stearn, A. E.; Eyring, H. Pressure and Rate Processes. Chem. Rev. 1941, 29, 509− 523. (c) Hamann, S. D. Physico-Chemical Effects of Pressure; Academic Press: New York, 1957. (d) Gonikberg, M. G. Chemical Equilibria and Reaction Rates at High Pressures; National Science Foundation: Washington, DC, 1963. (e) Isaacs, N. S. Liquid Phase High Pressure Chemsitry; Wiley: New York, 1981. (2) (a) le Noble, W. J. Kinetics of Reactions in Solutions Under Pressure. Prog. Org. Chem. 1967, 5, 207−330. (b) Asano, T.; le Noble, W. J. Activation Volumes and Reaction Volumes in Solution. Chem. Rev. 1978, 78, 407−489. (c) van Eldik, R.; Asano, T.; le Noble, W. J. Activation Volumes and Reaction Volumes in Solution 2. Chem. Rev. 1989, 89, 549−688. (d) Drljaca, A.; Hubbard, C. D.; van Eldik, R.; Asano, T.; Basilevsky, M. V.; le Noble, W. J. Activation Volumes and Reaction Volumes in Solution 3. Chem. Rev. 1998, 98, 2167−2289. (3) Grochala, W.; Hoffmann, R.; Feng, J.; Ashcroft, N. W. The Chemical Imagination at Work in Very Tight Places. Angew. Chem., Int. Ed. 2007, 46, 3620−3642. (4) (a) Hemley, R. J. Effects of High Pressure on Molecules. Annu. Rev. Phys. Chem. 2000, 51, 763−800. (b) Bini, R. Laser-Assisted HighPressure Chemical Reactions. Acc. Chem. Res. 2004, 37, 95−101. (c) Schettino, V.; Bini, R. Constraining Molecules at the Closest Approach: Chemistry at High Pressure. Chem. Soc. Rev. 2007, 36, 869−880. (d) Bini, R.; Ceppatelli, M.; Citroni, M.; Schettino, V. From Simple to Complex and Backwards. Chemical Reactions Under Very High Pressure. Chem. Phys. 2012, 398, 262−268. (5) (a) Classtone, S.; Laidler, K. J.; Eyring, H. The Theory of Rate Processes: the Kinetics of Chemical Reactions, Viscosity, Diffusion And Electrochemical Phenomena; McGraw-Hill: New York, 1941. (b) Hehre, W. J.; Radom, L.; Schleyer, P. V.; Pople, J. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (c) Shaik, S. S.; Schlegel, H. B.; Wolfe, S. Theoretical Aspects of Physical Organic Chemistry: The SN2 Mechanism; Wiley: New York, 1992. (6) (a) le Noble, W. J.; Kelm, H. Chemistry in Compressed Solutions. Angew. Chem., Int. Ed. 1980, 19, 841−856. (b) van Eldik, R.; Kelm, H. Interpretation of the Volume of Activation of Inorganic Reactions in Solution. Rev. Phys. Chem. Jpn. 1980, 50, 185−206. (c) van Eldik, R.; Hubbard, C. D. Elucidation of Inorganic Reaction Mechanisms. Application of High Pressure. Afr. J. Chem. 2000, 53, 139−168. (d) van Eldik, R.; Hubbard, C. D. The Interpretation and Mechanistic Significance of Activation Volumes for Organometallic Reactions. Adv. Phys. Org. Chem. 2006, 41, 1−78. (e) Hubbard, C. D.; van Eldik, R. Mechanistic Information on Some Inorganic and Bioinorganic Reactions from Volume Profile Analysis. Inorg. Chim. Acta 2010, 363, 2357−2374. (7) (a) le Noble, W. J. The Volume Profile of Chemical Reactions. Rev. Phys. Chem. Jpn. 1980, 50, 207−216. (b) Basilevsky, M. V.; Weinberg, N. N.; Zhulin, V. M. Pressure-Dependence of Activation and Reaction Volumes. J. Chem. Soc., Faraday Trans. 1 1985, 81, 875− 884. 776



Article

(21) Leach, A. R. Molecular Modelling: Principles and Applications, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 2001. (22) Burden, R. L.; Faires, J. D. Numerical Analysis, 9th ed.; Brooks/ Cole: Boston, MA, 2011. (23) (a) Marcus, R. Theoretical Relations Among Rate Constants Barriers and Bronsted Slopes of Chemical Reactions. J. Phys. Chem. 1968, 72, 891−899. (b) Marcus, R. Energetic and Dynamical Aspects of Proton Transfer Reactions in Solution. Faraday Symp. Chem. Soc. 1975, 10, 60−68. (24) (a) Hooper, J.; Zurek, E. Rubidium Polyhydrides Under Pressure: Emergence of the Linear H3− Species. J. Phys. Chem. C 2012, 116, 13322. (b) Hooper, J.; Zurek, E. Chem.Eur. J. 2012, 18, 5013− 5021. (c) Shamp, A.; Hooper, J.; Zurek, E. Compressed Cesium Polyhydrides: Cs+ Sublattices and H3− Three-Connected Nets. Inorg. Chem. 2012, 51, 9333−9342. (d) Hooper, J.; Altintas, B.; Shamp, A.; Zurek, E. Polyhydrides of the Alkaline Earth Metals: A Look at the Extremes under Pressure. J. Phys. Chem. C 2013, 117, 2982−2992. (25) (a) Ayouz, M.; Dulieu, O.; Guérout, R.; Robert, J.; Kokoouline, V. Potential Energy and Dipole Moment Surfaces of H3− Molecule. J. Chem. Phys. 2010, 132, 194309. (b) Matta, C. F.; Huang, L.; Massa, L. Characterization of a Trihydrogen Bond on the Basis of the Topology of the Electron Density. J. Phys. Chem. A 2011, 115, 12451−12458. (26) (a) Drickamer, H. G. Pi Electron Systems at High Pressure. Science 1967, 156, 1183−1189. (b) Bastron, V. C.; Drickamer, H. G. Solid State Reactions in Organic Crystals at Very High Pressure. J. Solid State Chem. 1971, 3, 550−563. (c) Pruzan, P.; Chervin, J. C.; Thiery, M. M.; Itie, J. P.; Besson, J. M.; Forgerit, J. P.; Revault, M. Transformation of Benzene to a Polymer After Static Pressurization to 30 GPa. J. Chem. Phys. 1990, 92, 6910−6915. (d) Santoro, M.; Ciabini, L.; Bini, R.; Schettino, V. High-Pressure Polymerization of Phenylacetylene and of the Benzene and Acetylene Moieties. J. Raman Spectrosc. 2003, 34, 557−566.

777


Reaction Profiles and Energy Surfaces of Compressed Species - The

Recommend Documents