Position-Specific and Clumped Stable Isotope Studies: Comparison of

Article pubs.acs.org/JPCA

Position-Specific and Clumped Stable Isotope Studies: Comparison of the Urey and Path-Integral Approaches for Carbon Dioxide, Nitrous Oxide, Methane, and Propane Michael A. Webb and Thomas F. Miller, III* Department of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, California 91125, United States S Supporting Information *

ABSTRACT: We combine path-integral Monte Carlo methods with high-quality potential energy surfaces to compute equilibrium isotope effects in a variety of systems relevant to ‘clumped’ isotope analysis and isotope geochemistry, including CO2, N2O, methane, and propane. Through a systematic study of heavy-atom isotope-exchange reactions, we quantify and analyze errors that arise in the widely used Urey model for predicting equilibrium constants of isotope-exchange reactions using reduced partition function ratios. These results illustrate that the Urey model relies on a nontrivial cancellation of errors that can shift the apparent equilibrium temperature by as much as 35 K for a given distribution of isotopologues. The calculations reported here provide the same level of precision as the best existing analytical instrumentation, resolving the relative enrichment of certain isotopologues to as little as 0.01‰. These findings demonstrate path-integral methods to be a rigorous and viable alternative to more approximate methods for heavy-atom geochemical applications.

1. INTRODUCTION Stable isotope analyses are essential for the understanding of many atmospheric, environmental, biological, geochemical, and astrochemical processes.1−10 Recently developed methods for experimental analysis of isotopic enrichment can detect the enhanced thermodynamic stability of a specific placement or ‘clumping’ of rare isotopes.11−16 Position-specific or clumped isotope enrichment results from homogeneous isotopic fractionation, which is quantified by ⎤ ⎡ (xi /x0)eq Δi = 1000⎢ − 1⎥ ⎢⎣ (xi /x0)r ⎦⎥

Though the effects of homogeneous isotopic fractionation may be subtle, modern instruments can determine Δi to precisions of 0.01−0.02‰.13,18−20 Relating the enrichment of isotopologues to their thermal populations (eq 1) enables the correlation of measurements of isotopic enrichment with temperatures in the geological record, thus facilitating applications that include the reconstruction of ancient marine ocean temperatures, the determination of body temperatures of extinct vertebrates, and the assessment of hydrocarbon deposits and organic matter.11,21−25 Clumped and position-specific isotope analyses are also used to characterize the production and consumption mechanisms of greenhouse gases.5,12,14,15,26,27 Theoretical methods are essential for both understanding and predicting stable isotope fractionation.17,24,28−33 In particular, theoretical predictions regarding isotopic clumping can be used to establish an absolute reference frame for standardizing isotope ratio measurements, which assist in interlaboratory reproducibility.18 For decades, the primary theoretical framework for

(1)

the degree to which a given isotopologue is enriched in the equilibrium ensemble due to its relative thermodynamic stability.17 Here, xi/x0 is the abundance of a particular isotopologue, i, relative to that with no rare isotope substitutions; (···)eq indicates quantities obtained from the equilibrium thermal distribution; and (···)r indicates quantities obtained from the stochastic distribution, in which the abundances of isotopologues are strictly determined by the abundances of their composite isotopes. © 2013 American Chemical Society

Received: November 12, 2013 Revised: December 20, 2013 Published: December 27, 2013 467

dx.doi.org/10.1021/jp411134v | J. Phys. Chem. A 2014, 118, 467−474

The Journal of Physical Chemistry A

Article

substitution. The relative enrichment of (Xm−1X′Yn)i for a set of t possible isotopomers is given by the expression33

characterizing equilibrium isotope effects has been the Urey model.17,24,33−35 In this approach, however, effects associated with vibrational anharmonicity and ro-vibrational coupling are ignored, which may become problematic as analytical instruments resolve increasingly smaller temperature differences. Although anharmonic and other corrections to the Urey model have been developed,36−38 the correction terms are generally complicated and have not been widely employed. The Feynman path-integral (PI) formulation of quantum statistical mechanics39 provides a useful framework for calculating equilibrium isotope effects. By exploiting the mathematical isomorphism between the quantum Boltzmann statistics of a given system and the classical Boltzmann statistics of its ring-polymer representation,40−43 PI methods have been applied to study equilibrium isotope effects, particularly for hydrogen-containing molecules, in both gas-phase44−50 and condensed-phase51−64 systems. In this study, PI Monte Carlo (PIMC) is used to predict the temperature-dependence of equilibrium constants in isotopeexchange reactions featuring heavy-atom isotope exchange between isotopologues of CO2, N2O, methane, and propane. The reported calculations resolve the enrichment of isotopologues to the same level as experimental resolution, while employing accurate, isotopically independent potential energy surfaces for CO2, N2O, and methane. These simulations are used to identify and examine breakdowns in the assumptions of the Urey model. The impact of these potential errors on the determination of apparent equilibrium temperatures is further explored.

⎡ ⎤ mSi ⎥ Δ(X m − 1X′Yn)i = 1000⎢ − 1 ⎢ 1 + ∑t K ij ⎥ ⎣ ⎦ j≠i

where Si = σi/σ0 is the ratio of the symmetry number for the isotopologue (Xm−1X′Yn)i and for that with no rare isotopes; the symmetry number is defined as the number of configurations that are indistinguishable by rotation.65 This relation has a leadingorder error that is proportional to the relative abundances of minority isotopes,33 which are negligible for the specific reactions considered here. Equations 3 and 5 have a clear connection to the temperature of the thermal ensemble through the various equilibrium constants. Differences between PI and Urey model calculations of the equilibrium constant can thus be used to identify errors in the determination of apparent equilibrium temperatures using the Urey model. 2.1. Urey Model. The Urey model uses the rigid-rotor and harmonic-oscillator approximations to compute partition function ratios (PFRs) that determine the relative abundances of isotopologues. By treating the rotational motions classically, the Teller−Redlich product rule65,66 is applied to avoid explicitly calculating the molecular moments of inertia, such that the total PFR between two isotopologues is N ⎛ m′ ⎞ Q′ σ = e−β(E0′− E0) ∏ ⎜ i ⎟ Q m σ′ i=1 ⎝ i ⎠

2. METHODOLOGY We aim to investigate the relative enrichment of isotopologues that are of primary interest in clumped and position-specific isotope analysis. For clumped isotope analysis, we assess the balance of thermodynamic equilibrium in isotopic clumping reactions that principally influence the abundance of clumped isotopologues, such as K

XY + X′Y′ ⇌ XY′ + X′Y

α

∏ j=1

ωj′ 1 − exp[−β ℏωj] ωj 1 − exp[−β ℏωj′]

where σ and σ′ again indicate the rotational symmetry numbers, β = 1/(kBT) is the inverse temperature, E0 is the zero-point energy, mi is the mass of the ith atom in a molecule of N atoms, ωj is the harmonic frequency of the jth normal mode, and α is the total number of normal vibrational modes (α = 3N − 5 for linear molecules and 3N − 6 for nonlinear molecules). The utility of eq 6 is that it reduces the problem of calculating PFRs to that of determining the relative harmonic vibrational frequencies for different isotopologues. The zero-point energy is typically calculated from harmonic vibrational contributions, such that E0 = Σjα= 1ℏωj/2, although anharmonic contributions to the zero-point energy have also been approximately included in the context of heavy-atom isotope exchange.17,33 In the following, calculations that employ the Urey model with a harmonic treatment of the zero-point energy are denoted UreyHO; Urey model calculations that include anharmonic corrections to the zero-point energy are denoted Urey-AHO. We note that it is possible to obtain higher-order corrections to the Urey model36−38 or to perform explicit eigenstate summation67,68 to calculate statistical mechanical properties while including effects like vibrational anharmonicity and rovibrational coupling; however, these approaches typically require the availability of extensive spectroscopic data and are not considered here. 2.2. Path Integral Statistics. The exact quantum mechanical canonical partition function for a system that exhibits Boltzmann statistics is given by the trace of the equilibrium ̂ density operator, Q = Tr(e−βH).39 The primitive, discretized PI representation of the partition function for a system of N distinguishable particles takes the form of a classical configuration integral,

(2)

(3)

where Kr is the equilibrium constant in eq 2 for the stochastic distribution; this relation has a leading-order error of 6 (ΔX′Y + ΔXY′), which is generally small when there are no structural isotopomers for singly substituted species (i.e., isomers with the same number of each isotope), as is the case for the specific reactions considered here.17,33 For position-specific isotope analysis, we assess the balance of thermal equilibrium in isotopomerization reactions that influence the relative isotopic enrichment of specific positions in a molecule, such as K ij

(X m − 1X′Yn)i ⇌ (X m − 1X′Yn)j

3/2

(6)

where K = (QXY′QX′Y)/(QXYQX′Y′). Here, X′ and Y′ are rare isotopes of X and Y, respectively, and X′Y′ is the isotopically clumped species. The equilibrium constant for eq 2 can be connected to the experimentally measurable enrichment of X′Y′ as given by eq 1, ΔX′Y′, via the relationship17 ΔX ′ Y ′ = −1000ln(K /K r)

(5)

(4)

where Kij = (Q(Xm−1X′Yn)j)/(Q(Xm−1X′Yn)i). Here, (Xm−1X′Yn)i and (Xm−1X′Yn)j indicate isotopomers of XmYn with a single X′ 468



Article

Table 1. Number of Beads and MC steps for PIMC Simulations CO2a

N2O

T (K)

P

MC steps/10

300 400 500 600

720/46 480/30 360/22 290/20

5.0/25.0 5.0/25.0 5.0/25.0 5.0/25.0

9

propaneb

methane

P

MC steps/10

720 450 350 250

2.5 2.5 2.5 2.5

9

P

MC steps/10

128 96 80 70

5.0 5.0 5.0 5.0

9

P

MC steps/109

90/246 58/126 41/78 32/55

5.0/5.0 5.0/5.0 5.0/5.0 5.0/5.0

a

Numbers on the left are for the HB simulations; numbers on the right are for the LB simulations. bNumbers on the left are for the exchange in eq 12; numbers on the right are for the exchange in eq 13.

1 Q (N , β) = lim P →∞ σ e

⎛ miP ⎞3P /2 ∏⎜ 2⎟ i = 1 ⎝ 2πβ ℏ ⎠ N

N

P

∫ ∏∏

12

(7)

Here, P indicates the number of ring-polymer beads, r(k) j indicates the Cartesian position of the jth atom in the kth ring-polymer bead, and βP = β/P. The effective ring-polymer potential is N

UP({r (j k)}) =

P

⎤

⎡

∑ ∑ ⎣⎢ 1 mjωP 2(r(j k) − r(j k− 1))2 ⎦⎥ j=1 k=1

12

∑

...,

rN(k))

(8)

k=1

where ωP = 1/(βPℏ) is the intrabead vibrational frequency, r(0) = r(P), and U(r1,...,rN) is the Born−Oppenheimer PES for the system. To enable direct comparison with the Urey model, the indistinguishability of identical nuclei in the PI calculations is treated in eq 7 using the classical, rotational symmetry number, σ. We note that PI descriptions that explicitly account for nuclear exchange statistics have been previously developed;43,69 regardless, such effects are expected to be negligible for the temperatures (T ≥ 300 K) and molecular systems considered here.

12

K

O C O + 16O13C18O ⇌ 16O13C16O + 16O12C18O

(9)

K

N15 N16O ⇌ 15 N14 N16O

(10)

The equilibrium constant for this reaction is K = R N N16O/ R14N14→15N16O, where R14→15N14N16O = Q(15N14N16O)/Q(14N14N16O) and R14N14→15N16O = Q(14N15N16O)/Q(14N14N16O). In a third application, the enrichment of 13CH3D, a clumped isotopologue of methane with relevance to isotope studies of natural gas,24 is investigated using the isotopic-clumping reaction 14→15

(13)

The equilibrium constant for this reaction is K = (R1→2Ht)/ (R 1 → 2 H c ), where R 1 → 2 H t = Q( 1 2 CH 2 D 1 2 CH 2 1 2 CH 3 )/ Q(12CH312CH212CH3) and R1→2Hc = Q(12CH312CHD12CH3)/ Q(12CH312CH212CH3). 3.2. Potential Energy Surfaces. The simulations for CO2 and N2O use intramolecular potential energy surfaces from Zúñiga and co-workers70−72 for which the coefficients of a fourth-order Morse-cosine expansion are determined using spectroscopic data. The simulations for methane utilize a PES computed at the CCSD(T) level of theory.73,74 The simulations for propane use the empirical CHARMM PES with general force field parameters;75 although the propane PES is significantly lower in quality than those used for the other systems, it nonetheless enables us to identify deviations between PI calculations and those based on the Urey model. 3.3. PI Calculations. PIMC sampling trajectories are performed in Cartesian coordinates with an explicit staging transformation.76 The staging length, j, is set such that 35−45% of all proposed staging moves are accepted. Prior to any data collection, each sampling trajectory is equilibrated for 107 MC steps, with P/j staging moves (rounded up to the nearest integer) attempted per MC step. Thereafter, ring-polymer configurations are sampled every 25 MC steps to reduce correlation among successive samples from a given sampling trajectory. The number of MC steps run for each sampling trajectory is detailed in Table 1. The equilibrium constants and PFRs defined in eqs 9−12 are computed at T = 300, 400, 500, and 600 K using the weighted histogram analysis method (WHAM).77 Each PFR calculation requires several independent sampling trajectories. The ringpolymer potentials employed for these sampling trajectories are

The equilibrium constant for this reaction is K = R16O12→13C16O/ R16O12→13C18O, where R16O12→13C16O = Q(16O13C16O)/Q(16O12C16O) and R16O12→13C18O = Q(16O13C18O)/Q(16O12C18O). In a second application, the enrichment of 14N15N16O, which is a primary contributor to the overall enrichment of 15N at the central position of N2O and is thus a potentially useful tool for characterizing origin processes,4,27 is investigated using the isotopomerization reaction 14

K

CH312CHD12CH3 ⇌ 12CH 2D12CH 212CH3

3. CALCULATION DETAILS 3.1. Systems. We perform calculations on four molecules that are of current interest in stable isotope studies: CO2, N2O, methane, and propane. In a first application, the enrichment of 16 13 18 O C O, which dominates the mass-47 experimental signature in clumped isotope studies of CO2,5 is investigated using the isotopic-clumping reaction 16 12 16

(12)

The equilibrium constant for this reaction is K = (R12→13Ct)/ (R 1 2 → 1 3 C c ), where R 1 2 → 1 3 C t = Q( 1 3 CH 3 1 2 CH 2 1 2 CH 3 )/ Q(12CH312CH212CH3) and R12→13Cc = Q(12CH313CH212CH3)/ Q(12CH312CH212CH3). In a final application, the enrichment of 12CH312CHD12CH3, a deuterium-substituted isotopomer of propane, is investigated using the isotopomerization reaction

P

+

K

CH313CH 212CH3 ⇌ 13CH312CH 212CH3

2

U (r1(k) ,

(11)

The equilibrium constant for this reaction is K = (R12→13CH4)/ (R12→13CH3D), where R12→13CH4 = Q(13CH4)/Q(12CH4) and R12→13CH3D = Q(13CH3D)/Q(12CH3D). In a fourth application, the enrichment of 12CH313CH212CH3, 13 a C-substituted isotopomer of propane, is investigated using the isotopomerization reaction

j=1 k=1

−βP UP({r (j k)})

K

CH4 + 13CH3D ⇌ 13CH412CH3D

dr (j k)

14

469



Article

4. RESULTS 4.1. Carbon Dioxide, Nitrous Oxide, and Methane. We begin by considering the isotopic enrichment of 16O13C18O, 14 15 16 N N O, and 13CH3D. Figure 1 presents Δ16O13C18O, Δ14N15N16O,

identical except for the isotope masses that appear in eq 8. For each PFR calculation, the independent sampling trajectories employ atomic masses that correspond to the two ‘end point’ isotopologues, as well as fractional atomic masses that lie intermediate to these end point isotopologues; the end point isotopologues are those that define the PFR. For example, to compute R14→15N14N16O, the different sampling trajectories utilize masses in the ring-polymer potential that are evenly spaced from 14 N to 15N. Nine sampling trajectories are used to compute the PFRs detailed in eqs 9, 10, and 13; five sampling trajectories are used to compute the PFRs in eqs 11 and 12. Statistical uncertainties for the PIMC calculations are reported as the standard error of the mean obtained from block-averaging the sampling trajectories.78 The convergence of the WHAM calculations is tested with multiple sensitivity checks. In particular, the PFR calculations are repeated after removing/ adding entire sampling trajectories at some of the intermediate isotope masses to check that the PFRs are converged with respect to the overlap of neighboring probability distributions. Additionally, the PFR calculations are repeated with different histogram resolutions to confirm that the bin sizes do not bias the results. To ensure that the number of ring-polymer beads is sufficiently large to converge the PI calculations, tests are performed in which the analytical expression for the primitive PI discretization of the partition function for a simple harmonic oscillator79 is used to compute PFRs and equilibrium constants. The harmonic oscillator frequencies that are employed in these convergence tests correspond to the normal mode vibrational frequencies for each molecule, and the number of beads, P, that are required to converge the PFRs or equilibrium constants to within a specified tolerance of the exact result is determined. In general, PI calculations of equilibrium constants converge more rapidly than PFRs as a function of P, due to a systematic cancellation of errors in the equilibrium constant calculation. These tests indicate that equilibrium constants can be computed using fewer ring-polymer beads in the following applications, resulting in less statistical variance and computational cost. The bead number for each PIMC sampling trajectory is detailed in Table 1. Two sets of sampling trajectories are employed for calculating the equilibrium constant for the isotopic-clumping reaction for CO2 (eq 9). The first, denoted HB, uses values of P needed to converge the PFRs to within 5 × 10−6; the second set, denoted LB, uses values of P needed to converge the equilibrium constant to within 10−5. PFR convergence tolerances of 5 × 10−6, 6 × 10−5, 10−4, and 4 × 10−3 are used for the PIMC calculations pertaining to eqs 10, 11, 12, and 13, respectively. These parameters are based on the limiting experimental precision of 0.01−0.2‰ in the measurement of Δi.13,15,19,24 The Urey-HO calculations employ vibrational frequencies obtained from the potential energy surfaces described in section 3.2. Anharmonic corrections to the zero-point energy are obtained from the literature17,33,72 and used to compute UreyAHO results for CO2, N2O, and methane. The anharmonic zeropoint energies reported for the CO2 and N2O isotopologues were obtained with the same potential energy surfaces as those described in section 3.2; the anharmonic zero-point energies reported for the methane isotopologues were obtained using B3LYP with the 6-311+G(3df,2p) basis set in contrast to the PES reported in section 3.2. The harmonic frequencies for the UreyHO calculations and the anharmonic zero-point energies for the Urey-AHO calculations are provided in the Supporting Information.

Figure 1. Δ16O13C18O, Δ14N15N16O, and Δ13CH3D as functions of 1000/T. PIMC results are indicated by circles; the error bars are smaller than the symbol size. PIMC results for Δ16O13C18O correspond to the set of calculations that employ a lower number of ring-polymer beads (LB). The solid and dashed lines correspond to the Urey-HO and Urey-AHO results, respectively. Note that the scaling of the y-axis differs above and below Δ = 6‰.

and Δ13CH3D as functions of temperature, computed with both PIMC and Urey model calculations; detailed numerical results are reported in Table 2. With the exception of the HB simulations for CO2, the statistical errors of the PIMC calculations are better than or comparable to the resolving power of analytical instrumentation, demonstrating the capability to precisely determine Δi using PI methods. Figure 1 illustrates that the Urey-HO results agree well with PIMC calculations for some of the systems considered, but not for others. For Δ16O13C18O, the results for both the Urey-HO and Urey-AHO methods are in good agreement with the PIMC calculations, with Urey-AHO results deviating from the PIMC calculations by less than 0.01‰. In contrast, the Urey-HO and Urey-AHO calculations are clearly different for both Δ14N15N16O and Δ13CH3D. For Δ14N15N16O, only the Urey-AHO results agree with the PIMC calculations to within statistical error; the Urey-HO results are in error by 0.3−0.6‰ For Δ13CH3D, the reverse is true; only the Urey-HO results are within statistical error of the PIMC calculations. The uneven performance of the Urey model can be understood by examining errors in the Urey-HO calculation of the PFRs. For a given PFR, R, the relative error (per mil) in the Urey-HO calculation with respect to the corresponding PIMC calculation is ⎛ R(Urey‐HO) ⎞ MR = 1000⎜⎜ (PIMC) − 1⎟⎟ ⎝ R ⎠ (PIMC)

(14)

(Urey‑HO)

where R and R are the PFRs computed using PIMC and the Urey-HO method, respectively. Figure 2 reports this quantity for a variety of PFRs. Interestingly, it is evident from the 470



Article

Table 2. Δi Values for Carbon Dioxide, Nitrous Oxide, and Methanea Δ16O13C18O

a

Δ14N15N16O

Δ13CH3D

T (K)

PIMC-HB

PIMC-LB

Urey-HO

Urey-AHO

PIMC

Urey-HO

Urey-AHO

PIMC

Urey-HO

Urey-AHO

300 400 500 600

0.96(8) 0.58(5) 0.36(4) 0.25(6)

0.960(9) 0.571(9) 0.357(7) 0.232(5)

0.938 0.556 0.350 0.228

0.949 0.564 0.357 0.234

22.85(4) 15.33(3) 11.06(3) 8.44(3)

23.54 15.85 11.53 8.79

22.81 15.31 11.10 8.43

5.73(5) 3.50(4) 2.29(4) 1.51(4)

5.71 3.46 2.23 1.50

5.33 3.18 2.01 1.31

Statistical errors in the final digit of the PIMC calculations are indicated in parentheses. All Δi values are reported in units of ‰.

associated with Urey-HO calculations are negligible for Δ13CH3D and small for Δ16O13C8O. In contrast, temperature errors as large as 15 K are observed for the Urey-HO calculation of Δ14N15N16O. For the Urey-AHO calculations, the predicted values of Δ16O13C18O and Δ14N15N16O do not result in statistically significant temperature errors; however, the temperature errors for Δ13CH3D are found to be as high as 35 K. We emphasize that the errors for the Urey-HO calculations of Δ14N15N16O and the Urey-AHO calculations of Δ13CH3D are both clearly larger than the statistical error of the PIMC calculations and the associated resolution of experimental measurements.13,15,24 Taken together, Figures 1−3 illustrate potential pitfalls of inherent approximations in Urey model calculations for predicting Δi values even in simple molecules, irrespective of whether the anharmonic correction to the zeropoint energy is applied. Furthermore, these results demonstrate that the use of approximate theoretical methods for computing relative isotopic enrichment can lead to mis-calibration of experimental isotope ratio measurements,18 leading to errors of as large as 35 K in the determination of apparent equilibrium temperatures. 4.2. Propane. To test the assumptions of the Urey model in a molecule that exhibits torsional motions, we investigate the isotopic enrichment of the propane isotopologues 12 CH312CHD12CH3 and 12CH313CH212CH3. Figure 4 presents PIMC and Urey-HO calculations for Δ12CH313CH212CH3 and Δ12CH312CHD12CH3 as functions of temperature; detailed numerical results are reported in Table 3. For the 12CH2D12CH212CH3

Figure 2. Relative error (per mil) of Urey-HO calculations of PFRs for the isotope-exchange reactions given by eqs 9−11. The error bars are smaller than the symbol size.

figure that the predictions for the Urey-HO method have errors that range from 1−5‰. These PFR errors are very large in comparison to the experimental resolution on Δi and the statistical uncertainty of the PIMC calculations. For CO2 and methane, these relatively large PFR errors precisely cancel when calculating the equilibrium constant, giving rise to the agreement between the Urey-HO and PIMC calculations of Δi in Figure 1. In contrast, the errors for R14→15N14N16O and R14N14→15N16O differ in magnitude such that the errors only partially cancel when computing the equilibrium constant in eq 10, resulting in the residual difference between the Urey-HO and PIMC calculations of Δ14N15N16O. Inaccurate estimates of Δi values based on Urey model calculations result in errors in determination of the apparent equilibrium temperature. Figure 3 quantifies these temperature errors by comparing the apparent equilibrium temperature for a given Δi, as determined from Urey model calculations versus the PIMC results. As expected from Figure 1, the temperature errors

Figure 4. Δ12CH312CHD12CH3 and Δ12CH313CH212CH3 as functions of 1000/T, with PIMC results given by circles and Urey-HO results given by the solid lines. Note that the scaling of the y-axis differs above and below Δ = 11‰.

Figure 3. Errors in the apparent equilibrium temperature obtained for Δ16O13C18O, Δ14N15N16O, and Δ13CH3D using (a) Urey-HO calculations and (b) Urey-AHO calculations. 471



Article

Table 3. Δi Values for Propanea Δ12CH312CHD12CH3

temperature errors for the Urey-HO calculations of both Δ12CH313CH212CH3 and Δ12CH312CHD12CH3 are observed to be near or less than 5 K. This is partially due to the fact that large changes in Δ12CH313CH212CH3 and Δ12CH312CHD12CH3 correspond to small changes in the apparent equilibrium temperature when the temperature is low. Consequently, the errors in Δ 1 2 CH 3 1 3 CH 2 1 2 CH 3 and Δ12CH312CHD12CH3 at low temperatures do not drastically change the apparent equilibrium temperature.

Δ12CH313CH212CH3

T (K)

PIMC

Urey

PIMC

Urey

300 400 500 600

72.1(3) 40.5(2) 25.4(1) 16.7(1)

69.1 39.6 24.9 16.7

9.74(2) 5.42(1) 3.09(1) 1.90(1)

9.58 5.29 3.07 1.84

a Statistical errors in the final digit of the PIMC calculations are indicated in parentheses. All Δi values are reported in units of ‰.

5. CONCLUSIONS In this work, we utilize path-integral methods to accurately characterize the equilibrium enrichment of certain isotopologues, Δi, in four molecules: CO2, N2O, methane, and propane. It is shown that PIMC methods combined with high-quality potential energy surfaces enable the determination of Δi to the same level of precision as the best analytical instrumentation currently available. These capabilities are used to demonstrate that the Urey-HO model relies upon a substantial cancellation of errors in partition function ratios to calculate Δi. Errors in Δi are observed when partition function ratio errors do not precisely cancel, as shown for Δ14N15N16O. We additionally find that applying anharmonic corrections to the zero-point energy does not reliably improve results of the Urey model. Using Urey model predictions for Δi, with or without anharmonic corrections, is found to lead to experimentally resolvable errors in the apparent equilibrium temperature of up to 35 K for the isotopologues studied here. The use of path-integral methods neither relies on uncontrolled cancellation of errors nor requires any a priori assumptions about the relative importance of effects such as anharmonicity and ro-vibrational coupling. These results demonstrate that PIMC is an accurate and feasible method for clumped and position-specific isotope analyses, as well as other heavy-atom geochemical applications.

isotopologue, deuterating the terminal methyl group can result in trans, gauche-plus, and gauche-minus rotamers with respect to the carbon backbone. Although these rotamers have identical equilibrium geometries, the trans and gauche rotamers have different normal vibrational frequencies. Consequently, the reported Urey-HO results for Δ12CH312CHD12CH3 are obtained by Boltzmann-averaging R1→2Ht over these rotamers.80 In Figure 4, the Urey-HO results for both Δ12CH313CH212CH3 and Δ12CH312CHD12CH3 display very little deviation from PIMC for T ≥ 500 K; however, larger deviations are found at lower temperatures. Figure 5 illustrates that the good agreement in Figure 4 is

■

Figure 5. Relative error (per mil) of the Urey-HO calculations of PFRs for the isotope-exchange reactions given by eqs 12 and 13. The error bars are smaller than the symbol size.

ASSOCIATED CONTENT

S Supporting Information *

once again due to a cancellation of errors during the Urey-HO calculation of the equilibrium constant. For the calculation of Δ12CH312CHD12CH3, the errors in R1→2Ht and R1→2Hc are determined to be nearly 20‰. Nevertheless, the error in the Urey-HO calculation of Δ12CH312CHD12CH3 is determined to be only as large as 2‰, which is smaller than might be anticipated. For R12→13Ct and R12→13Cc, the associated errors in the Urey-HO model are less than 2‰, and they too mostly cancel when calculating the equilibrium constant for eq 12. Errors in the apparent equilibrium temperature when determined with the Urey model are shown in Figure 6. The

PIMC results for the partition function ratios and equilibrium constants used to compute Δi estimates are provided. Additionally, the harmonic frequencies and anharmonic zero-point energies necessary to compute the Urey-HO and Urey-AHO results for CO2, N2O, and methane are included. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*(T.F.M.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the Resnick Sustainability Institute and the Department of Energy (DE-SC0006598). We acknowledge computing resources from the National Energy Research Scientific Computing Center (DE-AC02-05CH11231) and a Director’s Discretionary Allocation from the Argonne Leadership Computing Facility. We thank José Zúñiga and Tim Lee for providing information and subroutines for their potential energy surfaces. We also thank John Eiler, Alex Sessions, and Adam Subhas for helpful discussions.

Figure 6. Errors in the apparent equilibrium temperature obtained for Δ12CH312CHD12CH3 and Δ12CH313CH212CH3 using Urey-HO calculations. 472



■

Article

(22) Finnegan, S.; Bergmann, K.; Eiler, J. M.; Jones, D. S.; Fike, D. A.; Eisenman, I.; Hughes, N. C.; Tripati, A. K.; Fischer, W. W. The Magnitude and Duration of Late Ordovician-Early Silurian Glaciation. Science 2011, 331, 903−906. (23) Eagle, R. A.; Tütken, T.; Martin, T. S.; Tripati, A. K.; Fricke, H. C.; Connely, M.; Cifelli, R. L.; Eiler, J. M. Dinosaur Body Temperatures Determined from Isotopic (13C-18O) Ordering in Fossil Biominerals. Science 2011, 333, 443−445. (24) Ma, Q.; Wu, S.; Tang, Y. Formation and Abundance of Doublysubstituted Methane Isotopologues (13CH3D) in Natural Gas Systems. Geochim. Cosmochim. Acta 2008, 72, 5446−5456. (25) Tsuji, K.; Teshima, H.; Sasada, H.; Yoshida, N. Spectroscopic Isotope Ratio Measurement of Doubly-Substituted Methane. Spectrochim. Acta, Part A 2012, 98, 43−46. (26) Eiler, J. M.; Schauble, E. 18O13C16O in Earth’s Atmosphere. Geochim. Cosmochim. Acta 2004, 68, 4767−4777. (27) Toyoda, S.; Yano, M.; Nishimura, S.; Akiyama, H.; Hayakawa, A.; Koba, K.; Sudo, S.; Yagi, K.; Makabe, A.; Tobari, Y.; et al. Characterization and Production and Consumption Processes of N2O Emitted from Temperate Agricultural Soils Determined via Isotopomer Ratio Analysis. Global Biogeochem. Cycles 2011, 25, GB2008. (28) Schauble, E. A.; Rossman, G. R.; Taylor, H. P. Theoretical Estimates of Equilibrium Chlorine-Isotope Fractionations. Geochim. Cosmochim. Acta 2003, 67, 3267−3281. (29) Anbar, A. D.; Jarzecki, A. A.; Spiro, T. G. Theoretical Investigation 3+ of Iron Isotope Fractionation between Fe(H2O)3+ 6 and Fe(H2O)6 : Implications for Iron Stable Isotope Geochemistry. Geochim. Cosmochim. Acta 2005, 69, 825−837. (30) Schauble, E. A.; Ghosh, P.; Eiler, J. M. Preferential Formation of 13 C−18O Bonds in Carbonate Minerals, Estimated Using FirstPrinciples Lattice Dynamics. Geochim. Cosmochim. Acta 2006, 70, 2510−2529. (31) Guo, W.; Mosenfelder, J. L.; Goddard, W. A., III; Eiler, J. M. Isotopic Fractionations Associated with Phosphoric Acid Digestion of Carbonate Minerals: Insights from First-Principles Theoretical Modeling and Clumped Isotope Measurements. Geochim. Cosmochim. Acta 2009, 73, 7203−7225. (32) Ni, Y.; Ma, Q.; Ellis, G. S.; Dai, J.; Katz, B.; Zhang, S.; Tang, Y. Fundamental Studies on Kinetic Isotope Effect (KIE) of Hydrogen Isotope Fractionation in Natural Gas Systems. Geochim. Cosmochim. Acta 2011, 75, 2696−2707. (33) Cao, X.; Liu, Y. Theoretical Estimation of the Equilibrium Distribution of Clumped Isotopes in Nature. Geochim. Cosmochim. Acta 2012, 77, 292−303. (34) Urey, H. C. The Thermodynamic Properties of Isotopic Substances. J. Chem. Soc. 1947, 562−581. (35) Bigeleisen, J.; Mayer, M. G. Calculation of Equilibrium Constants for Isotopic Exchange Reactions. J. Chem. Phys. 1947, 15, 261−267. (36) Richet, P.; Bottinga, Y.; Jayvoy, M. Review of Hydrogen, Carbon, Nitrogen, Oxygen, Sulfur, and Chlorine Stable Isotope Fractionation among Gaseous Molecules. Annu. Rev. Earth Planet. Sci. 1977, 5, 65− 110. (37) Barone, V. Vibrational Zero-Point Energies and Thermodynamic Functions beyond the Harmonic Approximation. J. Chem. Phys. 2004, 120, 3059−3065. (38) Liu, Q.; Tossell, J. A.; Liu, Y. On the Proper Use of the Bigeleisen−Mayer Equation and Corrections to It in the Calculation of Isotopic Fractionation Equilibrium Constants. Geochim. Cosmochim. Acta 2010, 74, 6965−6983. (39) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill Companies: New York, 1965. (40) Chandler, D.; Wolynes, P. G. Exploiting the Isomorphism between Quantum-Theory and Classical Statistical-Mechanics of Polyatomic Fluids. J. Chem. Phys. 1981, 74, 4078−4095. (41) Parninello, M.; Rahman, A. Study of an F Center in Molten KCl. J. Chem. Phys. 1984, 80, 860. (42) De Raedt, B.; Sprik, M.; Klein, M. L. Computer Simulation of Muonium in Water. J. Chem. Phys. 1984, 80, 5719.

REFERENCES

(1) Schoell, M. Recent Advances in Petroleum Isotope Geochemistry. Org. Geochem. 1984, 6, 645−663. (2) Whiticar, M. J. A Geochemical Perspective of Natural-Gas and Atmospheric Methane. Org. Geochem. 1990, 16, 531−547. (3) Conny, J. M.; Currie, L. A. The Isotopic Characterization of Methane, Non-Methane Hydrocarbons and Formaldehyde in the Troposphere. Atmos. Environ. 1996, 30, 621−638. (4) Stein, L. Y.; Yung, Y. L. Production, Isotopic Composition, and Atmospheric Fate of Biologically Produced Nitrous Oxide. Annu. Rev. Earth Planet. Sci. 2003, 31, 329−356. (5) Affek, H. P.; Eiler, J. M. Abundance of Mass 47 CO2 in Urban Air, Car Exhaust, and Human Breath. Geochim. Cosmochim. Acta 2006, 70, 1−12. (6) Eiler, J. M. Clumped-Isotope” Geochemistry: The Study of Naturally-Occurring, Multiply-substituted Isotopologues. Earth Planet. Sci. Lett. 2007, 262, 309−327. (7) Sturup, S.; Hansen, H. R.; Gammelgaard, B. Application of Enriched Stable Isotopes as Tracers in Biological Systems: A Critical Review. Anal. Bioanal. Chem. 2008, 390, 541−554. (8) Dlugokencky, E. J.; Nisbet, E. G.; Fisher, R.; Lowry, D. Global Atmospheric Methane: Budget, Changes and Dangers. Philos. Trans. R. Soc., A 2011, 369, 2058−2072. (9) Tilley, B.; Muehlenbachs, K. Isotope Reversals and Universal Stages and Trends of Gas Maturation in Sealed, Self-Contained Petroleum Systems. Chem. Geol. 2013, 339, 194−204. (10) Webster, C. R.; Mahaffy, P. R.; Flesch, G. J.; Niles, P. B.; Jones, J. H.; Leshin, L. A.; Atreya, S. K.; Stern, J. C.; Christensen, L. E.; Owen, T.; et al. Isotope Ratios of H, C, and O in CO2 and H2O of the Martian Atmosphere. Science 2013, 341, 260−263. (11) Oba, Y.; Naraoka, H. Site-Specific Carbon Isotope Analysis of Aromatic Carboxylic Acids by Elemental Analysis/Pyrolysis/Isotope Ratio Mass Spectrometry. Rapid Commun. Mass Spectrom. 2006, 20, 3649−3653. (12) Mohn, J.; Tuzson, B.; Manninen, A.; Yoshida, N.; Toyoda, S.; Brand, W. A.; Emmenegger, L. Site Selective Real-Time Measurements of Atmospheric N2O Isotopomers by Laser Spectroscopy. Atmos. Meas. Tech. 2012, 5, 1601−1609. (13) Eiler, J. M.; Clog, M.; Magyar, P.; Piasecki, A.; Sessions, A.; Stolper, D.; Deerberg, M.; Schlueter, H.-J.; Schwieters, J. A HighResolution Gas-Source Isotope Ratio Mass Spectrometer. Int. J. Mass Spectrom. 2013, 335, 45−56. (14) Well, R.; Eschenbach, W.; Flessa, H.; von der Heide, C.; Weymann, D. Are Dual Isotope and Isotopomer Ratios of N2O Useful Indicators for N2O Turnover During Denitrification in NitrateContaminated Aquifers? Geochim. Cosmochim. Acta 2012, 90, 265−282. (15) Köster, J. R.; Well, R.; Tuzson, B.; Bol, R.; Dittert, K.; Giesemann, A.; Emmenegger, L.; Manninen, A.; Cárdenas, L.; Mohn, J. Novel Laser Spectroscopic Technique for Continuous Analysis of N2O Isotopomers: Application and Intercomparison with Isotope Ratio Mass Spectrometry. Rapid Commun. Mass Spectrom. 2013, 27, 216−222. (16) Eiler, J. M.; Jeanloz, R. Annu. Rev. Earth Planet. Sci. 2013, 41, 411− 441. (17) Wang, Z. G.; Schauble, E. A.; Eiler, J. M. Equilibrium Thermodynamics of Multiply Substituted Isotopologues of Molecular Gases. Geochim. Cosmochim. Acta 2004, 68, 4779−4797. (18) Dennis, K. J.; Affek, H. P.; Passey, B. H.; Schrag, D. P.; Eiler, J. M. Defining an Absolute Reference Frame for ‘Clumped’ Isotope Studies of CO2. Geochim. Cosmochim. Acta 2011, 75, 7117−7131. (19) Yoshida, N.; Vasilev, M.; Ghosh, P.; Abe, O.; Yamada, K.; Morimoto, M. Precision and Long-Term Stability of Clumped-Isotope Analysis of CO2 Using a Small-Sector Isotope Ratio Mass Spectrometer. Rapid Commun. Mass Spectrom. 2013, 27, 207−215. (20) Grauel, A.-L.; Schmid, T. W.; Hu, B.; Bergami, C.; Capotondi, L.; Zhou, L.; Bernasconi, S. M. Calibration and Application of the ‘Clumped Isotope’ Thermometer to Foraminifera for High-Resolution Climate Reconstructions. Geochim. Cosmochim. Acta 2013, 108, 125−140. (21) Eiler, J. M. Paleoclimate Reconstruction Using Carbonate Clumped Isotope Thermometry. Quat. Sci. Rev. 2011, 30, 3575−3588. 473



Article

(43) Ceperley, D. M. Path-Integrals in the Theory of Condensed Helium. Rev. Mod. Phys. 1995, 67, 279−355. (44) Böhm, M. C.; Schulte, J.; Hernández, E.; Ramírez, R. Electrons and Nuclei of Ethylene Isomers; a Feynman Path Integral-Ab Initio Study. Chem. Phys. 2001, 264, 371−400. (45) Miller, T. F., III; Clary, D. C. Torsional Path Integral Monte Carlo Method for the Quantum Simulation of Large Molecules. J. Chem. Phys. 2002, 116, 8262−8269. (46) Miller, T. F., III; Clary, D. C. Quantum Free Energies of the Conformers of Glycine on an Ab Initio Potential Energy Surface. Phys. Chem. Chem. Phys. 2004, 6, 2563−2571. (47) Lynch, V. A.; Mielke, S. L.; Truhlar, D. G. Accurate VibrationalRotational Partition Functions and Standard-State Free Energy Values for H2O2 from Monte Carlo Path-Integral Calculations. J. Chem. Phys. 2004, 121, 5148−5162. (48) Zimmermann, T.; Vaníček, J. Path Integral Evaluation of Equilibrium Isotope Effects. J. Chem. Phys. 2009, 131, 024111. (49) Pérez, A.; von Lilienfeld, O. A. Path Integral Computation of Quantum Free Energy Differences Due to Alchemical Transformations Involving Mass and Potential. J. Chem. Theory Comput. 2011, 7, 2358− 2369. (50) Mielke, S. L.; Truhlar, D. G. Accelerating the Convergence and Reducing the Variance of Path Integral Calculations of Quantum Mechanical Free Energies by Using Local Reference Potentials. J. Chem. Theory Comput. 2012, 8, 1589−1596. (51) Balog, E.; Hughes, A. L.; Martyna, G. J. Constant Pressure Path Integral Molecular Dynamics Studies of Quantum Effects in the Liquid State Properties of N-Alkanes. J. Chem. Phys. 2000, 112, 870−880. (52) Chen, B.; Ivanov, I.; Klein, M. L.; Parrinello, M. Hydrogen Bonding in Water. Phys. Rev. Lett. 2003, 91, 215503. (53) Hernández de la Peña, L.; Kusalik, P. Quantum Effects in Light and Heavy Liquid Water: A Rigid-Body Centroid Molecular Dynamics Study. J. Chem. Phys. 2004, 121, 5992−6002. (54) Lynch, V. A.; Mielke, S. L.; Truhlar, D. G. High-Precision Quantum Thermochemistry on Nonquasiharmonic Potentials: Converged Path-Integral Free Energies and a Systematically Convergent Family of Generalized Pitzer-Gwinn Approximations. J. Phys. Chem. A 2005, 109, 10092−10099. (55) Miller, T. F., III; Clary, D. C. Quantum Simulation of a Hydrated Noradrenaline Analog with the Torsional Path Integral Method. J. Phys. Chem. A 2006, 110, 731−740. (56) Paesani, F.; Iuchi, S.; Voth, G. A. Quantum Effects in Liquid Water from an Ab Initio-Based Polarizable Force Field. J. Chem. Phys. 2007, 127, 074506. (57) Mielke, S. L.; Truhlar, D. G. Improved Methods for Feymnan Path Integral Calculations of Vibrational-Rotational Free Energies and Application to Isotopic Fractionation of Hydrated Chloride Ions. J. Phys. Chem. A 2009, 113, 4817−4827. (58) Azuri, A.; Engel, H.; Doron, D.; Major, D. T. Path-Integral Calculations of Nuclear Quantum Effects in Model Systems, Small Molecules, and Enzymes via Gradient-Based Forward Corrector Algorithms. J. Chem. Theory Comput. 2011, 7, 1273−1286. (59) Herrero, C.; Ramírez, R. Isotope Effects in Ice Ih: A Path-Integral Simulation. J. Chem. Phys. 2011, 134, 094510. (60) Markland, T. E.; Berne, B. J. Unraveling Quantum Mechanical Effects in Water Using Isotopic Fractionation. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 7988−7991. (61) Liu, J.; Andino, R. S.; Miller, C. M.; Chen, X.; Wilkins, D. M.; Ceriotti, M.; Manolopoulos, D. E. A Surface-Specific Isotope Effect in Mixtures of Light and Heavy Water. J. Phys. Chem. C 2013, 117, 2944− 2951. (62) Ceriotti, M.; Markland, T. E. Efficient Methods and Practical Guidelines for Simulating Isotope Effects. J. Chem. Phys. 2013, 138, 014112. (63) Mielke, S. L.; Dinpajooh, M.; Siepmann, J. I.; Truhlar, D. G. Efficient Methods for Including Quantum Effects in Monte Carlo Calculations of Large Systems: Extension of the Displaced Points Path Integral Method and Other Effective Potential Methods to Calculate Properties and Distributions. J. Chem. Phys. 2013, 138, 014110.

(64) Buchowiecki, M.; Vaníček, J. Monte Carlo Evaluation of the Equilibrium Isotope Effects Using the Takahashi−Imada Factorization of the Feynman Path Integral. Chem. Phys. Lett. 2013, 588, 11−16. (65) McQuarrie, D. A. Statistical Mechanics; University Science Books: South Orange, NJ, 2000. (66) Redlich, O. A General Relationship between the Oscillation Frequency of Isotropic Molecules (with Remarks on the Calculation of Harmonious Force Constants). Z. Phys. Chem. 1935, 28, 371−382. (67) Stein, E.; Rabinovitch, B. S. Accurate Evaluation of Internal Energy Level Sums and Densities including Anharmonic Oscillators and Hindered Rotors*. J. Chem. Phys. 1973, 58, 2438−2445. (68) Beyer, T.; Swinehart, D. F. Algorithm 448: Number of MultiplyRestricted Partitions. Commun. ACM 1973, 16, 379. (69) Marx, D.; Müser, M. H. Path Integral Simulations of Rotors: Theory and Applications. J. Phys.: Condens. Matter 1999, 11, R117− R155. (70) Zúñiga, J.; Alacid, M.; Bastida, A.; Carvajal, F. J.; Requena, A. Determination of a Potential Energy Surface for CO2 Using Generalized Internal Vibrational Coordinates. J. Mol. Spectrosc. 1999, 195, 137−146. (71) Zúñiga, J.; Alacid, M.; Bastida, A.; Carvajal, F. J.; Requena, A. Determination of Highly Excited Rovibrational States for N2O Using Generalized Internal Coordinates. J. Chem. Phys. 1999, 110, 6339−6352. (72) Zúñiga, J.; Bastida, A.; Requena, A. Theoretical Calculations of Vibrational Frequencies and Rotational Constants of the N2O Isotopomers. J. Mol. Spectrosc. 2003, 217, 43−58. (73) Lee, T. J.; Martin, J. M. L.; Taylor, P. R. An Accurate Ab-Initio Quartic Force-Field and Vibrational Frequencies for CH4 and Isotopomers. J. Chem. Phys. 1995, 102, 254−261. (74) Dateo, C. E.; Lee, T. J.; Schwenke, D. W. An Accurate Quartic Force-Field and Vibrational Frequencies for HNO and DNO. J. Chem. Phys. 1994, 101, 5853−5859. (75) Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. J. Comput. Chem. 2010, 31, 671−690. (76) Tuckerman, M. E.; Berne, B. J.; Martyna, G. J.; Klein, M. L. Efficient Molecular-Dynamics and Hybrid Monte-Carlo Algorithms for Path-Integrals. J. Chem. Phys. 1993, 99, 2796−2808. (77) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. 1. The Method. J. Comput. Chem. 1992, 13, 1011−1021. (78) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: New York, 2002. (79) Schweizer, K. S.; Stratt, R. M.; Chandler, D.; Wolynes, P. G. Convenient and Accurate Discretized Path Integral Methods for Equilibrium Quantum Mechanical Calculations. J. Chem. Phys. 1981, 75, 1347−1364. (80) Wang, Y.; Sessions, A. L.; Nielsen, R. J.; Goddard, W. A., III. Equilibrium 2H/1H Fractionations in Organic Molecules: I. Experimental Calibration of Ab Initio Calculations. Geochim. Cosmochim. Acta 2009, 73, 7060−7075.

474


Position-Specific and Clumped Stable Isotope Studies: Comparison of

Recommend Documents