Computational Investigations of Position-Specific Vapor Pressure Isotope Effects in Ethanol—Toward More Powerful Isotope Models for Food Forensics

With the advent of new experimental techniques, measurements of individual, per-position, vapor pressure isotope effects (VPIEs) became possible. Frequently, they are in opposite directions (larger and smaller than unity), leading to the cancellation when only bulk values are determined. This progress has not been yet paralleled by the theoretical description of phase change processes that would allow for computational prediction of the values of these isotope effects. Herein, we present the first computational protocol that allowed us to predict carbon VPIEs for ethanol—the molecule of great importance in authentication protocols that rely on the precise information about position-specific isotopic composition. Only the model comprising explicit treatment of the surrounding first-shell molecules provided good agreement with the measured values of isotope effects. Additionally, we find that the internal vibrations of molecules of the model to predict isotope effects work better than the entire set of normal modes of the system.


■ INTRODUCTION
Ethanol isotopic composition is used in the authentication of alcoholic beverages using isotope-ratio mass spectroscopy (IRMS) 1 for carbon and oxygen isotopic composition and stable natural isotope fractionation (SNIF-NMR) 2 for hydrogen isotopic composition. The determination of carbon isotopic composition, both bulk and position-specific (i.e., isotopic composition of each at each of carbon atom positions separately), requires the isolation of pure ethanol in a manner that preserves isotopic composition. This is routinely achieved nowadays by spinning-band distillation. It is assumed that under these conditions, the isotopic fractionation is minimized to reflect only the vapor pressure isotope effect (VPIE) of the ethanol phase transfer from liquid to gaseous form, and sufficient recovery of ethanol should therefore allow us to refer the isotopic composition of the distillate to that of the original sample (wine or brandy with a yield of more than 96%). 3,4 However, neither determination nor interpretation of the VPIE of ethanol is straightforward.
The first theoretical assessment of the phase transition isotope effects was formulated by Lindemann, 5 followed by the work of Bigeleisen and Stern, Van Hook, and Wolfsberg. These early calculations of VPIEs used the harmonic approximation method. 6−8 Currently, well-established theories of VPIE are based on the existence of liquid and gas films at the border of phases that are at equilibrium at the interface. They provide the most complete explanation of isotope fractionation during the phase transition. Despite the lack of protocols allowing for the complete removal of the distillate, initial vaporization, an impact of theoretical plates, and difficulties with maintaining the system at perfect equilibration, distillation under rigorously controlled conditions is nowadays considered the best approach to study liquid−vapor isotope effects.
The development and improvement of accuracy and precision in stable carbon isotopic measurements at a specific position in a compound 9,10 allowed us to demonstrate that ethanol shows position-specific 13 C behavior during distillation and evaporation phenomenon processes. 11−13 Nevertheless, a more profound understanding of site-specific VPIE due to many intermolecular interactions in the liquid or condensed state is still a very challenging task. 14 Botosoa and co-workers 11 claim that the presence of 13 C on the methylenic carbon in ethanol has an impact on the energy amount needed to break noncovalent interactions. Isotope fractionation during volatilization as a result of VPIEs was reported to date to address the issue of environmental isotope study of contaminants in order to investigate their fate in groundwater or open water surfaces. As a result of these findings, normal (greater than unity) carbon isotope effects occurred in both pure-phase evaporation and volatilization of dissolved ethanol in water, which indicates preferential volatilization of lighter isotopes. 8,12,13 Table 1 summarizes the experimental data of VPIEs.
At this point, a comment on the terms and definitions is necessary as they are not uniformly used in the literature, especially across disciplines such as geochemistry, chemistry, and environmental chemistry. 15 From the physicochemical point of view, for an equilibrium process of substance S going from a condense, liquid phase (l) to a gaseous phase (g) presented in eq 1 the equilibrium isotope effect, in this case called VPIE, is given by the ratio of equilibrium constants, K L /K H , of light (L) and heavy (H) species. It can be obtained experimentally from the so-called isotopic ratio, R, which is defined as the ratio of the amount (e.g., concentration) of heavy species to the amount of light species, R = [H]/[L]. 6 Because K = [g]/[l], thus VPIE = K L /K H = R l /R g for isotopic compositions measured at equilibrium. This value has been also termed fractionation factor 4 or equilibrium isotopic fractionation factor 14 and marked as α. Unfortunately, the same symbol is used by very active in the field French groups to denote the so-called separation factor, 13,11,16 which is exactly inverse of the isotope effect (this is because they use inverse definition of the isotopic ratio). Sometimes, while the same definition as above of the isotopic ratio is used, α is reported as the ratio of equilibrium constants of gas to liquid leading again to the inversion of the reported values. 17 Because VPIEs are very small, they are usually presented as "enrichment factors", ε, which correspond to the exponent in the Rayleigh equation, 15 which originates from the mathematical formalism of the kinetic isotope effect on the first-order reaction. The enrichment factor, ε = 1/(K L /K H ) − 1, is traditionally reported in the so-called δ-notation, that is, after multiplication of the above expression by 1000, and termed "per-mil" [‰] (equivalent to mUr in the SI metric system).
To add to the confusion, the term ε is sometimes called directly an isotope effect. 18 The above discussion shows how confusing following the literature on phase-transfer isotope effect can be and how cautious the reader has to be while comparing results reported by different laboratories. In this contribution, we apply the physicochemical definition of isotope effects and isotopic ratios. Thus, normal isotope effect means the value is larger than unity, and the corresponding ε is negative. On the contrary, an isotope effect is called inverse, when it is smaller than unity, and the corresponding ε is positive.
The values of 13 C-VPIE obtained using distillation, reported by different laboratories, are consistently inverse (less than unity). We have, therefore, calculated enrichment factors, ε, based on the averaged values. The primary objective of this study was to give a proof-of-concept for using methods of theoretical chemistry to study ethanol VPIEs. Based on the new theoretical calculation results, the impact of using different models for calculating VPIEs on ethanol is discussed.

■ RESULTS
In the present work, we aimed at predicting experimentally observed inverse carbon isotopic fractionation during the ethanol distillation (Table 1). Additionally, we challenged selected computational models by making an attempt to capture isotopic distinction between two carbon positions present in the ethanol molecule. To this end, three types of models have been employed in the studies. The simplest included a single molecule of ethanol in each phase with liquid represented using the polarized continuum model (PCM). All these calculations were performed using the ωB97xD functional expressed in split valence Dunning or def2-family basis sets. Three variants of the PCM solvent model have been used as described in the Methods section. The single molecule models were then enlarged by adding ethanol molecules to the models mimicking liquid phase only and by including more    Table 2.
Because it was not possible to treat more complex models involving larger numbers of ethanol molecules computationally at the density functional theory (DFT) level (cost of full vibrational analysis necessary for the frequency and subsequent isotope effects calculations), we prepared larger spherical models starting from a box of 1165 ethanol molecules and containing 91 and 138 molecules. Next, we divided them into two or three parts (layers, following the terminology of ONIOM QM/MM calculations used here). The example of such a model is shown in Figure 1. The layer to which the highest theory level has been applied is labeled H (high) and contained from 1 to 8 molecules. For the description of the next layer (if present), a smaller basis set has been used. The remaining molecules were treated with molecular mechanics, mainly Amber (parm96.dat 21 as implemented in Gaussian) force fields. Dreiding force field 22 has been used to investigate the influence of the type of the force field, but only minimal differences have been obtained (see entries in Table 3).
The third group of models constitutes cluster models comprising 4, 8, or 17 ethanol molecules (exemplified in Figure 2). These models were prepared based on the ONIOM models as described in detail in the Methods section. Geometries of these cluster models were optimized in the gas phase. In one case, the IEF continuum solvent model has been included for comparison. That cluster for which Gibbs free energy rather than electronic energy was optimized yielded geometry very close to those optimized in the gas phase; therefore, this geometry was assumed as representative for the cluster. The results are collected in Table 4.

■ DISCUSSION
The results collected in Table 2 might lead to a general observation that inclusion of the reaction field can guarantee the right direction of isotopic fractionation (inverse effects, ε values greater than zero), however, not in all studied models. An agreement between experiment and theory is only qualitative, and in most of the cases, the CH 3 position-specific isotope effect is larger than the one at the CH 2 position, which is in contrast to the trend observed experimentally. Quite surprisingly, the agreement with the measured direction of isotope effects was achieved only when the liquid phase was represented by a single ethanol molecule immersed in the continuum solvent model. Enlarging the model to eight molecules resulted in the wrong direction of CH 3 position isotope effects (normal instead of inverse). The only exception was seen when both phases, liquid and gas, were mimicked by the same size model (eight molecules). The problem, in this case, was, however, the lack of distinction between the CH 3 and CH 2 carbon position isotope effects. Moreover, we found that using internal normal modes (of molecules from a cluster mimicking the liquid phase of the system) in the calculation of isotope effects is by far better than taking into account the entire (internal and collective) set of frequencies (for the whole cluster comprising eight molecules). Only minor temperature dependence of the calculated VPIE in the range of 25−78°C has been observed (see the first three entries in Table 2); therefore, all further isotope effects discussed in this work were calculated only at 25°C.
All hybrid protocols (QM/MM or QM/QM/MM) tested in this work failed to predict isotope effects even qualitatively.   Except one, all of them resulted in normal instead of inverse isotope effects. Part of the reason behind this poor performance of hybrid schemes may lie in inadequate description of ethanol by force fields that were used for a low layer of the system, which, in each model, comprised most of the overall model. The proper description of organic liquids by selected force fields and its consequence in magnitudes of predicted isotope effects have been also addressed in recent studies by some of us. 24,25 We have shown that substantial differences in predicted isotope effects can be found if a force field takes into account some important characteristic properties of the liquid from the interaction point of view. Furthermore, the failure of QM/MM approach might also come from overpolarization of the QM part, especially when the QM part is not large enough. 26 However, it does not seem to be the case of the models presented here as the atomic charges investigation did not provide any significant changes in the charges of high and medium layers (both represented by the QM levels of theory), and additionally, the atomic charges between the mechanical and the electronic embedding approaches almost did not differ at all. Substantial improvement of the predicted VPIEs was observed while using the largest cluster (17 molecules). Although it is not the quantitative agreement yet, the distinction in the position-specific isotope effects between methyl and methylene carbon atoms is captured (the 13 CH 3 VPIE is smaller than the 13 CH 2 VPIE). Size-wise, the 17molecule cluster is twice as large as the eight-molecule cluster that may have its consequence in apparently better architecture of the first solvation shell of the central molecule in the former model which evidently affects the predicted magnitudes and trends of VPIEs. In the case of the four-molecule cluster, it is even impossible to discuss the issue of a central molecule in the system. In the eight-molecule cluster, one could possibly find the one but its surrounding would not guarantee sufficient first solvation shell description ( Figure 2). Only the 17-molecule cluster provides adequate improvement. Very recent studies by some of us 41 and others 27 demonstrate that carbon−carbon and carbon−oxygen noncovalent interactions are important and should not be neglected when a proper description of liquid ethanol behavior is of interest. This particular feature was observed, regardless the level of theory used along with molecular dynamics simulations, and was with quite reasonable agreement with experimentally obtained radial distribution function using neutron scattering. 28 Ethanol, a linear monohydroxy alcohol, is a strongly hydrogen-bonded liquid; therefore, one might expect that the geometry as well as strength of these interactions might play a crucial role in isotopic fractionation observed during its evaporation, in particular, in specificity for distinct positions as previously demonstrated. 29 Because the ONIOM and cluster models used in this work resulted in two completely different prediction patterns, we used these two, the largest, models for comparative study with respect to the key hydrogen-bonding parameters such as its number in the model, the distance between the donor and the acceptor oxygen atom (r O···O ) as well as the OH···O angle. We found that although much less HBs are formed in the ONIOM model, there are basically no geometrical differences in the distance and the angle as the average r O···O is 2.84 ± 0.02 and 2.74 ± 0.06 Å and <OHO is 175 ± 3 and 171 ± 3°for the ONIOM and cluster model, respectively. Yet these two provide substantially different predictions of isotope effect (−0.37 and −0.97 and 1.08 and 2.33 for 13 CH 3 and 13 CH 2 VPIE, respectively). It is possible that a number of HBs present in a model may play a role; however, we do not have sufficient data to test it. Therefore, the reasons standing behind these differences must be sought somewhere else.
For cluster models, we also observe that calculating isotope effects based only on internal frequencies turn out to be beneficial, and it is noticed even for the smallest cluster treated at the ωB97xD/aug-cc-pVTZ level of theory (rows 2 and 3 in Table 4). In the case of the eight-molecule cluster, we also obtained the values of isotope effects more positive, thus closer   Table 4). The applicability of a Hessian calculated only for few atoms out of a larger system being isotopically substituted for the correct prediction of isotope effects has been addressed in a very recent study on the effect of dielectric environments on isotope effects. 30 Sẃiderek et al. in their work have demonstrated that a small set of atoms (a moiety isotopically substituted) instead of the entire system can be used for isotope effect prediction because as long as the influence of the whole environment on this small set is reflected from the IEs obtained from these two approaches are basically the same. This is in contrast to our observation of improved prediction when only a subset of total number of normal modes of a system is used for isotope effect prediction. What is different in our approach, and to the best of our knowledge has not been considered before, is the way how the selection of frequencies to calculate isotope effects is performed. In our work, we did not use "cut-off" procedures proposed previously 31−33 within which a series of smaller subset Hessians are generated by deleting rows and columns corresponding to particular atoms and are subsequently used for calculating the isotope effect. In this work, after computing full Hessian for each system and knowing the number of normal modes coming from all ethanol molecules in a model but treated separately, we calculated the number of collective modes. It is 6n − 6, where n is the number of ethanol molecules in a model. In order to obtain only internal normal modes, we projected 6n − 6 set of collective modes out of the full 3N − 6 set of modes for each analyzed model. Including all 3N − 6 set in the calculation of isotope effects turned out to give worse agreement with the experimental data, whereas using only internal modes proved to be superior (rows 9 and 10 in Table 4). Technically, all modes below 250 cm −1 were not included in the isotope effect computation. The predicted magnitudes within this particular model did improve even more by introducing additional continuum solvation (rows 10 and 11 in Table 4). In order to further explore the issue of inclusion or exclusion of various modes as well as distribution of isotope effects all over the molecules present in the largest cluster in this study, we performed additional analysis in which VPIEs were calculated for all CH 3 and all CH 2 carbon positions and then averaged (Table S1 in Supporting Information). We considered three scenarios. In the first one, all modes were included in the calculation of isotope effects; six translations and rotations were also taken into account (3N = 27n). In the second one, the six vibrational modes were projected out (3N − 6 = 27n − 6); however, it should be noted that they all were very close to zero. In the third one, only internal modes of each molecule in the 17molecule clusters were included (21n). The first striking observation is that the averaged magnitudes of VPIEs are substantially different from the values obtained for the central molecule. In the case of 3N and 3N − 6 approaches, calculating VPIEs based on average isotopic response coming from all molecules results in the change in the direction of position-specific isotope effects as compared to the values obtained for the central molecule only. The 27n case is different, though. The direction is not changed, and for both ways, we see inverse effects only as the relation between the two carbon positions changes to the one which does not agree with the experimental finding. Additionally, we calculated the variations in the isotope effect over the molecules present in the model. They turned out to be quite substantial. In the case of 3N and 3N − 6 models, their inclusion would lead to the change of the isotope effect direction. Only the 27n model does not seem to suffer from it. The results for the 17-molecule cluster also provide an interesting example of the positive effect of Def2 family of basis set, in particular the effect of its size. Def2-TZVP has been shown recently to perform well in calculations of isotopic fractionations. 34

■ CONCLUSIONS
In this work, we aimed at finding a theoretical level to model VPIEs of ethanol. To this end, we tested distinct theoretical models to describe the liquid phase such as implicit, explicit, and hybrid solvation models and different computational techniques including ONIOM protocol within the QM/MM method and the QM cluster approach.
Within the present study we demonstrated that in order to be on a right track with the predicted VPIEs, one has to take into account several factors: (i) proper solvation described using (ii) proper level of theory and (iii) reflecting adequate treatment of the first solvation shell (explicit representation is preferable and sometimes mandatory, (iv) selection of frequencies which are used for isotope effects computation surprisingly internal set provides better agreement with the experiment than the full set of frequencies. Furthermore, using larger models, where all molecules were equally treated in terms of the theory level applied to optimize the coordinates of the system and to perform normal mode analysis, may result in different sets of predicted isotope effects depending only on whether a central molecule or all molecules within this model are the source of isotopic substitution. For models having sufficient surroundings by other neighboring molecules, such as the 17-molecule cluster, the former approach seems to be much a better choice. Collecting isotopic fractionation from all molecules in a still reduced model or using all normal modes to calculate bulk VPIE results in large numerical noise in the predicted isotope effects which can be difficult (even impossible) to reduce because the VPIEs themselves represent very small deviations from unity. On the other hand, larger models might turn out to be too expensive to be treated quantum mechanically.
Based on the calculations and analysis carried out, we were able to provide, for the first time, the correct result in terms of direction (inverse) and distinction between the two carbon position 13 C VPIEs on evaporation of ethanol. Although the quantitative agreement with the experimental values is not there yet, our study sets some recommendations regarding the models and their theoretical treatment and offers a new approach for selecting frequencies used to compute isotope effects. It also highlights the necessity of having models for organic liquids such as alcohols to be constructed in such a way that not only hydroxyl part of the isotopically substituted molecule is surrounded by the proper number of other molecules but also the sufficient number of interactions at methyl and methylene groups is provided, if the prediction of carbon isotope effects is of interest.  Figure S1 in the Supporting Information). Geometries of 377 molecules contained within the selected sphere were initially optimized to a gradient of 0.3 kcal·mol −1 ·Å −1 using Amber99 37 force field constraining positions of the remaining molecules. Subsequently, the molecules outside the sphere were removed and a smaller sphere of 15 Å radius, containing 174 molecules was constructed (see the right panel of Figure S1 in the Supporting Information), and the optimization process was repeated to a gradient of 0.1 kcal·mol −1 ·Å −1 . Removal of molecules from the outer sphere and reoptimization did not cause noticeable puckering of the sphere, which was thus used in the preparation of hybrid (ONIOM) models by selecting atoms within 10 Å radius (91 molecule model) and 12.5 Å radius (138 molecule model). These models have been divided into a QM part (8 and 14 molecules, respectively) and a MM part using GaussView 38 and submitted to energy minimization to local minima using the ONIOM protocol 39 implemented in Gaussian. 40 A number of combinations of theory levels and protocols have been tested. These include two-layer QM/MM and QM/QM, and three layer QM/QM/MM protocols (as exemplified in Figure  1) with mechanical or electronic embedding. B3LYP, 41 M11, 42 B97D3, 43 and ωB97xD 44,45 functionals were used with def2-SVPP, def2-TZVP, 46 6-31+G(d,p), 47 6-311++G(d,p), 48 and aug-cc-pVTZ 49 basis sets. Polarized continuum models in the default IEF 50 form and Cosmo 51 or SMD 52 variants were used. The MM part of these models were modeled using available force fields: Amber and Dreiding as implemented in the Gaussian program.
So-called cluster models were also constructed and tested in this work. 4-, 8-, and 17-molecule clusters were prepared. Smaller cluster models were created from the QM/MM models by selecting four or eight molecules within the center of the two models described above. The 17-molecule cluster model was built based on the larger QM/MM model by selecting central molecules and all molecules within 5 Å radius from each carbon and oxygen atom of the central ethanol molecule. Cluster models were entirely treated by the same respective QM level of theory; the ωB97xD functional combined with the aug-cc-pVTZ or def2-family basis sets. In some cases, the continuum solvent model (IEFPCM) was additionally applied.
The gas phase was generally modeled by considering a single molecule optimized at the corresponding theory level. Only in two cases, a comparison of clusters (of 8 and 17 molecules) which geometries were optimized using the PCM model were used as models of the gas phase and compared with clusters of the same size but optimized without the continuum solvent model.
All cluster models including their liquid and gas phase counterparts were optimized to a local minimum using the very tight convergence criteria as well as ultrafine integrals.
Single structure representative was subsequently used for frequency and isotope effects calculations. In one of our previous studies, 24 we have explored this issue in more detail, and we observed that if the model and the theory level provide the prediction that is in agreement with the experimental value of an isotope effect, then the number of structures used for its computation does not matter that much.
Isotope Effect Calculations. VPIEs were calculated from the partition functions using the equation, introduced by Bigeleisen 53 where, u i = hν i /(k B T); ν i are isotopic frequencies of light (L) and heavy (H) species in liquid (l) and gaseous (g) phases, h and k B are Planck's and Boltzmann's constants, respectively. The simplification is valid under the Born-Oppenheimer approximation assuming pure harmonic frequencies. No scaling of frequencies was used. Two sets of frequencies were used for liquid ethanol models. In the first one, all frequencies obtained from the entire (supermolecule) respective (hybrid or cluster) system were used. All or either 3N or 3N − 6 modes were included in the calculation, where N denotes the number of atoms in a system. In the second set, only internal molecular modes of each molecule constituting the model were used, whereas any collective modes (i.e. ones associated with other molecules in a system) were omitted. Thus, a new model for deriving isotope effect based on normal modes of a system was implemented. We called it a "frequency-cut-off" model. Calculations were performed using our Isoeff program. 55