Can Path Integral Molecular Dynamics Make a Good Approximation

Subscriber access provided by Technical University of Munich University Library

B: Liquids, Chemical and Dynamical Processes in Solution, Spectroscopy in Solution

Can PIMD Make a Good Approximation for Vapor Pressure Isotope Effects Prediction for Organic Solvents? A Comparison to ONIOM QM/MM and QM Cluster Calculation Luis Vasquez, Michal Rostkowski, Faina Gelman, and Agnieszka Dybala-Defratyka J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b03444 • Publication Date (Web): 01 Jul 2018 Downloaded from http://pubs.acs.org on July 1, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Can PIMD Make a Good Approximation for Vapor Pressure Isotope Effects Prediction for Organic Solvents? A Comparison to ONIOM QM/MM and QM Cluster Calculation

Luis Vasquez,a Michal Rostkowski,a Faina Gelman,b and Agnieszka Dybala-Defratykaa,*

a

Institute of Applied Radiation Chemistry, Faculty of Chemistry, Lodz University of Technology, Żeromskiego 116, 90-924 Łódź, Poland, e-mail: [email protected]; Tel: +48 42 631 3198 b Geological Survey of Israel, Malkhei Israel St. 30, 95501, Jerusalem, Israel

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract Isotopic fractionation of volatile organic compounds (VOCs), that are under strict measures of control due to their potential harm to the environment and humans, has an important ecological aspect, as the isotopic composition of compounds may depend on the conditions in which such compounds are distributed in nature. Therefore, detailed knowledge on isotopic fractionation, not only experimental but also based on theoretical models, is crucial in order to follow conditions and pathways within which these contaminants are spread throughout the ecosystems. In this work we present carbon and, for the first time, bromine vapor pressure isotope effect (VPIE) on evaporation process from pure-phase systems - dibromomethane and bromobenzene, the representatives of aliphatic and aromatic brominated VOCs. We combine isotope effects measurements with their theoretical prediction using three computational techniques, namely PIMD, QM cluster, and hybrid ONIOM models. While evaporation of both compounds resulted in normal bromine VPIEs, the difference in direction of carbon isotopic fractionation is observed for the aliphatic and aromatic compounds, where VPIEs are inverse and normal, respectively. Even though theoretical models tested here turned out to be insufficient for quantitative agreement with the experimental values, cluster electronic structure calculations, as well as two-layer ONIOM computations, provided better reproduction of experimental trends.

2 ACS Paragon Plus Environment

Page 2 of 36

Page 3 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Volatile organic compounds (VOCs), defined in the European Union as organic compounds with a starting boiling point less than or equal to 250 °C at 1.013 MPa, are under strict measures of control due to their potential harm to human life1–3 as well as environmental concern regarding pollution4–6 and degradation.7 Among known VOCs, their brominated derivatives constitute a specific group of compounds as they have been identified as compounds which concentration directly determines treated water quality,8 toxicity,9,10 but most important, during the last three decades organobromine compounds have been pointed out to be likely cancerogenic11 and to possess large pollutant/harmful effects.12–15 Evaporation is one of the important pathways of VOCs attenuation driven by the vapor pressure that a compound possesses. One of the possible ways of quantifying this process is to study the isotopic composition of the compound phases, which may be enriched or depleted in either lighter or heavier isotopologue of the studied element as compared to the initial composition. When enrichment in the heavier isotopologue occurs in liquid phase isotopic composition variance is then classified as normal vapor pressure isotope effect (VPIE); in contrast, when the evaporation rate of the lighter isotope is slower then we deal with inverse VPIE. According to isotope effect theory, VPIEs are mainly governed by intermolecular forces of a different kind (hydrogen and halogen bonds, van der Waals forces, and any other non-covalent interactions between the species present in a system). Knowing the magnitudes and the direction of VPIE one can try to characterize the structure and nature of phases in more detail, in particular, a liquid one, and provide information on the intermolecular interactions as well as changes that the studied compound may undergo at the transition from one phase to the other. Such knowledge is particularly important for constructing more accurate models for studying organic species dissolved in water. However, before we are able to do it we need to make sure that we can predict isotope effect for pure compounds, capture any differences observed experimentally, and analyze them using information coming from the computational study and concomitant measurements. One has to, however, bear in mind that existing experimental data on various organic compounds either pure-phase or dissolved in water is generated using different experimental set-ups that may affect isotopic composition values,16–18 therefore direct comparison among the studies may not be so straightforward and additional analyses should be designed and performed.



Page 4 of 36

The aim of the presented study is three-fold. Firstly, testing some available computational approaches to predict VPIEs and compare them with the experimentally determined values of IEs. For this purpose, we selected two compounds; dibromomethane (DBM) and bromobenzene (BB) as representatives of brominated VOCs, clear examples of aliphatic and aromatic one, respectively. Secondly, with a use of carbon and bromine isotopic analysis, we hoped to explore possible differences in the direction of measured and predicted isotope effects. The isotope exchange reactions for DBM and BB and two elements: carbon and bromine can be written as follows: CH2Br2 (l) + 12CH2Br2 (g) ⇋ 12CH2Br2 (l) + 13CH2Br2 (g)

(1)

CH2Br81Br (l) + CH2Br2 (g) ⇋ CH2Br2 (l) + CH2Br81Br (g)

(2)

CH281Br2 (l) + CH2Br2 (g) ⇋ CH2Br2 (l) + CH281Br2 (g)

(3)

13

C6H5Br (l) + 12C6H5Br (g) ⇋ 12C6H5Br (l) + 13C6H5Br (g)

13

C6H581Br (l) + C6H5Br (g) ⇋ C6H5Br (l) + C6H581Br (g)

(4) (5)

where the subscripts l and g denote the liquid and vapour phase, respectively. Finally, with an insight coming from a computational study we made an attempt to investigate possible origins of observed differences. To this end, the path integral formalism of quantum chemistry, as well as ONIOM scheme,19 and QM cluster calculations, were used. The imaginary time path integral formalism20 of quantum chemistry has proven to be an effective method to quantize selected isotope effects on various processes previously.21–26 However, it has not been employed before to explore vapor pressure isotope effects in pure phase, except for water.27–30

Methodology I. Volatilization experiment Pure phase volatilization of dibromomethane (99%, Aldrich) and bromobenzene (99%, Aldrich) was performed in a fume hood at 22 °C. One mL of the brominated organic compound was placed into 2 mL open glass vial and allowed to evaporate until 5%. The samples of 5 µL were periodically taken from the vial. The vials were weighted before and after each sampling to determine the amount that had been evaporated and sampled from the vial. Carbon and bromine isotope analysis of the evaporated samples were performed by gas


Page 5 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


chromatography isotopic ratio mass spectrometry (GC-IRMS) and gas chromatography multiple ion collectors inductively coupled plasma mass spectrometry (GC-MC-ICPMS). II. Carbon and bromine isotope analysis For carbon and bromine isotope analysis the samples of the evaporated brominated organic compounds were dissolved in methanol. Operating parameters for GC-IRMS and GC-MCICPMS are listed in Table 1. Table 1. Operating parameters for GC-IRMS and GC-MC-ICPMS analysis. GC-IRMS

GC-MC-ICPMS

GC parameters: DB-5 column (30m, 0.25mm, 0.25µm); injector temperature 200°C, split 1/10; oven temperature: 40°C for 4 min, 8°C/min to 200°C; He carrier flow 1 ml/min. IRMS parameters:

MC-ICPMS parameters:

oxidation oven 940°C; reduction oven 650°C

RF power 120W; coolant flow 12 L/min; auxiliary flow 1.45 L/min; nebulizer gas flow 0.68ml/min; sweep gas (argon) 5.45 L/min; spray chamber temperature 75°C; m/Δm = 300

To determine carbon or bromine isotope enrichment factors (ε) along evaporation of dibromomethane and bromobenzene , a modified Rayleigh equation was applied: ோ

ln ቀோ೟ ቁ = ߝ݈n‫ܨ‬

(6)

೚

where R0 and Rt are isotopic ratios at time zero and time t of evaporation, respectively, and F is the remaining fraction of the pure compound.

Computational Part I. Geometry Optimization The geometries of BB and DBM were optimized using the Gaussian software package revision 09 E.01.31 Tight optimization convergence criteria were used along with three levels of theory: DFT with the B3LYP32–35 and the M06-2X36 functionals, the second order Møller– Plesset perturbation theory (MP2),37,38 and coupled cluster with double excitations (CCSD).39,40 In each of these calculations, the 6-311+G(d,p) basis set was used41,42 as its 5 ACS Paragon Plus Environment


relative accuracy and convergence time are an advantage for calculations using relatively high levels of theory such as MP2 and CCSD. Moreover, the ultrafine grid was used for computing the two-electron integrals as well as their derivatives. Once the geometry optimization of the molecules converged, the electrostatic partial charges (ESP) for given geometries were calculated using the Merz-Singh-Kollman scheme.43,44 ESP charges were obtained using the RESP ESP charge Derive Server Development scheme using the default settings,45 while general Amber force field (GAFF) parameters46 were taken from the database provided by the R.E.D Server Developments,47 and/or geometric optimizations if needed. II. Molecular Mechanics All molecular dynamics (MD) simulations were carried out in the AMBER simulation package, version 1448 by the SANDER module (Simulated Annealing with NMR-Derived Energy Restraints). Cubic-periodic solvent boxes of 40 angstroms (Å) were prepared based on optimized structures with four aforementioned methods for each of the compounds, respective GAFF parameters and ESPs in the LEaP module, in AMBER14. Then, the energies of the resulted boxes were minimized using the steepest descent method (25000 cycles) followed by the conjugate gradient method (55000 cycles). Temperature (298.15 K) was regulated by the Langevin dynamics thermostat49,50, while the pressure was kept constant at 0.1 MPa. Given the large size of the periodic boxes, a small time step 0.5 fs was used, the simulation time of the heating and equilibration was 1.1 ns, the production phase was run by the same amount of time. The long range cut-off for non-bonded interactions was 10 Å, however, accuracy to the calculations was added by decreasing the tolerance value of the Particle mesh Ewald (PME) to 10-6. Periodic boundary conditions were used throughout the simulations. III. Quantum Mechanics/Molecular Mechanics Molecular Dynamics In order to explore the systems behavior upon division into the two regions; one treated quantum mechanically and the other classically, we chose the PM6 Hamiltonian51 to represent a single, central molecule of each system. The rest of the system was treated using the aforementioned force field parameters. The protocol used during the QM/MM MD simulations was the same as the one used during the MD simulations except for the long range cut-off; here a switching function was set to 11 and 9 Å for an outer and inner radius, respectively. 6 ACS Paragon Plus Environment

Page 6 of 36

Page 7 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


IV. Path Integral Molecular Dynamics Path Integral Molecular Dynamics (PIMD) is one of the available methods to determine nuclear quantum effects and its theoretical and mathematical bases are known and available 20–23,25

.

Normal Mode Path-Integral Molecular Dynamics (NMPIMD)52 simulations were performed for each of the four periodic boxes. The equilibrium properties of the system (at 298.15 K, 0.1 MPa) were computed. Given that the size of the systems as for PIMD simulation are relatively large (1757 atoms for BB and 503 atoms in the case of DBM) and that the conservation of energy must be ensured along the simulation, a conservative (0.02 fs), timestep was selected which conveniently is one order of magnitude smaller than one used during the MM and QM/MM MD simulations. Each NMPIMD simulation had a duration of 0.8 ps, long enough to ensure the equilibration of the system by the NMPIMD, but also to observe possible variations/distortions along the simulation. The default thermostat was NoséHoover, with 8 chains, whereas the relevant quantities were calculated by the virial-like estimator.53 It is important to note that the PIMD simulation requires the system to be in equilibrium (system must reproduce the experimental density of the solvent) beforehand, thus, owing to the lack of MD simulations to include polarization and taking into account the inclusion of dipole-moment in QM calculations, the equilibrated systems were taken from the QM/MM MD simulations. The number of imaginary time slices in the path integral formalism chosen for DBM were 8, 16, 24, 32, and 64, whereas for BB 8, 12, 16, 20, 28, 32, 36, 40, 48, 64, and 128. Those numbers were selected using the most common settings from the literature and the size of the molecule in the study. We used the NMPIMD52 implementation in AMBER14. For the purpose of predicting vapor pressure isotope effects (VPIEs), path integral quantum transition state theory54 was used. Within this theory, a thermodynamic integration with respect to mass is computed (potential energy surface does not change, while mass does). During thermal integration calculations interpolation between the masses (λ) of light (0) - 12C and 79Br and heavy (1) - 13C and 81Br 55 isotopologues using 9 in-between values for DBM (0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 and 1) and 7 for BB (0, 0.01, 0.02, 0.1, 0.2, 0.4 and 1) were used.

13

C and

81

Br isotope effects were calculated for each position within the respective

molecule, and subsequently, the isotope effect was estimated by averaging in the case of the carbon VPIE for BB and bromine VPIE for DBM.



Page 8 of 36

V. MO:MM ONIOM and MO Clusters Calculations The estimation of VPIEs using two-layer MO:MM ONIOM computational scheme was also performed. For this purpose, nine structures from each MM MD simulation were chosen based on the lowest total energy criterion. Obtained snapshots where then decreased in size by extracting only molecules that had at least two atoms within a radius of 6 and 10 Å from a central molecule using VMD.56 These spherical systems were used as input geometries for ONIOM calculations performed in Gaussian. The QM region being one central molecule, either BB or DBM, was described by B3LYP, M06-2X, MP2, and CCSD model chemistries. To describe outer part of a system with the molecular mechanics, the same sets of atomic charges as well as the AMBER force field parameters were used as for the classical MD simulations described above. Electronic embedding technique was used to describe the electrostatic interactions between the QM and MM regions. Model for vapor phase consisted of only one molecule either of BB or DBM treated by the four levels of theory in the gas phase. Additionally, (DBM)n and (BB)n clusters, where n=18, 12, 10, 8, 6, 4, and 3 were constructed. They were optimized using the PM7 Hamiltonian57 implemented in the Gaussian 16 package.58 Each cluster was optimized to a minimum using tight convergence and optimization criteria. Subsequently, normal mode analysis was performed to confirm that the optimized structures were minima and to generate Hessian matrices necessary to compute isotope effects. Position specific carbon (k12/k13=C VPIE), bromine (k79/k81=Br VPIE), carbon average effect (bulk C VPIE), and bromine average effect (bulk Br VPIE) were calculated using the Bigeleisen equation59 and the ISOEFF07 program60 at 300 K. All predicted isotope effects throughout the paper are expressed in a form of enrichment factors (ε) based on the equation (7) relating IE with compound average ε. ଵ

ଵ

ߝ = ∙ ቀ − 1ቁ ௡ ூா

(7)

where n is the number of atoms of an element in a molecule. ε are usually expressed in numbers per mill (‰).


Page 9 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Results and Discussion Experimental Part Evaporation of pure-phase liquids from the open vials was accompanied by changes in carbon and bromine isotope composition of the residual fraction of the studied brominated organic compounds. As it can be seen in Figure 1, for both studied brominated organic compounds enrichment of the residual non-vaporized fraction by heavier bromine isotope (81Br) was observed and resulted in close values of bromine isotope enrichment factor (-0.8 and -0.9‰ for dibromomethane and bromobenzene, respectively). However, different trends in carbon isotope fractionation were observed for these compounds. Whereas evaporation of bromobenzene resulted in a little normal carbon isotope effect (ƐC= -0.4‰), an inverse carbon isotope

fractionation

with

ƐC=+0.6‰

was

observed

along

the

evaporation

of

dibromomethane.

Figure 1. Rayleigh plots for carbon and bromine isotope fractionation along bromobenzene (left) and dibromomethane (right) evaporation.

Computational Part I. Initial Systems Characteristics The geometrical parameters from ab initio optimizations for single DBM and BB molecules obtained with the B3LYP, M06-2X, MP2, and CCSD methods are provided in the Supporting Information (SI, Tables S1 and S2). All models resulted in almost the same geometry. Generated partial charges, as well as computed dipole moments for the two studied 9 ACS Paragon Plus Environment


molecules, are shown in Tables S3 and S4. Although none of the methods resulted in a perfect match with the experimentally determined values of dipole moments, all of them are similarly overestimated, with B3LYP DFT functional providing the closest agreement for both molecules among the four tested methods (1.79 vs experimental 1.43 and 1.97 vs experimental 1.70 D for DBM and BB, respectively). Next, the geometries resulting from electronic structure calculations were compared to the existing Amber force field parameters in order to identify any relevant differences. For example, in the worst of the cases, the HC – CT – HC angle in DBM was found to be overestimated by quantum mechanical calculations by 3 degrees (B3LYP), whereas in the case of the CA– CA – CA* angle (CA* carbon atom bonded to Br) in BB it was smaller than usual CA – CA – CA angle (by 0.7 degrees on average). Nevertheless, such differences are insignificant and none of the force field parameters except atomic charges was modified in either of the molecules prior to the subsequent simulations. All force field parameters used in the simulations are presented in Tables S5 and S6 in SI. II. Molecular Dynamics Simulations Four molecular dynamics simulations were carried out starting from the models constructed based on the compounds optimized using the B3LYP, M06-2X, MP2, and CCSD methods. With this in mind, for convenience, these notations are used exclusively for naming the models within which the single (base) molecule of either DBM or BB was initially optimized using either of these methods. In respective solvent boxes preparation, a special focus was put on solvent densities as it is one of the most important issues to take care of before performing the PIMD simulations. The systems constructed based on the CCSD and MP2-optimized geometries of DBM reached the stability earlier, at about 200 ps, whereas those based on the two DFT functionals around 100 ps later (Figure S1). In the case of BB, relatively all models equilibrated at about by 200 ps (Figure S2). Next, for the purpose of having the average density from each simulation under equal conditions, a statistical analysis was performed from 0.35 to 1.3 ns of simulation time. The former value was selected as the starting point since after this particular time all simulations were considered equilibrated. As the optimized geometries from the quantum mechanical calculations relatively did not differ from each other, consequently, any observed differences are mostly related to the computed ESP (Tables S3-S6).


Page 10 of 36

Page 11 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Statistical analysis of density showed that among all four DBM models M06-2X reproduced the experimental density of 2.477 g⋅cm-3

61

with great accuracy (2.478 g⋅cm-3), B3LYP

provided a number slightly underestimated (2.471 g⋅cm-3), and both MP2 and CCSD resulted in values slightly overestimated (2.484 and 2.485 g⋅cm-3, respectively). BB, on the other hand, has an experimental density equal to 1.4950 g⋅cm-3

55

, calculated density is

-3

underestimated by 8 % on average, B3LYP – 1.370 g⋅cm , CCSD – 1.380 g⋅cm-3, M06-2X – 1.376 g⋅cm-3, and MP2 – 1.376 g⋅cm-3. As the two systems under study are liquids at 25 ºC, it is essentially important to analyze the structural distribution of DBM molecules in the simulated box. For this purpose, an approach employing pair correlation functions was applied and radial distribution function (RDF) was calculated through particle number concentration at a radial distance from a determined (selected) particle (atom/residue number) center.62 As a center of the DBM system, the carbon atom (CT) was selected and its spatial correlation between carbon atoms belonging to the surrounding DBM molecules was calculated within a distance of 15 Å. By analogy the procedure was repeated for the Br – HC, and CT – Br pairs of atoms. RDF calculations were performed only for the data from the period within which the system had reached equilibrium. Figure 2 presents the RDFs for DBM, while Table 1 shows the data for BB. Figure reveals that in the case of DBM models share alike characteristics, i.e. the same shape of the RDF is preserved. RDF for bromine atoms (Figure 2a) around the central carbon atom exhibits 3 peaks. There are at least two main peaks for the carbon atoms from the neighboring molecules around the central carbon atom of each DBM molecule (Figure 2b) one of the peaks is found between 5 and 6 Å, and from 11 Å on, RDF converges to the values close to unity, which is also a common feature for all computed RDFs for the DBM system. In the case of the hydrogen atoms surrounding the bromine atom (Figure 2c) an interaction around 3 Å was observed.



Figure 2. Radial distribution function of a) the bromine atoms (Br) around the carbon atom (CT), b) carbon atoms (CT) around the carbon atom (CT), and c) the hydrogen atoms (HC) around the bromine atom (Br) in the DBM box.

In contrast to DBM, BB has been previously studied experimentally by others using the Xray, neutron diffraction, and extended X-ray absorption fine structure (EXAFS).63 Their results (Table 2, the second row) indicate that there are four shells around the bromine atom in BB, given by the different Br-CA# configurations, where # represents the number assigned to each of the carbon atoms that form the BB aromatic ring (CA1 to CA6). Within the third shell Bertagnolli et al. have found another distance of 3.70 Å which was interpreted as a combined effect of the interaction of the electron system of the aromatic ring and the bromine atom symmetrically positioned above the ring plane and van der Waals interaction which occur at longer distances. In order to compare the structures of our BB models to that experimental data, we calculated RDFs of the carbon atoms (CA1 to CA6) of the surrounding residues within the 15 Å range around the bromine atom. The four models agree unanimously, around the bromine atom four shells are being formed, specifically at 1.85 Å, 12 ACS Paragon Plus Environment

Page 12 of 36

Page 13 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.85 Å, 4.15 Å, and 4.65 Å. Beyond 11 Å the RDF values approach the unity. However, the interactions within 3.7 Å were not observed in the resulted configurations of the models. Computed Br−Br RDF did not reproduce Bertagnolli et al. experimental values showing that there are two distinguishable peaks (3.79 and 4.82 Å). In contrast, we found only one bromine RDF peak at 4.15, 4.25, 4.15, and 4.45 Å for B3LYP, CCSD, M06-2X, and MP2, respectively. These predicted values are close to the geometrical average (4.3 Å) of both shells recorded experimentally. Such geometrical flaw may be responsible for density issues, given that instead of having two bromine atoms in the surroundings (3.79 and 4.82 Å) only one atom (4.45 Å) is found within the interaction ranges (Figure S3). Also based on this comparison between the experimental and theoretical structures it is obvious that some of the configurations proposed by Bertagnolli et al. were not abundant in the models resulted from the performed MD simulations.

Table 2. Radial distribution function comparison between experimental and computed results, bromine – carbon density (Br-CA1 to Br-CA6) for BB molecule.

a

Parameter

Shell 1

Shell 2

Atom Distance R (Å)a B3LYP (Å) CCSD (Å) M06-2X (Å) MP2 (Å)

C 1.88 1.85 1.85 1.85 1.85

C 2.84 2.85 2.85 2.85 2.85

Shell 3 C 3.70

C 4.02 4.15 4.15 4.15 4.15

Shell 4 C 4.74 4.65 4.65 4.65 4.65

Liquid BB local order obtained from X-ray, neutron diffraction and extended X-ray absorption fine structure

(EXAFS)63

III. Quantum Mechanics / Molecular Mechanics Molecular Dynamics Simulations Subsequently, we used a QM/MM approach to allow for polarization of at least a central molecule of the system, represented by the semiempirical Hamiltonian, by surrounding remaining molecules treated by the same force field as previously. The QM/MM MD density results (Figure S4 and S5) for both DBM and BB relatively do not differ from that the ones obtained using classical MD. The slight difference observed for DBM are in the distances at which crucial interactions take place (Figure S6). Even though, within the QM/MM MD the atomic charges of the QM part have changed, and thus its interactions with surrounding 13 ACS Paragon Plus Environment


Page 14 of 36

molecules having fixed charges, comparison of charge between systems, models and even frames was not possible, owing to the fact that the Mulliken charges in the output are not comparable to the charges derived from the electrostatic potential,64 either does not meet requisites for ESP calculations.45 IV. Path Integral Molecular Dynamics Simulations It should be noted that PIMD is based on the principle of imaginary time slices and convergence of an estimator used to compute the derivative of the free energy efficiently, therefore one needs to test various numbers of time slices to obtain converged results. However, as convergence depends on the initial conditions and the system's nature, the number of imaginary time slices differs from system to system, and it should be obtained empirically.48 In the presented PIMD simulations a so-called virial estimator (VE)65,66 was applied and the relative differences between its values obtained with chosen time slices in increasing order were used as the criterion of convergence. Limitations on the computational resources can be a serious drawback while performing the PIMD calculations, notably in the case of isotopic fractionation by thermodynamic integration, as the interpolating mass parameter (λ) requires as many additional calculations as the total number of interpolation points, which are 9 and 7 in the case of DBM and BB, respectively. Despite the fact that all DBM models were initially equilibrated using PIMD in order to monitor their performance (8, 16, 24, 32, and 64 time slices), only two models along with the selected, based on test, number of time slices could have been selected for posterior thermodynamic integration and the final VPIE estimation. B3LYP (8, 32 and 64) and CCSD (16, 32 and 64 slices) models seemed to be suitable for the VPIE estimation. Nonetheless, virial estimator convergence as a function of imaginary time slices during thermodynamic integration is provided for the two selected models (Figures S7 and S8), where, for instance, substantial reduction of converged energy relative to 64 slices in

13

C (7.2 and 0.5% for

B3LYP, 8 and 32 slices respectively) is seen. Table 3 contains average 81Br and 13C VPIEs obtained for evaporation of DBM using B3LYP (8, 32, and 64 slices) and CCSD (16, 32, and 64 slices) models. In the case of bromine, two calculations were performed, one in which only one of the bromine atoms was substituted by heavier isotopologue (81Br), and the other in which both bromine atoms were substituted (81Br 81Br). For instance, the mean VPIE value for bromine (81Br – B3LYP 8 slices) is 3.5‰, while for carbon (13C) is -0.3‰, which gives the opposite trend to the experimental 14 ACS Paragon Plus Environment

Page 15 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


observation. Overall, bromine isotope effects were overestimated in both models, yet, their magnitudes are becoming moderately less inverse with a larger number of imaginary time slices. For instance, the B3LYP model passes from 3.5 with 8 slices via 1.9 with 32 to -0.1‰ with 64 slices. On the other hand, only in three (B3LYP 32 and CCSD 16 and 64) out of the four sets the observed trend based on the measurements was also reproduced by the predicted carbon VPIEs.

Table 3. Carbon and bromine VPIEs measured and predicted for evaporation of DBM at 298.15 K and 0.1 MPa. Method

B3LYP

CCSD Experimental a

81

8 32 64 16 32 64

Br 3.5 1.9 -0.1 2.5 1.3 1.1 -0.8

Isotope enrichment factor (ε,‰) 81 Br 81Br 2.5 1.9 2.3 1.3

13

C -0.3 3.2 -1.3 9. -0.5a 2.8 0.6

Value obtained using 8 mass interpolation points (λ).

In contrast, BB during both MD and QM/MM MD simulations required more effort for achieving equilibrium and reproducing properties such as density (note time differences in the Figure S2 and S7). For this reason, a wider range of imaginary slices was studied. Figure 3 and 4 contain mean virial-estimator convergence in a function of imaginary slice number for specific carbon position in BB and for bromine, respectively. Although the differences between the VE values are fractions of kcal⋅mol-1, remarkable reduction on values was observed with 8 and 28 slices as compared to 128 slices (the relative energy difference are 8.8 and 1.9% respectively, thus showing that likely calculated VPIE with less than 28 slices may not be sufficiently accurate prediction of

13

C VPIE). A rough estimation, although it

depends on a property of interest, for the required amount of slices in order to converge a PIMD simulation can be estimated by multiplying 2 or 3 times the term 2ߚℏ߱௠௔௫ ,67 where β represents the inverse temperature and ωmax gives the maximum frequency of interest, which gives ܲ ≈ 16 − 24 in the case of BB if ωmax is around 3400 cm-1.68,69 Position CA1 refers to a carbon atom forming a bond with bromine. Comparatively speaking for most of the slices it possesses the highest mean converged energy, while the adjoining CA2 and CA6 carbon 15 ACS Paragon Plus Environment


atoms, in contrast, keep the lowest converged energies (Figure 3). Nuclear quantum effects are supposed not to be large on bromine, mainly owing to its relatively large mass as compared to carbon, effectively conserved energy values are smaller than those for carbon. Surprisingly, conserved energy gradually decreases to 32 imaginary slices, from that point onwards, it only increases surpassing its initial value.

Figure 3. Mean estimator convergence as a function of imaginary slices number for specific carbon position within the BB molecule obtained using PIMD and CCSD model chemistry for the molecule initial optimization and ESP derivation.


Page 16 of 36

Page 17 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Mean estimator convergence as a function of imaginary slices number for bromine within the BB molecule obtained using PIMD and CCSD model chemistry for the molecule initial optimization and ESP derivation scheme.

It is important to realize that although the inset plot in Figure 4 might indicate that 32 and 128 imaginary slices provide better energy convergence along the PIMD simulations as the mean energy resulted from them is the lowest, it is not necessarily true as it should be taken into account that the plot scale covers differences of hundredths of cal⋅mol-1 at most, which is relatively low having in mind that the whole system has been quantized (more than 1500 molecules) during the simulation. Although it might be expected that carbon VPIE is more pronounced in comparison to bromine VPIE there are six carbon atoms within a molecule and the bulk carbon isotopic fractionation gets diluted. Figure 5 presents average carbon VPIEs for the BB evaporation as a function of number of imaginary time slices. The first set of slices (8) shows that except for the CA3 carbon position (Table S7) all computed isotope effects are inverse, which does not agree with the experimental trend (-0.4‰) and definitively indicates that 8 imaginary slices are not enough for predicting quantum nuclear effects on the BB evaporation. With a gradual increase of imaginary slices position specific carbon isotope effects equally tend to change toward normal values in general, for instance, with 36 slices, two carbon positions (CA4 and CA5) exhibit VPIE values smaller than 0 (-1.2 and -1.2‰, Table S7). With exception of 40 (17 ACS Paragon Plus Environment


0.2‰) and 64 (-3.6‰) imaginary slices, none other set of imaginary slices resulted in providing bulk carbon VPIE close to the experimental value even direction wise. Bromine VPIEs resulted from the PIMD simulations are given in Table 4. With the exception of 64 slices, all sets provided inverse isotope effect.

Figure 5. Average carbon VPIEs predicted for evaporation of BB as a function of imaginary slices at 298.15 K and 0.1 MPa..

Previously, Paesani et al.70 have revealed that in certain cases of using PIMD, parameterizations had to be done, but as force field parameters seemed to allow to predict liquid properties of studied solvent quite accurately, as demonstrated with a use of density and spacial molecules distribution in generated boxes, no further parametrization was done. However, one might consider it in any future study in order to tune up the convergence70 and thus decrease the required amount of slices for converging a determined property.67 In that case, it might be expected that given that convergence is stable, less interpolation points of mass will be needed for obtaining reliable energy differences for the system.

Table 4. 81Br isotope enrichment factors (ε) predicted for evaporation of BB at 298.15 K and 0.1 MPa using the PIMD and the CCSD model of the compound. № imaginary slices 8 12 16 20

Isotope enrichment factor (ε, ‰)

0.4 0.7 0.8 1.0 18

ACS Paragon Plus Environment

Page 18 of 36

Page 19 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1.1 1.6 2.6 1.3 1.5 0.0 0.8

28 32 36 40 48 64 128

V. VPIEs Prediction Based on the ONIOM and QM Cluster Calculations Conformations resulted from classical MD simulations were subsequently used for the ONIOM calculations using either B3LYP, M06-2X, MP2, and CCSD combined with the 6311+G(d,p) basis set for the higher (MO) layer which comprised the central either the BB or DBM molecule and the GAFF parameters for the remaining either solvent molecules (MM, lower layer). Two sizes of models were constructed, one comprising all molecules within 6 Å radius from the central molecule and the other comprising all molecules within 10 Å radius. Structure-specific isotope effects obtained for BB and DBM using the ONIOM scheme were collected and are presented in Tables S8-S11. Respective averaged values are given in Tables 5 and 6. In the case of BB, the qualitative trends (normal fractionation for carbon and bromine) were reproduced. However, the magnitudes of slightly overestimated, whereas

81

13

C enrichment factors seem to be

Br enrichment factors underestimated as compared to the

measured values. No significant differences were found neither among the models nor the model sizes. In contrast, for DBM only one model resulted in the matched prediction of carbon isotope effect with the experimental values (M06-2X). The rest of the models predicted the opposite trend. In the case of bromine, on the other hand, the same direction of enrichment factors was predicted but the values are again highly underestimated. We also subjected single structures of models of varying sizes to the pure QM calculations using the PM7 method - quite recently developed with a goal of being good and reasonably accurate semiempirical method for predicting non-covalent interactions such as hydrogen and halogen bonds as well as dispersion. It has been also shown to be superior to other well-known semiempirics such as PM6 or PM3.71 The isotopic fractionation in the above analysis was predicted using either only the central molecule of each cluster or by calculating isotope effects for each of the molecule comprising a cluster and then averaging them over all molecules present in it. In the case of smaller clusters (3-6 molecules), the choice of central 19 ACS Paragon Plus Environment


Page 20 of 36

molecules was rather arbitrary. The resulted magnitudes of isotopic fractionation factors are shown in Figures 6-9. In the case of carbon fractionation for BB volatilization (Figure 6) the values of epsilon obtained by averaging over all molecules in the system seem to be a better choice as they converged around -0.4‰ for clusters consisted of 8, 10, 12, and 18 molecules, and are in good agreement with the experiment. The values resulted from the prediction of isotope fractionation only for the central molecule in those clusters are mostly overestimated (~ -0.9 ‰). The results obtained for smaller clusters are either underestimated or provided almost a perfect reproduction but it could be purely accidental.

13

C enrichment factors for DBM

evaporation (Figure 8) show a similar trend for values averaged over all molecules. These values are slightly underestimated with one exception for the cluster containing 12 molecules, that happened to be a perfect match with the measured value. Also quite similar to the BB clusters containing 8, 10, and 12 are the ones that provided overestimated fractionation but to a lesser extent than in the case of BB. Prediction of 81Br enrichment factors (Figures 7 and 9) always led to underestimation. Values resulted from averaging over all molecules in a cluster turned out to be in a worse agreement with the experiment than the magnitudes obtained for one central molecule, in particular, for larger clusters (10, 12, 18). That is, a central one is mostly buried, thus is affected by the highest number of interactions, that elevates the value of isotope effect in the expected direction. Table 5.

13

C isotopic enrichment factors (ε, ‰) for evaporation of BB and DBM predicted using

various ONIOM=MO:AMBER approaches for two different sizes of the systems. MO

a

Model_size

B3LYP

M06-2X

MP2

CCSD

BB_6Å

-0.7 ± 0.2

-1.0 ± 0.4

-0.8 ± 0.4

n/ca

BB_10Å

-0.9 ± 0.2

-1.1 ± 0.4

-0.6 ± 0.7

n/c

DBM_6Å

-0.6 ± 0.5

0.4 ± 0.3

-0.9 ± 0.2

n/c

DBM_10Å

-1.0 ± 0.4

-0.1 ± 0.3

-1.0 ± 0.3

-0.9 ± 0.4

n/c - not calculated

Table 6.

81

Br isotopic enrichment factors (ε, ‰) for evaporation of BB and DBM predicted using

various ONIOM=MO:AMBER approaches for two different sizes of the systems.


Page 21 of 36

MO

a

Model_size

B3LYP

M06-2X

MP2

CCSD

BB_6Å

-0.1 ± 0.1

-0.2 ± 0.1

-0.1 ± 0.0

n/ca

BB_10Å

-0.1 ± 0.0

-0.2 ± 0.0

-0.1 ± 0.1

n/c

DBM_6Å

-0.1 ± 0.1

-0.1 ± 0.0

-0.1 ± 0.1

n/c

DBM_10Å

-0.1 ± 0.0

-0.1 ± 0.0

-0.1 ± 0.0

-0.1 ± 0.0

n/c - not calculated

number of molecules in a cluster 2

4

6

8

10

12

14

16

18

20

0.00 -0.10 -0.20 fractionation factor, ε (‰)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-0.30 -0.40 -0.50 -0.60 -0.70 -0.80 -0.90 -1.00

Figure 6. 13C enrichment factors (ε, ‰) for evaporation of BB predicted using different size clusters and the PM7 Hamiltonian. The dark grey diamonds denote averaged over all carbon atomic positions magnitudes for each molecule in a cluster, the red ones are the averaged values over all molecules in a cluster, and the blues ones are the enrichment factors predicted for the central molecule in a cluster. The horizontal dashed line represents the experimental magnitude of isotopic fractionation for the system.



number of molecules in a cluster 2

4

6

8

10

12

14

16

18

20

-0.05 -0.15 fractionation factor, ε (‰)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 36

-0.25 -0.35 -0.45 -0.55 -0.65 -0.75 -0.85 -0.95

Figure 7. 81Br enrichment factors (ε, ‰) for evaporation of BB predicted using different size clusters and the PM7 Hamiltonian. The dark grey diamonds denote separate magnitudes for each molecule in a cluster, the red ones are the averaged values over all molecules in a cluster, and the blues ones are the enrichment factors predicted for the central molecule in a cluster. The horizontal dashed line representss the experimental magnitude of isotopic fractionation for the system.


Page 23 of 36

1.00

0.80

0.60 fractionation factor, ε (‰)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.40

0.20

0.00 2

4

6

8

10

12

14

16

18

20

-0.20

-0.40

Figure 8.

13

number of molecules in a cluster

C enrichment factors (ε, ‰) for evaporation of DBM predicted using different size

clusters and the PM7 Hamiltonian. The dark grey diamonds denote position magnitudes for each molecule in a cluster, the red ones are the averaged values over all molecules in a cluster, and the blues ones are the enrichment factors predicted for the central molecule in a cluster. The horizontal dashed line represents the experimental magnitude of isotopic fractionation for the system.



2

4

6

8

10

12

14

16

18

20

-0.05 -0.15 fractionation factor, ε (‰)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 36

-0.25 -0.35 -0.45 -0.55 -0.65 -0.75 -0.85

Figure 9.

79

number of molecules in a cluster

Br enrichment factors (ε, ‰) for evaporation of DBM predicted using different size

clusters and the PM7 Hamiltonian. The dark grey diamonds denote averaged over all bromine atomic positions magnitudes for each molecule in a cluster, the red ones are the averaged values over all molecules in a cluster, and the blues ones are the enrichment factors predicted for the central molecule in a cluster. The horizontal dashed line represents the experimental magnitude of isotopic fractionation for the system.

When we compare VPIEs resulted from the PIMD simulations (the best match with the experimental data: 13C VPIE of -3.7 and 2.8‰ and 81Br VPIE of 0.0 and -0.1‰ for BB and DBM, respectively) with the predictions resulted from either ONIOM and QM cluster calculations one may wonder whether we should consider the former approach a failure. First of all, it is noteworthy that within the PIMD calculations all molecules comprising each system (1757 and 503) in the case of BB and DBM, respectively were subjected to quantization, whereas a potential energy function was computed using the classical force field. It can be the first striking difference between this approach and any other, also employed in this work, as they all involved the theory level allowing to describe changes in electronic structure. As it may not seem to be important for the process under study but when we discuss isotope effects in general, one should keep in mind that they arise from changes in force constants which can be caused by electron density redistribution due to charge transfer, 24 ACS Paragon Plus Environment

Page 25 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for instance, not necessarily associated with bond breaking/forming processes, or change in energy levels. Therefore, in any attempt to predict isotope effects a quantum mechanical method ought to be rather used to take into account the quantum nature of electrons. In the PIMD simulation performed in the presented study, it was noticed that energy converges at a different rate for studied isotopes (81Br and

13

C). Bromine mean estimator convergence as a

function of imaginary slices number in the case of BB (Figure 4) reflects that highly likely GAFF parameters need further parametrization because only in the case when energy increases one obtaines better VPIEs predictions. Considering that IEs are local but they might occur from electron density redistribution, within the capabilities of AMBER package, a test in which a central molecule (DBM and BB) and its first solvation shell (based on our simulations) are represented by the PM3 Hamiltonian72 was performed with the aim of finding ways to correct/improve VPIE estimations. VPIE predictions coming from the aforementioned calculations (Table S12) turned out to be poorer than those from the classical PIMD. A similar behavior can be observed comparing results from PIMD29 and ab initio PIMD73 In general, PIMD provides better results when force field parameters are parametrized for the propose of reproducing experimental properties,74 whereas ab initio PIMD performing electronic structure calculation on the fly is subjected to errors owing to the nature of the ab initio part, for instance errors coming from the way electron correlation is included73,75 and definitely requires more testing and tuning with respect to theory level and system size before applying it to real problems. Having this in mind, in the current study we decided to make a compromise between the accuracy of theory level used for describing electrons and system size required for predicting bulk properties of a solvent. Better description of the quantum mechanical nature of the nuclei in isotopic substitution in the systems under study seemed to us a more crucial issue. Last, but not least it is worth mentioning that computational cost of PIMD calculations as compared to two-layer ONIOM or QM cluster calculations is much more higher, approximately two orders of magnitude (Table S13). The magnitudes of VPIEs coming from the measurements as well as computational models clearly indicate that the difference in carbon isotopic fractionation depends on whether it is related to the aliphatic or aromatic compound. Small inverse isotope effect observed and predicted for DBM can be easily explained by invoking non-polarity of a compound and the fact that the molecules of pure non-polar organic liquid mostly interact by van der Waals forces.

,77

BB, on the other hand, being only slightly more polar than DBM (dipole moment 25 ACS Paragon Plus Environment


of 1.70 D vs 1.43) is characterized by normal carbon VPIE at ambient temperature. Computed dipole polarizability for these two compounds indicates that even the system containing a single molecule of DBM is able to induce half of the instantaneous dipoles less than BB, 3.96 vs 8.73 Å3, respectively. Since the most polarized bond in either of the two molecules is the one between carbon and bromine atoms we also inspected the absolute difference of atomic charges sitting on it (their) ends. It turned out that in the case of BB this difference is 0.58 au on average, whereas in DBM it is 0.35 au. For both compounds, however, small normal 81Br VPIEs were measured and predicted. The isotopic fractionation associated with the dibromethane evaporation might be compared to the carbon and chlorine fractionation of dichloromethane reported by Huang et al. The direction of isotope effects obtained in that study is the same to the one resulted from the presented study. Depending on the experiment design one can expect a different set of data as it was also pointed out by the authors. They demonstrated for the studied systems that depending whether the measurement was taken under non-equilibrium conditions, when rather kinetic processes were governing evaporation, or at least partial equilibrium was established, one may expect different contribution to the overall isotope effect as the either set of data (kinetic vs equilibrium) can differ in both magnitude and in direction. Now, taking into account the current experimental set-up as well as judging based on the magnitude of obtained carbon isotopic fractionation of 0.6‰ for DBM evaporation under equilibrium conditions should rather operate in the presented study. Theoretical prediction of isotope effects associated with evaporation of dichloromethane (DCM) using the cluster comprising 12 molecules and the PM7 Hamiltonian leads to the same conclusions as for DBM (13C VPIE for DBM: 0.5 and for DCM: 0.6‰). Chlorine VPIE of -0.6‰ predicted for DCM evaporation has also the same direction as isotope effects measured and calculated by Huang et al. (-0.5‰).

Conclusions The current study presents experimental and theoretical evidence for carbon and for the first time reported bromine VPIEs for pure-phase organics. The difference in direction in fractionation with respect to carbon fractionation is demonstrated for the aliphatic and aromatic compound. Evaporation of both compounds resulted in normal bromine VPIE of the very similar magnitude.


Page 26 of 36

Page 27 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Computational prediction of measured isotope effects was also motivated by the need for testing different approaches and choosing the one that provides the best agreement with the experimentally obtained values. Comparing the PIMD to ONIOM and QM cluster results, it is clear that in general MO:MM and MO calculation provided better reproduction of the experimental trends even though the quantitative agreement has not been obtained in most of the cases. Nevertheless, we were able to find a model and an approach that would allow predicting the right direction of isotopic fractionation for both elements, carbon and bromine. However, in the case of carbon isotope effects we were able to demonstrate much better predictive capability than for bromine isotopic fractionation. None of the models used in the current study allowed to get quantitatively closer to the experimental values. In contrast to the carbon fractionation prediction which resulted in over- and underestimation, the bromine effects are always either in the wrong direction or underestimated, even for the position specific values. Although the cluster models provided at least the right direction of the observed effects, in this theoretical description some important components most likely are not included by default which could account for the magnitude of the measured effects. It must be pointed out that the main aim of MD simulations was to obtain reasonable approximations of the systems (DBM and BB) density. Deviations up to 9% of experimental density were observed on the BB model, whereas DBM gave percentages lower than 0.3%. Apart from density, computed radial distribution functions agreed very well with available experimental data obtained for BB. Having the right models of the systems under study was important for any subsequent calculations (either PIMD or MO:MM or MO optimizations) as they were all based on the configurations resulted from the MD simulations. A preliminary test of PIMD for VPIEs estimations reveals that it can be potentially used for such task even though such isotope effects are small. As the software package used for PIMD simulations only contains the most rudimentary implementation for performing VPIEs estimations, therefore it is suggested to aim for a rather more specialized/sophisticated tools not only for isotope effect calculations, but path integral as well, where features like ab initio and higher-order approximations might be included in the simulation customization. In the current study a compromise between the theory level for treating electronic structure (or rather neglecting it) and the model size for full quantization in order to predict nuclear



quantum effects on one hand and for bulk solvent properties on the other was made. This of course opens up a range of possible directions any future studies within this topic might take.

Acknowledgements This work was supported by the National Science Center in Poland (Sonata BIS grant UMO2014/14/E/ST4/00041) and in part by PLGrid Infrastructure (Poland). Supporting Information Force field parameters, atomic charges, averaged density for DBM and BB, additional radial distribution functions and mean estimator convergence in a function of imaginary slices.

References (1)

European Parliament, C. o. t. E. U. Eur-Lex - 32004l0042 - En - Eur-Lex , Limitation of Emissions of Volatile Organic Compounds Due to the Use of Organic Solvents in Certain Paints and Varnishes and Vehicle Refinishing Products and Amending Directive 1999/13/Ec_2004. 2004.

(2)

Huang, L.; Sturchio, N. C.; Abrajano, T.; Heraty, L. J.; Holt, B. D. Carbon and Chlorine Isotope Fractionation of Chlorinated Aliphatic Hydrocarbons by Evaporation.

Org. Geochem. 1999, 30, 777–785. (3)

Kornilova, A.; Huang, L.; Saccon, M.; Rudolph, J. Stable Carbon Isotope Ratios of Ambient Aromatic Volatile Organic Compounds. Atmos. Chem. Phys. 2016, 16, 11755–11772.

(4)

Kelley, C. A.; Hammer, B. T.; Coffin, R. B. Concentrations and Stable Isotope Values of BTEX in Gasoline-Contaminated Groundwater. Environ. Sci. Technol. 1997, 31, 2469–2472.

(5)

Kelly, C. A.; Rudd, J. W. M.; Bodaly, R. A.; Roulet, N. P.; St.Louis, V. L.; Heyes, A.; Moore, T. R.; Schiff, S.; Aravena, R.; Scott, K. J.; et al. Increases in Fluxes of Greenhouse Gases and Methyl Mercury Following Flooding of an Experimental Reservoir. Environ. Sci. Technol. 1997, 31, 1334–1344.

(6)

Aelion, C. M. Environmental Isotopes in Biodegradation and Bioremediation; 28 ACS Paragon Plus Environment

Page 28 of 36

Page 29 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Höhener, P., Hunkele, D., Aravena, R., Eds.; CRC Press Boca Raton, 2010. (7)

Eggenkamp, H. The Geochemistry of Stable Chlorine and Bromine Isotopes; 2014.

(8)

Winid, B. Bromine and Water Quality - Selected Aspects and Future Perspectives.

Appl. Geochemistry 2015, 63, 413–435. (9)

Flury, M.; Papritz, A. Bromide in the Natural Environment: Occurrence and Toxicity.

J. Environ. Qual. 1993, 22, 747. (10)

Gao, Y.; Yang, C.; Ma, J.; Yin, M. Characteristics of the Trace Elements and Arsenic, Iodine and Bromine Species in Snow in East-Central China. Atmos. Environ. 2018,

174, 43–53. (11)

Vonderheide, A. P.; Mueller, K. E.; Meija, J.; Welsh, G. L. Polybrominated Diphenyl Ethers: Causes for Concern and Knowledge Gaps Regarding Environmental Distribution, Fate and Toxicity. Sci. Total Environ. 2008, 400, 425–436.

(12)

Weijs, L.; Dirtu, A. C.; Das, K.; Gheorghe, A.; Reijnders, P. J. H.; Neels, H.; Blust, R.; Covaci,

A.

Inter-Species

Differences

for

Polychlorinated

Biphenyls

and

Polybrominated Diphenyl Ethers in Marine Top Predators from the Southern North Sea: Part 1. Accumulation Patterns in Harbour Seals and Harbour Porpoises. Environ.

Pollut. 2009, 157, 437–444. (13)

Bartrons, M.; Grimalt, J. O.; de Mendoza, G.; Catalan, J. Pollutant Dehalogenation Capability May Depend on the Trophic Evolutionary History of the Organism: PBDEs in Freshwater Food Webs. PLoS One 2012, 7.

(14)

El Megdiche, Y.; Ameur, W. Ben; Bèchir, H.; Hassine, S. Ben; Badreddine, B.; Touil, S.; Driss, M. R.; Eljarrat, E.; Barceló, D. Anthropogenic (PBDE) and NaturallyProduced (MeO-PBDE) Brominated Compound Levels in Bizerte Lagoon Clams ( Ruditapes Decussatus ): Levels and Human Health Risk Assessment. Mar. Pollut.

Bull. 2017, 125, 176–185. (15)

Le Breton, M.; Bannan, T. J.; Shallcross, D. E.; Khan, M. A.; Evans, M. J.; Lee, J.; Lidster, R.; Andrews, S.; Carpenter, L. J.; Schmidt, J.; et al. Enhanced Ozone Loss by Active Inorganic Bromine Chemistry in the Tropical Troposphere. Atmos. Environ. 2017, 155, 21–28.

(16)

Julien, M.; Parinet, J.; Nun, P.; Bayle, K.; Höhener, P.; Robins, R. J.; Remaud, G. S. 29 ACS Paragon Plus Environment


Page 30 of 36

Fractionation in Position-Specific Isotope Composition during Vaporization of Environmental Pollutants Measured with Isotope Ratio Monitoring by

13

C Nuclear

Magnetic Resonance Spectrometry. Environ. Pollut. 2015, 205, 299–306. (17)

Horst, A.; Lacrampe-Couloume, G.; Sherwood Lollar, B. Vapor Pressure Isotope Effects in Halogenated Organic Compounds and Alcohols Dissolved in Water. Anal.

Chem. 2016, 88, 12066–12071. (18)

Julien, M.; Höhener, P.; Robins, R. J.; Parinet, J.; Remaud, G. S. Position-Specific 13C Fractionation during Liquid-Vapor Transition Correlated to the Strength of Intermolecular Interaction in the Liquid Phase. J. Phys. Chem. B 2017, 121, 5810– 5817.

(19)

Dapprich, S.; Komáromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J. A New ONIOM Implementation in Gaussian98. Part I. The Calculation of Energies, Gradients, Vibrational Frequencies and Electric Field Derivatives. J. Mol. Struct.

THEOCHEM 1999, 461–462, 1–21. (20)

Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; 1965; Vol. 13.

(21)

Vanícek, J.; Miller, W. H.; Castillo, J. F.; Aoiz, F. J. Quantum-Instanton Evaluation of the Kinetic Isotope Effects. J. Chem. Phys. 2005, 123, 054108.

(22)

Zimmermann, T.; Vaníček, J. Path Integral Evaluation of Equilibrium Isotope Effects.

J. Chem. Phys. 2009, 131. (23)

Zimmermann, T.; Vaníček, J. Three Applications of Path Integrals: Equilibrium and Kinetic Isotope Effects, and the Temperature Dependence of the Rate Constant of the [1,5] Sigmatropic Hydrogen Shift in (Z)-1,3-Pentadiene. J. Mol. Model. 2010, 16, 1779–1787.

(24)

Lin, L. Density Functional Theory and Nuclear Quantum Effects, Princeton University: New Jersey, USA, 2011.

(25)

Ceriotti, M.; Markland, T. E. Efficient Methods and Practical Guidelines for Simulating Isotope Effects. J. Chem. Phys. 2013, 138, 014112.

(26)

Wong, K.-Y.; Xu, Y.; York, D. M. Ab Initio Path-Integral Calculations of Kinetic and Equilibrium Isotope Effects on Base-Catalyzed RNA Transphosphorylation Models. J.

Comput. Chem. 2014, 35, 1302–1316. 30 ACS Paragon Plus Environment

Page 31 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(27)

Del Buono, G. S.; Rossky, P. J.; Schnitker, J. Model Dependence of Quantum Isotope Effects in Liquid Water. J. Chem. Phys. 1991, 95, 3728–3737.

(28)

De La Peña, L. H.; Kusalik, P. G. Quantum Effects in Light and Heavy Liquid Water: A Rigid-Body Centroid Molecular Dynamics Study. J. Chem. Phys. 2004, 121, 5992– 6002.

(29)

Markland, T. E.; Berne, B. J. Unraveling Quantum Mechanical Effects in Water Using Isotopic Fractionation. Proc. Natl. Acad. Sci. 2012, 109, 7988–7991.

(30)

Cheng, B.; Ceriotti, M. Direct Path Integral Estimators for Isotope Fractionation Ratios. J. Chem. Phys. 2014, 141.

(31)

Frisch, M. J.; Trucks, W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision E.01. Gaussian, Inc.: Wallingford, CT. 2009.

(32)

Vosko, S. H.; Wilk, L.; Nusair, M. Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: A Critical Analysis. Can. J.

Phys. 1980, 58, 1200–1211. (33)

Lee, C.; Yang, W.; Parr, R. Development of the Colle- Salvetti Correlation Energy Formula into a Functional of the Electron Density. Phys Rev B 1988, 37, 785–789.

(34)

Becke, A. Density Functional Thermochemistry III The Role of Exact Exchange. J.

Chem. Phys. 1993, 98, 5648–5652. (35)

Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623–11627.

(36)

Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Function. Theor. Chem. Acc. 2008, 120, 215– 241.

(37)

Møller, C.; Plesset, M. S. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618–622.



(38)

Head-Gordon, M.; Pople, J. A.; Frisch, M. J. MP2 Energy Evaluation by Direct Methods. Chem. Phys. Lett. 1988, 153, 503–506.

(39)

Čížek, J. On the Correlation Problem in Atomic and Molecular Systems. Calculation of Wavefunction Components in Ursell‐Type Expansion Using Quantum‐Field Theoretical Methods. J. Chem. Phys. 1966, 45, 4256–4266.

(40)

Purvis, G. D.; Bartlett, R. J. A Full Coupled-Cluster Singles and Doubles Model: The Inclusion of Disconnected Triples. J. Chem. Phys. 1982, 76, 1910–1918.

(41)

Petersson, G. A.; Bennett, A.; Tensfeldt, T. G.; Al-Laham, M. A.; Shirley, W. A.; Mantzaris, J. A Complete Basis Set Model Chemistry. I. The Total Energies of ClosedShell Atoms and Hydrides of the First-Row Elements. J. Chem. Phys. 1988, 89, 2193– 2218.

(42)

Petersson, G. A.; Al-Laham, M. A. A Complete Basis Set Model Chemistry. II. OpenShell Systems and the Total Energies of the First-Row Atoms. J. Chem. Phys. 1991,

94, 6081–6090. (43)

Singh, U. C.; Kollman, P. A. An Approach to Computing Electrostatic Charges for Molecules. J. Comput. Chem. 1984, 5, 129–145.

(44)

Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic Charges Derived from Semiempirical Methods. J. Comput. Chem. 1990, 11, 431–439.

(45)

Dupradeau, F.-Y.; Pigache, A.; Zaffran, T.; Savineau, C.; Lelong, R.; Grivel, N.; Lelong, D.; Rosanski, W.; Cieplak, P. The R.E.D. Tools: Advances in RESP and ESP Charge Derivation and Force Field Library Building. Phys. Chem. Chem. Phys. 2010,

12, 7821. (46)

Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157–1174.

(47)

Vanquelef, E.; Simon, S.; Marquant, G.; Garcia, E.; Klimerak, G.; Delepine, J. C.; Cieplak, P.; Dupradeau, F. Y. R.E.D. Server: A Web Service for Deriving RESP and ESP Charges and Building Force Field Libraries for New Molecules and Molecular Fragments. Nucleic Acids Res. 2011, 39.

(48)

Case, D. A.; Darden, T.; Iii, T. E. C.; Simmerling, C.; Brook, S.; Roitberg, A.; Wang, J.; Southwestern, U. T.; Duke, R. E.; Hill, U.; et al. Amber 14. University of 32 ACS Paragon Plus Environment

Page 32 of 36

Page 33 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


California, San Francisco. University of California 2014. (49)

Uberuaga, B. P.; Anghel, M.; Voter, A. F. Synchronization of Trajectories in Canonical

Molecular-Dynamics

Simulations:

Observation,

Explanation,

and

Exploitation. J. Chem. Phys. 2004, 120, 6363–6374. (50)

Sindhikara, D. J.; Kim, S.; Voter, A. F.; Roitberg, A. E. Bad Seeds Sprout Perilous Dynamics:

Stochastic

Thermostat

Induced

Trajectory

Synchronization

in

Biomolecules. J. Chem. Theory Comput. 2009, 5, 1624–1631. (51)

Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods V: Modification of NDDO Approximations and Application to 70 Elements. J. Mol.

Model. 2007, 13, 1173–1213. (52)

Berne, B. J.; Thirumalai, D. On the Simulation of Quantum Systems: Path Integral Methods. Annu. Rev. Phys. Chem. 1986, 37, 401–424.

(53)

Vaníček, J.; Miller, W. H. Efficient Estimators for Quantum Instanton Evaluation of the Kinetic Isotope Effects: Application to the Intramolecular Hydrogen Transfer in Pentadiene. J. Chem. Phys. 2007, 127.

(54)

Voth, G. A.; Chandler, D.; Miller, W. H. Rigorous Formulation of Quantum Transition State Theory and Its Dynamical Corrections. J. Chem. Phys. 1989, 91, 7749–7760.

(55)

Haynes, W. M. CRC Handbook of Chemistry and Physics, 95th Edition, 2014-2015; 2014; Vol. 54.

(56)

Humphrey, W.; Dalke, A.; Schulten, K. {VMD} -- {V}isual {M}olecular {D}ynamics. J. Mol. Graph. 1996, 14, 33–38.

(57)

Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods VI: More Modifications to the NDDO Approximations and Re-Optimization of Parameters. J.

Mol. Model. 2013, 19, 1–32. (58)

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; et al. Gaussian 16, Revision A.03. Gaussian, Inc.: Wallingford, CT. 2016.

(59)

Bigeleisen, J. The Relative Reaction Velocities of Isotopic Molecules. J. Chem. Phys. 1949, 17, 675–678. 33 ACS Paragon Plus Environment


(60)

Anisimov, V.; Paneth, P. ISOEFF98. A Program for Studies of Isotope Effects Using Hessian Modifications. J. Math. Chem. 1999, 26, 75–86.

(61)

Chemicalbook. Dibromomethane. November 2016.

(62)

Gotoh, K. Particulate Morphology: Mathematics Applied to Particle Assemblies; 2012.

(63)

Bertagnolli, H.; Engelhardt, T.; Lengeler, B. Local Order in Liquid Bromobenzene EXAFS Measurement and Its Combination with X-Ray and Neutron Diffraction Results. Berichte der Bunsengesellschaft für Phys. Chemie 1987, 91, 1287–1292.

(64)

Martin, F.; Zipse, H. Charge Distribution in the Water Molecule -A Comparison of Methods. J. Comput. Chem. 2005, 26, 97–105.

(65)

Cao, J.; Berne, B. J. On Energy Estimators in Path Integral Monte Carlo Simulations: Dependence of Accuracy on Algorithm. J. Chem. Phys. 1989, 91, 6359–6366.

(66)

Ceperley, D. M. Path Integrals in the Theory of Condensed Helium. Rev. Mod. Phys. 1995, 67, 279–355.

(67)

Markland, T. E.; Manolopoulos, D. E. An Efficient Ring Polymer Contraction Scheme for Imaginary Time Path Integral Simulations. J. Chem. Phys. 2008, 129.

(68)

Nakanishi, K.; Solomon, P. H. Infrared Absorption Spectroscopy; 1977; Vol. 26.

(69)

J. Gutow, A. H. Bromobenzene. 2014.

(70)

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.

(71)

Hostaš, J.; Řezáč, J.; Hobza, P. On the Performance of the Semiempirical Quantum Mechanical PM6 and PM7 Methods for Noncovalent Interactions. Chem. Phys. Lett. 2013, 568–569, 161–166.

(72)

Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods I. Method. J.

Comput. Chem. 1989, 10, 209–220. (73)

Wang, L.; Ceriotti, M.; Markland, T. E. Quantum Fluctuations and Isotope Effects in Ab Initio Descriptions of Water. J. Chem. Phys. 2014, 141.

(74)

Marsalek, O.; Chen, P. Y.; Dupuis, R.; Benoit, M.; Méheut, M.; Bačić, Z.; Tuckerman, 34 ACS Paragon Plus Environment

Page 34 of 36

Page 35 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


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. (75)

Machida, M.; Kato, K.; Shiga, M. Nuclear Quantum Effects of Light and Heavy Water Studied by All-Electron First Principles Path Integral Simulations. J. Chem. Phys. 2018, 148, 102324.

(76)

Jancso, G.; Van Hook, W. A. Condensed Phase Isotope Effects. Chem. Rev. 1974, 74, 689–750.

(77)

Wolfsberg, M.; Alexander Van Hook, W.; Paneth, P.; Rebelo, L. P. N. Isotope Effects:

In the Chemical, Geological, and Bio Sciences; 2009.

TOC Graphic




Page 36 of 36

Can Path Integral Molecular Dynamics Make a Good Approximation

Recommend Documents