Importance of a Fully Anharmonic Treatment of Equilibrium Isotope

Article pubs.acs.org/accounts

Importance of a Fully Anharmonic Treatment of Equilibrium Isotope Fractionation Properties of Dissolved Ionic Species As Evidenced by Li+(aq) Romain Dupuis,† Magali Benoit,‡ Mark E. Tuckerman,¶,§,∥ and Merlin Méheut*,⊥ †

DIPC, Paseo Manuel de Lardizabal, 4, 20018 San Sebastiàn, Spain CEMES CNRS UPR8011, 29 rue Jeanne Marvig, 31055 Cedex Toulouse, France ¶ Department of Chemistry, New York University, New York, New York 10003, United States § Courant Institute of Mathematical Sciences, New York University, New York, New York 10012, United States ∥ NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China ⊥ GET, CNRS UMR 5563, IRD UR 154, Université Toulouse III, OMP, 14 avenue Edouard Belin, 31400 Toulouse, France ‡

S Supporting Information *

CONSPECTUS: Equilibrium fractionation of stable isotopes is critically important in fields ranging from chemistry, including medicinal chemistry, electrochemistry, geochemistry, and nuclear chemistry, to environmental science. The dearth of reliable estimates of equilibrium fractionation factors, from experiment or from natural observations, has created a need for accurate computational approaches. Because isotope fractionation is a purely quantum mechanical phenomenon, exact calculation of fractionation factors is nontrivial. Consequently, a severe approximation is often made, in which it is assumed that the system can be decomposed into a set of independent harmonic oscillators. Reliance on this often crude approximation is one of the primary reasons that theoretical prediction of isotope fractionation has lagged behind experiment. A class of problems for which one might expect the harmonic approximation to perform most poorly is the isotopic fractionation between solid and solution phases. In order to illustrate the errors associated with the harmonic approximation, we have considered the fractionation of Li isotopes between aqueous solution and phyllosilicate minerals, where we find that the harmonic approximation overestimates isotope fractionation factors by as much as 30% at 25 °C. Lithium is a particularly interesting species to examine, as natural lithium isotope signatures provide information about hydrothermal processes, carbon cycle, and regulation of the Earth’s climate by continental alteration. Further, separation of lithium isotopes is of growing interest in the nuclear industry due to a need for pure 6 Li and 7Li isotopes. Moving beyond the harmonic approximation entails performing exact quantum calculations, which can be achieved using the Feynman path integral formulation of quantum statistical mechanics. In the path integral approach, a system of quantum particles is represented as a set of classical-like ring-polymer chains, whose interparticle interactions are determined by the rules of quantum mechanics. Because a classical isomorphism exists between the true quantum system and the system of ring-polymers, classical-like methods can be applied. Recent developments of efficient path integral approaches for the exact calculation of isotope fractionation now allow the case of the aforementioned dissolved Li fractionation properties to be studied in detail. Applying this technique, we find that the calculations yield results that are in good agreement with both experimental data and natural observations. Importantly, path integral methods, being fully atomistic, allow us to identify the origins of anharmonic effects and to make reliable predictions at temperatures that are experimentally inaccessible yet are, nevertheless, relevant for natural phenomena.

1. INTRODUCTION

are present and often play a crucial role in important fractionation processes. This research cross-pollinates other areas of science such as electrochemistry3 and medicine.4 Broadly, the measured isotopic variations can be used either to identify the source of the measured material (tracing) or to assess

Since the discovery of isotopes in the early years of the last century, geochemists have invested a considerable effort in documenting variations in the isotopic signatures of natural phases.1,2 Currently, due to significant progress in mass spectrometry, isotope geochemistry has become a forerunner in the exploration of numerous new isotopic systems and their applications. In many of these new isotopic systems (e.g., Fe, Li, Cu, Zn, Ca), dissolved ionic species © 2017 American Chemical Society

Received: December 4, 2016 Published: June 23, 2017 1597

DOI: 10.1021/acs.accounts.6b00607 Acc. Chem. Res. 2017, 50, 1597−1605

Article

Accounts of Chemical Research

mechanics20,21 in the calculation of Li fractionation properties between solid and solution. The results of the two methods are compared for the fractionation properties of Li2O, of Li+ in water, and the value of Δ7LiA−B between aqueous solution and the phyllosilicates. We found that the HA does not hold in the case of dissolved Li+ while the PI treatment of the aqueous solution phase improves the agreement with experiment significantly. The PI simulations additionally yielded important insights into the origin of the anharmonicities that cause the discrepancy. We elucidate the atomistic origin of the anharmonicity via an analysis of relevant spatial distribution functions from the PI configurations. Finally, we speculate on the consequences of these findings for other isotopic systems of geochemical interest.

the mechanism or conditions of formation of the material. In all cases, a basic understanding of the fractionation process is necessary. In particular, the isotopic fractionation occurring in thermodynamic equilibrium is a key aspect. Theoretical calculations of equilibrium fractionation effects have accompanied isotopic measurements since their beginning. In most systems, equilibrium isotope fractionation results solely from a change in the thermovibrational energy associated with the substitution of one isotope for another and is a purely quantum mechanical phenomenon. The theoretical framework for computing this thermovibrational effect is traditionally based on the widely used harmonic approximation (HA).5 In this scheme, it is assumed that atomic motions can be approximated as a set of coupled harmonic oscillators. The thermovibrational energies are then obtained as the set of normal-mode frequencies. The HA has been widely employed for computing fractionation properties of solids and gases, and the calculated gas/gas, gas/solid, and solid/solid fractionation factors have shown reasonable agreement with experiment. These computations originally employed experimentally measured vibrational frequencies. While such measurements are relatively easy in the solid and gas phases, obtaining vibrational frequencies of liquids is considerably more challenging. Consequently, the approach was not applied to liquids until approximately two decades ago, when significant developments in electronic structure methods allowed the prediction of vibrational properties from first principles. In a first attempt to use these techniques to assess fractionation properties of dissolved species, some groups developed simplified models of liquids such as small cation/water molecular clusters6,7 and hydrated solids.8 More recently, some groups carried out calculations based on a more rigorous treatment of liquids by molecular dynamics. 9−11 However, in all these studies, fractionation properties were computed in the HA, the validity of which is theoretically questionable in a liquid. Unfortunately, validation of such calculations is rendered difficult by the lack of reliable experimental estimates of equilibrium fractionation. In one recent study of Pinilla et al.,12 it was shown that the fractionation of Mg2+ between aqueous solution and carbonate materials could be predicted accurately using the HA. Based on this result, these authors suggested that the HA might generally extend to other ions in solution. Recent developments in liquid simulation methods allow for a proper fully anharmonic treatment of thermovibrational properties of liquids, providing detailed atomistic insights into the fractionation process that can help understand experimental and natural observations and permit a rigorous test of the validity of the HA to fractionation properties of dissolved species.12−14 In light of these important advances, we sought to carry out a study of the equilibrium fractionation of Li isotopes between phyllosilicates (here, smectite and mica) and Li+ in solution. The fractionation of lithium is particularly interesting, as Li is a light, rapidly diffusing atom, suggesting that isotopic equilibrium is more easily reachable experimentally on observable time scales, allowing for a reliable comparison between theory and experiment. Moreover, because of its considerable geochemical interest for tracking fluid−mineral interactions,15−17 this equilibrium fractionation has benefited from high-quality experimental estimates over an extended temperature domain.18,19 As an important part of our study and in light of recent claims,12 we sought to compare the HA to an explicit treatment of quantum and thermal fluctuations of all nuclei via the Feynman path integral (PI) formulation of quantum statistical

2. THEORETICAL APPROACH Equilibrium fractionation properties of condensed phases are expressed theoretically by their β-factor, which is related to the change in free energy associated with isotopic substitution in a particular phase. By contrast, isotopic substitution measurements yield the variation in δ-values (δ7Li), which expresses the per mill (‰) variations in the isotopic signature of a given phase with respect to a recognized reference material (LSVEC for Li). (Natural Li isotopic signatures vary typically between −15‰ and +45‰.22) The connection between the δ-values and the β-factors is expressed via the isotopic fractionation between phases A and B: δ 7 Li − δ 7 Li ≡ Δδ 7 LiA − B = 1000 ln α 7 LiA − B AB measured by analysts

= 1000 ln β 7 LiA − 1000 ln β 7 LiB  calculated by theorists

(1)

where Δδ7LiA−B and 1000 ln α7LiA−B are referred to as the “fractionation between phases A and B” by analysts and theorists, respectively. The ln β7LiA(B) are the so-called β-factors of Li in phase A (B). The Li β-factor is related to a free energy difference associated with a change of the Li mass as ln β 7 Li =

m ΔF 3 − ln 7 kBT 2 m6

(2)

where ΔF = F7 − F6 is the free energy difference and m7(6) is the mass of the heavy (light) isotope. In this work, β-factors were calculated both for a solid and for a dissolved species. For these two cases, we compare the standard HA with the aforementioned PI approach, which we now briefly describe. The β-factor of a species at temperature T is, in this study, computed via a path-integral based molecular dynamics (PIMD) scheme combined with thermodynamic integration (TI). The path integral formalism specifies that each atom in an N-atom system having coordinate operators r̂1, ..., r̂N, masses m1, ..., mN, and interaction potential U(r̂1, ..., r̂N) be replaced by a ring-polymer containing P points or “beads” connected to their two nearest neighbors by harmonic springs of frequency ωP = √P/(βℏ), where β = 1/kBT. While exact quantum statistics are recovered in the limit that P → ∞, in practice, a finite value of P must be employed and convergence with respect to P should be checked. Under the assumption of Boltzmann statistics, the expression for the quantum canonical partition function, QP(N, V, T), at volume V and temperature T 1598


Article

Accounts of Chemical Research for finite P is 3P /2 ⎡ N ⎛ mP ⎞ Q P(N , V , T ) = ⎢∏ ⎜ i 2 ⎟ ⎢⎣ i = 1 ⎝ 2πβ ℏ ⎠

For solid Li2O, PIMD was used to compute the β-factor employing the same protocol as for Li+(aq), while the harmonic β-factor was obtained directly from the lattice dynamics. For the complex polylithionite structures, β-factors are computed only within the HA, which is assumed to be accurate for solids. In order to facilitate efficient TI-PIMD calculations, empirical potentials were employed for Li2O and Li+(aq). For phyllosilicate materials, accurate empirical potentials do not exist; hence we employed a density functional theory (DFT)-based approach. When combining our mica-type and Li+(aq) calculations, the consistency of employing different models of atomic interactions was thoroughly checked. In particular, in section 4.2 of SI, we show that in the HA, the DFT-based and empirical approaches give similar results for Li+(aq) fractionation properties.

⎤

(P)⎥ ∫ dr(1) i ... dri ⎥

⎦

P ⎡ N ⎧ ⎪ ⎢∑ 1 miωP 2(r(i k) − r(i k + 1))2 × exp⎨ − β ∑ ⎪ ⎣ i=1 2 ⎩ k=1 ⎢

+

⎤⎫ ⎪ 1 U (r1(k), ..., r(Nk))⎥⎬ ⎪ ⎥⎦⎭ P

(3)

where r(k) denote the coordinates of the kth bead of the ith i atom, with r(P+1) = r(1) i i . Thermodynamic averages are computed from a canonical distribution of these ring-polymers generated using MD. Finally, the free energy ΔF needed to determine the β-factor is obtained by interpolating the mass of Li from 6 amu to 7 amu by an external switching parameter λ, λ ∈ [0,1] and an interpola-

3. RESULTS 3.1. Anharmonic Effects in Li+(aq) vs Li2O

tion function 1/ m(λ) = (1 − λ)/ m6 + λ / m7 .13,23 The motivation for this inverse-square switching scheme is based on the fact that at low temperature, the leading term in F for a harmonic oscillator of mass m and force constant k is (ℏ/2) k /m , which is linear in 1/√m. The mass interpolation causes the partition function and free energy to acquire a λ dependence, FP(N, V, T, λ) = −kBT ln QP(N, V, T, λ), and ΔFP is then computed by thermodynamic integration ΔFP =

∫0

1

∂FP dλ = − ∂λ

Fractionation properties are generally approximately linear in 1/T or 1/T2,30 and recognizing this dependence, in Figure 1, the

1 N

dm ∫0 ∑ m 1(λ) dλi ⟨Tî⟩P ,λ dλ i=1

i

(4)

where ⟨T̂ i⟩P,λ is the average of the kinetic energy of the ith atom at a particular value of λ and for a given choice of P. In the limit that P → ∞, ΔFP → ΔF. Efficient path integral estimators for the kinetic energy have been the subject of a large body of literature.13,21,24−29 In the case of dissolved Li+, the sum over the atoms disappears since the mass of only one atom is changed. The dissolved species here is Li+ in pure water solution (Li+(aq)). For the mineral phase, first Li2O was chosen since it is a structurally simple solid, for which complete vibrational properties have been measured over the entire Brillouin zone (see SI), thus enabling a comparison of different approaches at low computational costs. Li2O has the disadvantage of not being a natural phase, and there are no solid−liquid fractionation data with which to compare. Thus, in order to validate our approach, we also computed fractionation properties of mica-type minerals: the polylithionite and a Mg-type mica model. For Li+(aq), path integral calculations at temperatures of 300, 365, 400, 500, and 600 K were carried out. Details of our free-energy convergence criteria and convergence tests are provided in the SI. The harmonic β-factor of Li+ was calculated following the approach of ref 10. Classical MD trajectories were generated at the same temperatures, from which individual configurations (“snapshots”) were extracted and quenched to 0 K, giving the so-called inherent structures (ISs). For each IS, the harmonic vibrational frequencies and harmonic β-factor were then calculated. The harmonic β-factor of the dissolved species is taken as the average of these IS β-factors. The uncertainty inherent to this method, measured as the standard error, results from the variability of the IS β-factors. We consider this approach as rigorous, as it takes into account the configurational disorder of a solution (see SI).

Figure 1. Comparison of fractionation properties (β-factors) calculated using TI-PIMD (◆) and using the harmonic approximation (HA) (straight lines), for Li2O (black) and Li+ in solution (red) as a function of temperature. For the harmonic calculations of Li+(aq), the red squares represent the average over 10 IS extracted from MD trajectories performed at the corresponding temperatures. The red straight line is based on the configurations quenched from the 365 K MD trajectory. For TI-PIMD calculations (Li2O and Li+(aq)), error bars equal to 1‰ (symbol size). Inset: Difference between the β-factors computed using the two approaches, Δanh = [ln β7Li]TI‑PIMD − [ln β7Li]harm; red dashed line, Δanh (Li+) = − 3.38

103 . T

β-factors obtained using the two approaches, HA and TI-PIMD, for Li2O and Li+(aq) are presented as a function of the inverse square temperature (1/T2). For Li2O (black symbols and lines), the two β-factors are indistinguishable, whereas for Li+(aq) (red symbols and lines), the difference is large. This emphasizes the fact that although anharmonic effects might be negligible in solid Li2O, they cannot be neglected in Li+(aq) where they lead to a significant decrease of ln β7Li. The relative importance of these anharmonic effects in Li2O and Li+(aq) is depicted in the inset of Figure 1, which plots the difference between the β-factors computed within the two approaches, Δanh ≡ [ln β7Li]TI−PIMD − [ ln β7Li]harm. 1599


Article

Accounts of Chemical Research The value of Δanh for Li+(aq) remains close to −12‰ from 300 to 400 K and then drops to approximately −4‰ at 600 K. Given the uncertainty of our calculation, we chose to describe the variation of Δanh for Li+(aq) as a function of 1/T. The magnitude of Δanh is significant compared to the natural variations in δ7Li. It demonstrates that, for Li+(aq), the harmonic calculations performed on the IS are insufficient to capture primary quantum anharmonic effects since only the configurational disorder is taken into account. Indeed, by averaging over the IS configurations, only metastable states corresponding to local minima of the potential energy surface are taken into account in the calculation of the β-factor. The significant difference between the β-factors thus obtained and those computed using TI-PIMD sheds light on the importance of exploring the free energy landscape beyond just the local minima. Use of TI-PIMD at finite temperature is, therefore, necessary to capture quantum effects that give rise to fractionation properties. Since the work of Bigeleisen31 and Urey,5 calculations of fractionation properties have overwhelmingly employed the HA. Some effort has been devoted to the development of methods for adding anharmonic corrections,32 either from experimental data33 or from atomistic calculations.34 More recently, anharmonic effects have also been evaluated on molecules through PIMD approaches.23,35−37 Among all of the isotopic systems studied, the H/D isotope fractionation exhibited the greatest degree of anharmonicity. While it was shown that anharmonic effects on H/D β-factors of individual phases are large (e.g., 150‰ for the H/D β-factor of water vapor33), the ability of these anharmonic calculations to improve the agreement with experimental isotope fractionation values has not been clearly demonstrated thus far, with the notable exception of the H/D internal equilibrium in H2O or H2S.34 One possible reason for this is that these anharmonic effects are similar in both phases and largely compensate when their β-factors are subtracted.37,38 The anharmonic effects observed here are different in several respects. First, for Li+(aq), the anharmonic contribution to the β-factor is significantly more important than for hydrogen in relative value (13% vs 6% of the harmonic β-factor at 25 °C, respectively33). Second, this anharmonic behavior is not observed in Li2O, thus justifying the use of the HA for this solid. A notable consequence is that, contrary to hydrogen, the contribution of anharmonicity to the liquid β-factor will contribute significantly to the mineral-liquid fractionation. However, comparing water vapor to Li+, we observe similar behavior regarding the variation of the anharmonic contributions to their β-factors with temperature.33 From these results, it is clear that solvated Li+ exhibits specific anharmonic effects, as evidenced by its fractionation properties.

Figure 2. Li−O pair correlation functions computed from PIMD trajectories using an empirical potential at 300 K (black lines), 365 K (red lines), 400 K (green lines), 500 K (pink lines), and 600 K (blue lines) for (a) the Li2O solid and (b) Li+ in solution.

more striking at 500 and 600 K in Li+(aq) where the Li−O distance distributions extend toward large distances. This can be attributed to water molecule exchanges taking place around Li+ in the liquid at these temperatures. Between 300 and 400 K, one can see, however, that the distributions for Li2O and Li+(aq) are very similar. This suggests that the different behavior of these two systems with respect to anharmonicity (inset of Figure 1) cannot be explained on the basis of the Li−O distance fluctuations alone. In Figure 3, we show the two-dimensional distributions of the O−Li−O angle and the Li−O distance in Li2O (top panels) and Li+(aq) (bottom panels) for all simulated temperatures. One can clearly see that the distance-angle distributions are quite different in Li2O and Li+(aq) at all temperatures. At 300 K, if the distance distributions are similar for the solid and the liquid, the angular distribution in the liquid extends over a larger domain, from 80° to almost 180°. The points appearing at large distances (∼2.8 Å) and at small (∼90°) and large (∼170°) angles correspond to configurations with 5 oxygen neighbors around the Li atom. These configurations become increasingly numerous at high temperatures (500 and 600 K) and give rise to bigger ”tails” in the large-angle region of the plots, reflecting strong deformations of the Li polyhedron. This result underlines the influence of effects beyond Li−O distance fluctuations, specifically, the deformation of the coordination polyhedron of Li in the liquid. Since important differences exist between the solid and the liquid regarding the Li coordination cage, one can assume that the anharmonic contributions to the β-factor originate mainly from the flexibility of this cage rather than unidirectional Li−O distance fluctuations. These strong deformations of the Li coordination cage disappear in the IS configurations that are used to compute the β-factor using the HA (Figure 1). Indeed, the IS are obtained from atomic relaxations to 0 K that bring the system to local minima of the potential energy surface, removing all the liquid configurations that give rise the to “tails” in the 2D-plots of Figure 3.

3.2. Local Environment around Li

Anharmonic effects likely originate from the deformation of the coordination polyhedron around Li, which is reflected in the fluctuations of Li−O distances and O−Li−O angles. By analyzing the PIMD trajectories of the two systems, Li2O and Li+(aq), it is possible to determine which of these fluctuations dominate the anharmonic contribution to the β-factors. The Li−O distance fluctuations can be correlated with the width of the first peak of the Li−O pair correlation function. Figure 2 shows these first peaks for the two systems, Li2O and Li+(aq), from the PIMD trajectories carried out at the five temperatures. These distributions exhibit an asymmetry about the peak in both systems, which reflects the intrinsic anharmonicity of the Li−O potential. This anharmonicity is even

3.3. Comparison to Experimental and Natural Data

Up to now, we have studied Li2O as a model solid for Li fractionation properties. In order to compare our approach with experimental and natural data, we consider the fractionation of Li isotopes between dissolved Li and phyllosilicates. For this family of minerals, experimental estimates and natural data are available for various temperatures, pressures, and compositions. In particular, smectite clays were considered in ref 19, while Li-micas were studied in ref 18. These two minerals are structurally alike. As shown in Figure 4, both materials consist of TOT stacking in which, for smectite, the octahedral cations are Mg2+ and Li+, whereas it is Al3+ and Li+ for polylithionite. 1600


Article


Figure 3. Contour plots of the two-dimensional distributions of the O−Li−O angle and the Li−O distance from the PIMD trajectories at all simulated temperatures. All oxygen neighbors around a Li atom within a cutoff distance of 2.9 Å were taken into account. Note that, with this definition, configurations with more than 4 oxygen neighbors around the Li are also considered. In these plots, one O−Li−O triplet gives rise to 2 points, (d1, α) and (d2, α), where d1 and d2 are the two Li−O distances in the triplet, and α is the O−Li−O angle. Upper panels, Li2O; lower panels, Li+(aq).

“all-harmonic” and “rigorous” in the subsequent discussion. Note that for the liquid computed within the HA (Figure 5a, open symbols), we performed two calculations: one with an empirical potential (squares and triangles) and one using a theoretical approach fully consistent with the solid (DFT-GGA, diamonds and line). The calculated fractionation does not significantly depend on the solid composition: the difference between Al and Mg micas is only 0.8‰ or 0.4‰ at 300 or 400 K, respectively (red squares vs red triangles). The two aforementioned approaches for Li+(aq) give very different fractionation values. In relative values, the all-harmonic calculations are in excess by amounts varying approximately between +30% at 300 K and +100% at 600 K (open vs closed red symbols). Comparing the calculated values with experimental data (Figure 5b), the all-harmonic values are significantly offset whereas very good agreement is obtained with the rigorous TI-PIMD, which is remarkable given the large diversity of conditions. The pressure in these experiments varies between 2.5 GPa18 and saturation pressure.19 In this range, experimental studies have shown a negligible effect of pressure on the Li isotope fractionation properties.45 Also, the minerals formed are of varying composition and partly disordered (see discussion in SI). In comparison, calculations are performed using a fixed density (here, 1.006 g/cm−3), and the models correspond to simple and well-defined ordered minerals. The influence of these differences on the fractionation values appears to be relatively small, in particular on the scale of the difference between rigorous path integral and all-harmonic results. Note that, due to an underestimation of the harmonic vibrational frequencies with respect to experimental values, our ab initio approach (DFT-GGA) tends to underestimate fractionation properties calculated in the HA by 5% at low temperature and 10% at high temperature.10 Taking the Li−Al mica/Li+(aq) fractionation ratio calculated using the all-harmonic method with the DFT-GGA approach for both phases (green open triangles in Figure 5), −5% of ln α = 27.3‰ corresponds to an error of −1.4‰ at 300 K. The corresponding calculation with empirical potentials for the liquid phase (red open triangles in Figure 5) is indistinguishable within the error bars. We will further assume that our description of atomic interactions leads to the same error for all calculations and that the error is within the experimental uncertainty. Finally, our calculations can be compared to natural data. These correspond to complex natural processes, such as oceanic or continental alteration, in which clay formation is expected to

Figure 4. Micas (here, polylithionite) consist of stackings of TOT layers. One TOT layer consists of two tetrahedral layers (T-layers, blue tetrahedra) separated by a layer of octahedral cations (O-layers, green Li+ and blue Al3+). In these materials, octahedral layers count two different sites: the M1 site (1 per unit formula) and the M2 site (2 per unit formula), slightly smaller in size.

To model these materials, we assumed simple setups, corresponding to a 1 M mica polymorph of composition K(Mg2Li)Si4O10(OH)2 (Mg−Li mica model) and K(Li2Al)Si4O10(OH)2 (Li−Al mica model). Within these systems, octahedral cations can occupy two different positions, denoted M1 and M2 (see Figure 4). Different possible orderings have been considered (see SI), and the structural and vibrational properties of these systems were computed using DFT (see Details). A discussion of the validity of this modeling procedure and of the structural and compositional variabilities of these compounds and their potential impact on fractionation properties is available in the SI. Figure 5 presents the results of our calculation (Figure 5a) and compares these to experimental and natural estimates of the Li fractionation factors (Figure 5b). The experimental results presented here have been selected (see SI) to correspond, as strictly as possible, to our calculations, that is, to a true isotopic equilibrium between a phyllosilicate and Li+(aq). Theoretical values have been obtained by combining the solid (Li−Al mica, red triangles; Mg−Li mica, red squares) harmonic β-factors with β-factors for the liquid computed either within the HA or using full TI-PIMD. Assuming that harmonic β-factors are rigorous for solids, we term these two approaches, respectively, as 1601


Article


Figure 5. (a) Theoretical fractionation between phyllosilicate minerals and Li+(aq). Open black symbols and solid line, harmonic fractionation between Li−Al/Mg−Li mica model and Li+(aq); here, the liquid β-factors are averaged over 10 IS, and the error bar is given by the standard error.10 Squares and triangles, Li+(aq) calculated with empirical potentials. Diamonds and line, Li+(aq) calculated with BLYP, consistently with the solid (Li−Al mica model). Red solid triangles/squares, rigorous fractionation between Li−Al/Mg−Li mica models and Li+(aq). (b) Comparison with selected experimental estimates of phyllosilicate-solution Li isotopes equilibrium (green symbols), and Li isotope fractionations associated with natural processes (blue symbols). Green squares, ref 19. Green triangles, ref 18. Green circle, ref 39. Blue diamonds, refs 40−42. Blue triangle, ref 43. Blue square, ref 44. Blue circle, ref 39.

a short oxygen residence time is likely to be surrounded by a loose solvation shell and to show strong thermovibrational anharmonicity. The contrast between Mg2+(aq), which shows no thermovibrational anharmonicity and a relatively long residence time51 (τres ≈ 10−5 to 10−6 s), and Li+(aq), which shows a strong vibrational anharmonicity together with a short residence time51 (τres ≈ 10−11 to 10−12 s) is consistent with this idea. In this context, cations with short residence times, such as alkali and alkaline earth metals, Cu+/2+, Cr2+, Cd2+, or Hg2+, are expected to show particularly large thermovibrational anharmonicity. Extending this idea to dissolved halogen anionic species, such as F−, Cl−, Br−, or I−, which show particularly rapid water exchange, these species should also exhibit significant anharmonicity. On the other hand, dissolved species such as H4SiO4 or B(OH)3 can be considered essentially as nondissociating. For these species, the HA should be valid. In order to confirm these assumptions, classical and path-integral MD could be carried out for these systems. In this work, quantum anharmonicity exhibited by dissolved species has been demonstrated for the isotopic equilibrium between a solid and a dissolved species. Another important case, which is the subject of many studies, is the isotopic equilibrium between two coexisting species in solution. One might wonder about the extent to which anharmonic effects on the β-factors of each species might cancel when computing fractionation factors. Such a cancellation would be unlikely when considering two different species with coordination shells of differing degrees of rigidity. Finally, the overall agreement of our theoretical curve with experimental results taken in various conditions and with natural data and the small difference between our different mineral models strongly support a negligible impact on Li isotope fractionation of the structure, composition or cation ordering of the mineral phase (clay19 or mica), of water composition or pressure (see also discussion in SI). This agreement supports the use of many clay isotope compositions for paleoenvironmental reconstruction, for example.

play an important role. Indeed these data show good agreement with our rigorous calculations.

4. DISCUSSION Previously calculations of the fractionation properties of dissolved species were mostly based on the HA. Our results indicate that this approach is inaccurate for Li+(aq). The question then arises as to whether the failure of the HA can be generalized to all solvated species. The only way to answer this question definitively is to perform rigorous nuclear quantum simulations via path integration for each solvation situation. There have previously been a number of studies using nuclear PI methods to investigate the structural properties of dissolved species including F− and Li+, Cl−, I−, H3O+, OH−.46−50 However, to our knowledge, applications to isotopic fractionation have been limited, thus far, to Mg2+(aq),12 where it was reported that the HA method retains validity provided that the liquid configurational disorder is taken into account. The contrast between Li+(aq) and Mg2+(aq) might be due to differences in the flexibility of their respective solvation shells. Along this line of reasoning, the analysis of the local environment of the dissolved species could give a first indication of whether anharmonic contributions to the thermovibrational properties are significant. This could be done via a simple structural analysis of MD trajectories. Large fluctuations of the coordination polyhedra around the solute could be a first indication of significant anharmonic contributions to the thermovibrational properties. For instance, analysis of MD trajectories of H4SiO4 in water10 reveals very limited angular fluctuations of the Si coordination polyhedra at 300 K (109 ± 3°), suggesting that full PI calculations may not be necessary. As a first indication of the flexibility of coordination polyhedra, the residence time τres of water molecules in the first hydration shell of dissolved ions51 can enable a comparison between different ionic systems and suggest which species are likely to exhibit strong vibrational anharmonicity. Indeed, an ion with 1602


Article


5. DETAILS Classical MD and PIMD simulations were performed with the PINY_MD package52 using empirical potentials. The different parameters used in these simulations can be found in SI. For Li2O, a 2 × 2 × 2 supercell of the Li2O fcc lattice (of nominal composition Li64O32) was employed. Simulations were performed using empirical potentials of the Buckingham type, which were obtained by the authors of ref 53, to fit to the experimental vibrational modes. Accuracy of these potentials with respect to experimental data are shown in SI. For Li+(aq), simulations were performed on a system of 62 water molecules, 1 Li+ cation, and 1 Cl− anion at a density of 1.006 g/cm−3. Since these properties are very local, a small box of 62 water molecules appeared to be sufficient. The q-SPC/fw model54 was used to treat the water molecules, and the model of ref 55 was chosen for the Li−water interactions. The performance of the potentials for Li+(aq) was evaluated with respect to ab initio MD (AIMD) simulations at the DFT-BLYP level of electronic structure theory (see SI and Figure 5). To compute fractionation properties in the HA, both for the mica-type minerals and for Li+(aq), DFT calculations were employed. The calculation of the fractionation properties was based on the phonon modes of quenched configurations extracted from AIMD in the case of the liquid, following the standard procedure.10,56 The calculation was based on the BLYP functional,57,58 a plane-wave basis set, and atomic pseudopotentials as implemented in the Quantum Espresso package. Pseudopotentials used for Al, Si, O, and H are described in ref 10. The K, Li, and Mg pseudopotentials were taken from the PSlibrary.59 Planewave cutoffs of 80 and 320 Ry were used for the electronic orbitals and density, respectively.

■

strongly interested in isotopic fractionation and its relation to the structural properties of matter. Magali Benoit obtained her B.S. in particle physics from the Université Paris XI in 1993 and her Ph.D. in physics from the Université Paris XI in 1996. During her Ph.D. research, she worked at the IBM Forschungslaboratorium in Rüschlikon, Switzerland, for 6 months, at the Max-Planck-Institut für Festkörperforschung in Stuttgart, Germany, for 1 year, and at the Laboratoire des Verres in Montpellier, France, for 18 months. From 1997 to 2006, she was a permanent CNRS research fellow at the Laboratoire des Verres in Montpellier, France. In 1999, she obtained a NATO grant for a four-month visit to the Chemistry Department at New York University. Since 2006, she has been a permanent CNRS research fellow in the Centre d’Elaboration des Matériaux et d’Etudes Structurales (CEMES) in Toulouse, France. Her research interests include investigation of nuclear quantum effects on the transitions between high-pressure phases (ice and hydrous minerals), elucidation of the relationship between structural and thermodynamic properties in silicate glasses and melts, exploration of the effects of growth and morphology on the functionalization of metallic nanoparticles, and study of isotopic fractionation between liquids and minerals. Mark Tuckerman obtained his B.S. in physics from U.C. Berkeley in 1986 and his Ph.D. in physics from Columbia University in 1993. From 1993 to 1994, he was an IBM postdoctoral fellow at the IBM Forschungslaboratorium in Rüschlikon, Switzerland, and from 1995 to 1996, he held an NSF Postdoctoral Fellowship in Advanced Scientific Computing at the University of Pennsylvania. He is currently Professor of Chemistry and Mathematics at New York University and is affiliated with the NYU-ECNU Center for Computational Chemistry at NYU Shanghai. He is the recipient of an NSF CAREER award, an NYU Whitehead Fellowship in Biomedical and Biological Sciences, an NYU Golden Dozen Award for Excellence in Teaching, the Camille Dreyfus Teacher-Scholar award, the Friedrich Wilhelm Bessel Research Award from the Alexander von Humboldt Foundation, and a Japan Society for the Promotion of Science Research Fellowship. His research interests include proton and hydroxide transport phenomena, crystal structure prediction, the development of enhanced configurational sampling and free energy prediction algorithms for oligopeptide conformations and phase transitions, and explorations of nuclear quantum effects.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.accounts.6b00607. Molecular dynamics and path integral molecular dynamics simulation parameters, accuracy of the empirical potentials for Li2O and Li+(aq), and description of the different phyllosilicate models and discussion on their relevance for comparison with experimental and natural data (PDF)

■

Merlin Méheut is an Assistant Professor of Mineralogy at Paul Sabatier University, Toulouse (France). He obtained his Ph.D. in solid state physics from University Pierre et Marie Curie, Paris (France), in 2008. From 2008 to 2009, he carried out postdoctoral research in the Department of Earth and Space Sciences at the University of California, Los Angeles. His research focuses on the exploration of geochemical problems using methods of theoretical solid state physics (statistical physics and quantum modeling of materials). He is particularly interested in the understanding of natural isotope separation processes at atomic and molecular scales. His recent research has focused on Si and Li isotope fractionation between rock-forming minerals and aqueous solutions. He is also investigating the connection between mineral structures and compositions and isotope fractionation properties.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Romain Dupuis: 0000-0001-9451-1132 Mark E. Tuckerman: 0000-0003-2194-9955 Merlin Méheut: 0000-0003-2037-2365 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Calculations were performed using HPC resources from CALMIP (Grant 2015-P1037) and CINES (Grant c2014096831). R.D. and M.M. acknowledge support from the interdisciplinary program PEPS Theoretical Physics and Its Interfaces of French CNRS. M.E.T. acknowledges support from the National Science Foundation Grants CHE-1566085 and CHE-1565980.

Biographies Dr. Romain Dupuis is a young scientist who earned his Ph.D. in 2014 at the Université Paul Sabatier Toulouse III in France. Since then, he has been working at the Donostia International Physics Center as a postdoctoral fellow studying thermodynamical properties of cementbased materials. His research is dedicated to the simulation of solids and liquids at the atomic scale by the means of molecular dynamics. He is 1603


Article


■

(22) Schmitt, A.-D.; Vigier, N.; Lemarchand, D.; Millot, R.; Stille, P.; Chabaux, F. Processes controlling the stable isotope compositions of Li, B, Mg and Ca in plants, soils and waters: A review. C. R. Geosci. 2012, 344, 704−722 Erosion-Alteration: from fundamental mechanisms to geodynamic consequences (Ebelmen’s Symposium).. (23) Ceriotti, M.; Markland, T. E. Efficient methods and practical guidelines for simulating isotope effects. J. Chem. Phys. 2013, 138, 014112. (24) Parrinello, M.; Rahman, A. Study of an F center in molten KCl. J. Chem. Phys. 1984, 80, 860. (25) Herman, M. F.; Bruskin, E. J.; Berne, B. J. On path integral Monte Carlo simulations. J. Chem. Phys. 1982, 76, 5150. (26) 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. (27) Jang, S.; Jang, S.; Voth, G. Applocations of higher order composite factorization schemes in imaginary time path integral simulations. J. Chem. Phys. 2001, 115, 7832−7842. (28) Takahashi, M.; Imada, M. Monte Carlo Calculation of Quantum Systems. II. Higher Order Correction. J. Phys. Soc. Jpn. 1984, 53, 3765− 3769. (29) Li, X.-P.; Broughton, J. Q. High-order correction to the Trotter expansion for use in computer simulation. J. Chem. Phys. 1987, 86, 5094. (30) Criss, R. E. Temperature dependance of isotopic fractionation factors. In Stable isotope geochemistry: a tribute to Samuel Epstein; H. P. Taylor, J., O’Neil, J. R., Kaplan, I. R., Eds.; The Geochemical Society, 1991; pp 11−16. (31) Bigeleisen, J.; Mayer, M. G. Calculation of Equilibrium Constants for Isotopic Exchange Reactions. J. Chem. Phys. 1947, 15, 261−267. (32) Wolfsberg, M. Correction to the effect of anharmonicity on isotopic exchange equilibria-application to polyatomic molecules. J. Chem. Phys. 1969, 50, 1484. (33) Richet, P.; Bottinga, Y.; Javoy, M. A review of hydrogen, carbon, nitrogen, oxygen, and chlorine stable isotope fractionation among gaseous molecules. Annu. Rev. Earth Planet. Sci. 1977, 5, 65−110. (34) Liu, Q.; Tossell, J. A.; Liu, Y. On the proper use of the BigeleisenMayer equation and corrections to it in the calculation of isotopic fractionation equilibrium constants. Geochim. Cosmochim. Acta 2010, 74, 6965−6983. (35) Mielke, S. L.; Truhlar, D. G. Improved Methods for Feynman 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. (36) 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. (37) Webb, M. A.; Miller, T. F. Position-Specific and Clumped Stable Isotope Studies: Comparison of the Urey and Path-Integral Approaches for Carbon Dioxide, Nitrous Oxide, Methane, and Propane. J. Phys. Chem. A 2014, 118, 467−474. (38) Méheut, M.; Lazzeri, M.; Balan, E.; Mauri, F. First-principles calculation of H/D isotopic fractionation between hydrous minerals and water. Geochim. Cosmochim. Acta 2010, 74, 3874−3882. (39) Zhang, L.; Chan, L.-H.; Gieskes, J. M. Lithium isotope geochemistry of pore waters from ocean drilling program Sites 918 and 919, Irminger Basin. Geochim. Cosmochim. Acta 1998, 62, 2437− 2450. (40) Chan, L.; Edmond, J.; Thompson, G.; Gillis, K. Lithium isotopic composition of submarine basalts: implications for the lithium cycle in the oceans. Earth Planet. Sci. Lett. 1992, 108, 151−160. (41) Chan, L.-H.; Edmond, J. M.; Thompson, G. A Lithium Isotope Study of Hot Springs and Metabasalts From Mid-Ocean Ridge Hydrothermal Systems. J. Geophys. Res. 1993, 98 (B6), 9653−9659. (42) Lui-Heung, C.; Gieskes, J. M.; Chen-Feng, Y.; Edmond, J. M. Lithium isotope geochemistry of sediments and hydrothermal fluids of the Guaymas Basin, Gulf of California. Geochim. Cosmochim. Acta 1994, 58, 4443−4454.

REFERENCES

(1) Hoefs, J. Stable isotope geochemistry; Springer Science & Business Media: Berlin Heidelberg, 2008. (2) Johnson, C. M., Beard, B. L., Albarede, F., Eds. Geochemistry of nontraditional stable isotopes; Reviews in Mineraloly Vol 55; Mineralogical Society of America, Geochemical Society, 2004. (3) Black, J. R.; Umeda, G.; Dunn, B.; McDonough, W. F.; Kavner, A. Electrochemical Isotope Effect and Lithium Isotope Separation. J. Am. Chem. Soc. 2009, 131, 9904−9905. (4) Telouk, P.; Puisieux, A.; Fujii, T.; Balter, V.; Bondanese, V. P.; Morel, A.-P.; Clapisson, G.; Lamboux, A.; Albarede, F. Copper isotope effect in serum of cancer patients. A pilot study. Metallomics 2015, 7, 299−308. (5) Urey, H. C. The thermodynamic properties of isotopic substances. J. Chem. Soc. 1947, 562−581. (6) Yamaji, K.; Makita, Y.; Watanabe, H.; Sonoda, A.; Kanoh, H.; Hirotsu, T.; Ooi, K. Theoretical Estimation of Lithium Isotopic Reduced Partition Function Ratio for Lithium Ions in Aqueous Solution. J. Phys. Chem. A 2001, 105, 602−613. (7) Fujii, T.; Moynier, F.; Blichert-Toft, J.; Albaréde, F. Density functional theory estimation of isotope fractionation of Fe, Ni, Cu, and Zn among species relevant to geochemical and biological environments. Geochim. Cosmochim. Acta 2014, 140, 553−576. (8) Schauble, E. A. First-principles estimates of equilibrium magnesium isotope fractionation in silicate, oxide, carbonate and hexaaquamagnesium(2+) crystals. Geochim. Cosmochim. Acta 2011, 75, 844−869. (9) Kowalski, P. M.; Jahn, S. Prediction of equilibrium Li isotope fractionation between minerals and aqueous solutions at high P and T: An efficient ab initio approach. Geochim. Cosmochim. Acta 2011, 75, 6112−6123. (10) Dupuis, R.; Benoit, M.; Nardin, E.; Méheut, M. Fractionation of silicon isotopes in liquids: The importance of configurational disorder. Chem. Geol. 2015, 396, 239−254. (11) Rustad, J. In Annual Reports in Computational Chemistry; Dixon, D. A., Ed.; Elsevier: Amsterdam, 2016; Vol. 12, pp 117−156. (12) Pinilla, C.; Blanchard, M.; Balan, E.; Natarajan, S. K.; Vuilleumier, R.; Mauri, F. Equilibrium magnesium isotope fractionation between aqueous Mg2+ and carbonate minerals: Insights from path integral molecular dynamics. Geochim. Cosmochim. Acta 2015, 163, 126−139. (13) Marsalek, O.; Chen, P.-Y.; Dupuis, R.; Benoit, M.; Méheut, M.; Bačić, Z.; Tuckerman, M. E. Efficient Calculation of Free Energy Differences Associated with Isotopic Substitution Using Path-Integral Molecular Dynamics. J. Chem. Theory Comput. 2014, 10, 1440−1453. (14) 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. (15) Raymo, M.; Ruddiman, W.; Froelich, P. Influence of late cenozoix mountain building on ocean geochemical cycles. Geology 1988, 16, 649− 653. (16) Lui-Heung, C.; Edmond, J. M. Variation of lithium isotope composition in the marine environment: A preliminary report. Geochim. Cosmochim. Acta 1988, 52, 1711−1717. (17) Tang, Y.-J.; Zhang, H.-F.; Ying, J.-F. Review of the Lithium Isotope System as a Geochemical Tracer. Int. Geol. Rev. 2007, 49, 374− 388. (18) Wunder, B.; Meixner, A.; Romer, R. L.; Feenstra, A.; Schettler, G.; Heinrich, W. Lithium isotope fractionation between Li-bearing staurolite, Li-mica and aqueous fluids: An experimental study. Chem. Geol. 2007, 238, 277−290. (19) Vigier, N.; Decarreau, A.; Millot, R.; Carignan, J.; Petit, S.; FranceLanord, C. Quantifying Li isotope fractionation during smectite formation and implications for the Li cycle. Geochim. Cosmochim. Acta 2008, 72, 780−792. (20) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill Companies: New York, 1965. (21) Tuckerman, M. E. Statistical Mechanics: Theory and Molecular Simulation; Oxford University Press: Oxford, 2010. 1604


Article

Accounts of Chemical Research (43) Millot, R.; Scaillet, B.; Sanjuan, B. Lithium isotopes in island arc geothermal systems: Guadeloupe, Martinique (French West Indies) and experimental approach. Geochim. Cosmochim. Acta 2010, 74, 1852− 1871. (44) Huh, Y.; Chan, L.-H.; Zhang, L.; Edmond, J. M. Lithium and its isotopes in major world rivers: implications for weathering and the oceanic budget. Geochim. Cosmochim. Acta 1998, 62, 2039−2051. (45) Wunder, B.; Meixner, A.; Romer, R. L.; Jahn, S. Li-isotope fractionation between silicates and fluids: Pressure dependence and influence of the bonding environment. Eur. J. Mineral. 2011, 23, 333− 342. (46) Wilkins, D. M.; Manolopoulos, D. E.; Dang, L. X. Nuclear quantum effects in water exchange around lithium and fluoride ions. J. Chem. Phys. 2015, 142, 064509. (47) Schenter, G. K.; Garrett, B. C.; Voth, G. A. The quantum vibrational dynamics of Cl-(H2O) n clusters. J. Chem. Phys. 2000, 113, 5171−5178. (48) Gai, H.; Schenter, G. K.; Dang, L. X.; Garrett, B. C. Quantum statistical mechanical simulation of the ion-water cluster I-(H2O) n: The importance of nuclear quantum effects and anharmonicity. J. Chem. Phys. 1996, 105, 8835−8841. (49) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. The nature of the hydrated excess proton in water. Nature 1999, 397, 601−604. (50) Tuckerman, M. E.; Marx, D.; Parrinello, M. The nature and transport mechanism of hydrated hydroxide ions in aqueous solution. Nature 2002, 417, 925−929. (51) Ohtaki, H.; Radnai, T. Structure and Dynamics of Hydrated Ions. Chem. Rev. 1993, 93, 1157. (52) Tuckerman, M.; Yarne, D.; Samuelson, S.; Hughes, A.; Martyna, G. Exploiting multiple levels of parallelism in Molecular Dynamics based calculations via modern techniques and software paradigms on distributed memory computers. Comput. Phys. Commun. 2000, 128, 333. (53) Goel, P.; Choudhury, N.; Chaplot, S. Lattice dynamics of lithium oxide. Pramana 2004, 63, 409−412. (54) Paesani, F.; Zhang, W.; Case, D. A.; Cheatham, T. E.; Voth, G. A. An accurate and simple quantum model for liquid water. J. Chem. Phys. 2006, 125, 184507. (55) Lyubartsev, A. P.; Laasonen, K.; Laaksonen, A. Hydration of Li+ ion. An ab initio molecular dynamics simulation. J. Chem. Phys. 2001, 114, 3120−3126. (56) Méheut, M.; Lazzeri, M.; Balan, E.; Mauri, F. Equilibrium isotopic fractionation in the kaolinite, quartz, water system: predictions from first-principles density-functional theory. Geochim. Cosmochim. Acta 2007, 71, 3170−3181. (57) Becke, A. Density-Functional exchange-energy approximation with correct assymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098. (58) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785. (59) Dal Corso, A. Pseudopotentials periodic table: From H to Pu. Comput. Mater. Sci. 2014, 95, 337−350.

1605


Importance of a Fully Anharmonic Treatment of Equilibrium Isotope

Recommend Documents