Reaction ensemble Monte Carlo simulation of xylene isomerization in

Subscriber access provided by UNIVERSITY OF THE SUNSHINE COAST

Article

Reaction ensemble Monte Carlo simulation of xylene isomerization in bulk phases and under confinement Ryan G. Mullen, and Edward J. Maginn J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00498 • Publication Date (Web): 27 Jul 2017 Downloaded from http://pubs.acs.org on July 29, 2017

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 free 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 accessible to all readers and 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.

Journal of Chemical Theory and Computation 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 31

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

Journal of Chemical Theory and Computation

Reaction ensemble Monte Carlo simulation of xylene isomerization in bulk phases and under confinement Ryan Gotchy Mullen and Edward J. Maginn∗ Department of Chemical and Biomolecular Engineering, University of Notre Dame, Notre Dame, IN E-mail: [email protected]

Abstract The original reaction move for the reaction ensemble Monte Carlo (RxMC) method is adapted to align both the position and orientation of inserted product molecules and deleted reactant molecules. The accuracy and efficiency of this move is demonstrated for xylene isomerization in vapor, liquid, and supercritical phases. Classical RxMC requires the ideal gas free energy of reaction ∆Gideal rxn as an input. We compare three methods for computing ∆Gideal rxn : using tabulated enthalpies and entropies of formation, using the harmonic oscillator and rigid rotor approximations, and using QM/MM alchemical transformation combined with multistate Bennett acceptance ratio. We find that the tabulated free energies of reaction give the best agreement with experimental equilibrium compositions in bulk fluids. RxMC simulations in a carbon nanotube with an inner diameter of approximately 6 Å show that p-xylene becomes the dominant isomer under confinement, an effect consistent with the production of p-xylene in the zeolite ZSM-5. We also show that o-xylene becomes the dominant isomer in nanotubes with an inner diameter of 7-8 Å. We find that both m- and p-xylene exhibit a loss of

1

ACS Paragon Plus Environment


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

rotational entropy in nanotubes of this diameter, effectively allowing o-xylene to fit into cavities inaccessible to the other isomers.

1

Introduction

Chemists and chemical engineers are often interested in the equilibrium behavior of complex chemical systems, including multiphase systems, non-ideal liquid solutions, and systems under nanoscale confinement. Measuring equilibrium compositions can be a time-consuming laborious process even under mild conditions. At high temperature or high pressure or for reactions involving hazardous chemicals, specialized equipment and training may be required. In other cases it may not be possible to collect fluid samples, e.g. from inside porous catalysts. Robust simulation methods would greatly facilitate the prediction of chemical equilibria under these circumstances. Including chemical reactivity in molecular simulations of condensed phases faces a number of challenges. In simulations employing traditional classical force fields, covalent bonds are modeled harmonically and so cannot be formed or broken. Quantum mechanical methods, on the other hand, can explicitly model the changing electronic structure as molecules react but are typically limited to systems of O(10 − 100) atoms. 1 A solvent environment of O(100 − 1000) classical molecules can be used to surround the quantum cluster in a QM/MM simulation. 2 However, since the QM region must be specified a priori, only the few QM molecules can react throughout the simulation. With a classical, reactive force field like ReaxFF, 3 we can model O(100 − 1000) molecules without any restrictions on reactivity. Even so, in a brute force molecular dynamics simulation observing even a single reactive event would be rare (except at elevated temperatures 4 ) due to the free energy barrier separating reactants and products. The waiting time to observe several hundred reactive events would clearly exceed current computational resources. Using reaction ensemble Monte Carlo (RxMC), 5–7 we can avoid sampling configurations

2


Page 2 of 31

Page 3 of 31

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


at the top of a free energy barrier—so called transition states—by using carefully designed reactive moves that stochastically transform reactant molecules into product molecules (and vice versa). The reaction ensemble is defined by specifying the temperature, either the volume or the pressure, and an initial composition, including the excess amount of any nonlimiting reagents. For classical RxMC simulations, the set of reactions also needs to be prespecified. For each reaction the stoichiometry and either the change in the free energy of ideal an ideal gas due to the reaction ∆Gideal , need to rxn or the ideal gas equilibrium constant K

be supplied. The ideal gas parameter is needed because, using a classical model, information about the relative energy and entropy of each species has been omitted. When the reactions being modeled involve well-studied chemicals, ∆Gideal rxn can be computed from tabulated gas phase free energies of formation. 8 Otherwise, ∆Gideal rxn can be estimated by computing the ab initio free energy of each species from the potential energy of an optimized configuration and estimating the free energy contributions using the harmonic oscillator and rigid rotor approximations. 9 When the potential energy surface is anharmonic in the region near the minimum energy configuration or when multiple minima make significant contributions to the molecular partition function—as with longer chain hydrocarbons, polymers, or proteins that have multiple stable conformations—it is difficult or impractical to compute the partition function directly. The free energy difference between reactants and products can instead be estimated using molecular simulations that reversibly connect the reactant and product states. One such method, QM/MM alchemy, has recently been combined with the multistage Bennett Acceptance Ratio (MBAR) to compute differences in free energies of solvation between methanol and ethane in water. 10,11 All three methods will be used and compared in this work. In the original reaction move, each product molecule is inserted at the same position from which a reactant is deleted but with a random conformation and orientation. 5–7,12–14 If the reaction changes the total number of molecules, unpaired molecules are inserted at a random location or deleted without replacement. The efficiency of reaction moves that

3



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

require random insertions is increased by using configurational bias Monte Carlo. 13,14 There are two alternatives to inserting and deleting molecules in a single step. In continuous fractional component (CFC) MC, molecules are inserted in stages and the surrounding fluid is allowed to adjust as the molecule is gradually turned on. 15–18 However, the free energy penalty of solvent reorganization can be so large that the number of stages required to insert a molecule becomes prohibitive. This is the case, for example, if a reacting solute is much larger than the solvent molecules or if the solvent exhibits a nanostructure 19 that is disrupted by the inserted molecule. Additionally, inserting molecules via CFC will fail for confined fluids, as studied here, since fractional molecules can move through the walls that confine the fluid. In reactive first-principles MC, molecules are clusters of independent atomic nuclei that can be broken and reformed through a specialized set of MC moves. 1 No chemical reactions need to be specified a priori and since atomic clusters (i.e. molecules) are not inserted or deleted, no ideal gas parameters are required. However, reactive first-principles Monte Carlo faces two challenges that currently limit its utility: (1) developing specialized MC moves that can transform reactants into products can be a significant challenge for polyatomic molecules, and (2) the computational time required to calculate the ab initio potential energy limits the simulation size and duration. Rather than use CFC or first-principles MC, we employ a simple adaptation of the original move. In addition to inserting a product molecule at the same location as a deleted reactant, we also align its orientation to that of the reactant molecule. We demonstrate this adapted move for xylene isomerization, shown in Fig. 1, in bulk liquid, vapor, and supercritical phases. In particular, p-xylene is a valuable commodity chemical precursor for the production of PETE synthetic fibers and plastics. However, p-xylene is only a minority component in an equilibrium mixture of xylenes across a wide range of temperatures. 20 Using the zeolite ZSM-5 to catalyze toluene methylation or toluene disproportionation selectively produces p-xylene. 21–25 The shape selective properties of ZSM-5 stem from its unique pore structure: straight elliptical channels 5.4-5.6 Å in diameter criss-cross with similarly sized

4


Page 4 of 31

Page 5 of 31

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


sinusoidal channels to form cavities approximately 8 Å in diameter at the intersection. 26 It is hypothesized that an equilibrium mixture of xylenes forms at the intersections, but that m- and o-xylenes are not able to exit the cavity through the connecting channels. 27 We propose to study the effect of confinement on the equilibrium composition of xylene mixtures inside carbon nanotubes. Carbon nanotubes are an idealized model of confinement easily created with varying diameters, whereas adjusting the complex three-dimensional structure of zeolites is more difficult. Previous studies have shown that RxMC can capture the shift in equilibrium that confinement induces towards the state with the fewest molecules. 13,28 We will demonstrate that confinement can also shift equilibrium for unimolecular reactions in complex ways.

Figure 1: Molecular structures of para-, meta- and ortho-xylene.

2

Theory

The reaction ensemble is an open ensemble in which fluctuations in the number of molecules of each species are related by the stoichiometry of a reaction. We define the stoichiometric coefficients νs for species s to be negative for reactants and positive for products. The change P in the total number of molecules for a forward reaction move is therefore ∆ν = νs . In

addition to the stoichiometry, the reaction ensemble is defined by specifying the ideal gas change in free energy due to the reaction ∆Gideal rxn at standard pressure, the temperature T , either the volume V or pressure P , and the excess number of molecules of non-limiting reactants. 5



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 6 of 31

Thermal equilibrium is maintained by moves that translate, rotate, and regrow molecules using configurational bias. For simulations at constant pressure, mechanical equilibrium is maintained by changing the volume and scaling the center-of-mass position of each molecule. For simulations with multiple phases (i.e. multiple boxes), chemical equilibrium is maintained by swapping molecules between the boxes. Acceptance rules for each of these moves are the same as previously published. 29 The Metropolis acceptance criterion for a reaction move can be derived within the osmotic ensemble. 16,30,31 The probability of a microstate m in the osmotic ensemble is

pm ∝ exp −βUm − βP Vm + β

X

µs Ns

s

!

dV

Y

dq

s

Z

Ns exp(−βKs )dp

(1)

where β is the inverse temperature 1/kB T , Um is the potential energy of the microstate, Vm is the volume of the microstate, Ns is the number of molecules of species s in the microstate, Ks is the kinetic energy of a single molecule, and q and p are vectors of the generalized coordinates and momenta of each atom within a molecule. The differential elements have been included to make pm dimensionless and the momenta integrals have been decoupled. For a microstate n generated from microstate m by attempting a forward reaction move, the number of molecules of each species differs by the reaction stoichiometry, ∆Ns = νs . The use of forward and back reaction moves allows the system to relax to reaction equilibrium, P such that µs νs = 0. Lastly, the momentum integrals for each molecule are multiplied by

the single molecule configurational partition function Zsideal to give the ideal gas partition

function Qideal . The product of partition functions are replaced by the ideal gas equilibrium s constant K ideal or equivalently the ideal gas free energy of reaction ∆Gideal rxn according to Y Qideal νs s

s

V

=K

ideal

= exp

−β∆Gideal rxn

P◦ kB T

∆ν

(2)

where P o is the standard state pressure (typically 1 atm or 1 bar). The ratio of probabilities

6


Page 7 of 31

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


between microstates n and m is therefore pn = exp −β(Un − Um ) − β∆Gideal rxn pm

P oV kB T

∆ν Y s

dq Zsideal

ν s

(3)

Zsideal is separable into translational, rotational, and intramolecular components, Zsideal = V Zsrot Zsintra , where Zsrot is 4π for a linear molecule or 8π 2 for a non-linear molecule. At liquid densities the change in energy Un − Um from deleting reactants and inserting products in a random location is often much larger than kB T . If ∆ν = 0, then each product molecule can be inserted into a cavity left by deleting a reactant. Otherwise, any unpaired molecules need to be inserted randomly or deleted without replacement. In the following semi-grand algorithm, we will refer to product molecule i with configuration j and location and orientation k and reactant molecule i′ with configuration j ′ at location and orientation k ′ . 1. Select |νs | unique molecules of reactant species s with uniform probability. The probability of picking a collection of |νs | molecules is Psmol =

(Ns − |νs |)! Ns !

(4)

If there are less than |νs | molecules in the box, the move is rejected. 2. For each product molecule i, choose an ideal gas configuration j with Boltzmann probability exp(−βUijideal )/Ziintra . 3. Insert molecule i: • If molecule i is replacing reactant molecule i′ , align the coordinates of atom 1 on each molecule by setting r1i = r1i′ with probability Piins = 1. • Otherwise, choose κins trial locations at which to insert the molecule’s center of mass. Each location is chosen with probability Piins = dr/V . 4. If molecule i has 2 or more atoms, rotate it: 7



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 8 of 31

• If molecule i is replacing molecule i′ and i′ has at least two atoms, rotate molecule i to align vectors r2i − r1i and r2i′ − r1i′ with Pirot1 = 1. • Otherwise, choose trial orientations randomly each with probability Pirot1 = dθdφ/4π. 5. If molecule i is non-linear, rotate it about axis r2i − r1i : • If molecule i is replacing molecule i′ and i′ is non-linear, rotate i such that the atoms 1, 2, and 3 on molecule i are coplanar with atoms 1, 2, and 3 on molecule i′ with Pirot2 = 1. • Otherwise, rotate trial orientations randomly each with probability Pirot2 = dψ/2π. If the rotation is random, choose a total of κrot trials at each location. 6. If there is more than one trial, compute the change in energy ∆Uijk′′ from inserting molecule i at each trial (location and orientation) k ′′ . Choose one trial location and orientation k with probability exp(−β∆Uijk ) bias Pijk = PN trials i k′′ =1 exp (−β∆Uijk′′ )

(5)

The forward attempt probability αmn is

αmn

Y exp −βUijideal dqi bias Piins Pirot1 Pirot2 Nitrials Pijk = Psmol intra Z i i∈P s∈R Y

(6)

where the first product is over reactant species and the second product is over product molecules. The probability αnm for a reverse reaction move is computed similarly to ensure that the reaction move follows detailed balance. For a unimolecular reaction A ⇀ ↽ B, the forward move consists of aligning the inserted B

8


Page 9 of 31

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


with a deleted A. The acceptance rule is

acc(m → n) =

min 1, exp −β (Un − Um ) −

β∆Gideal A⇀ ↽B

−β

Uiideal ′j′

−

Uijideal

NA NB + 1

(7)

Metropolis acceptance criteria for additional reactions are derived in the SI.

3

Methods

3.1

Free energy calculations

The change in free energy of an ideal gas due to the reaction, ∆Gideal rxn , was computed three ways. First, the free energy was computed from tabulated thermodynamic properties. The enthalpy H and entropy S of formation of each isomer in the gas phase at 298 K and 1 atm are given in Table S1. 32,33 The gas phase heat capacities Cp as a function of temperature are given in Table S2. 34,35 The free energy of reaction was then calculated from 300 to 1000 K, at 100 K intervals, ∆Gtab rxn

=

◦

∆Hrxn +

Z

T

∆Cp dT 298K

−T

◦

∆Srxn +

Z

T 298K

∆Cp dT T

(8)

where the superscript tab indicates the free energy is computed from tabulated data. Second, a geometry optimization and frequency analysis was performed on each xylene isomer using Gaussian 09 36 at the MP2 level of theory and using a 6-311+G(2d,2p) basis set. The ideal gas free energy of each molecule was computed at 1 atm and temperatures up to 600 K using the harmonic oscillator and rigid rotor approximations. The harmonic contribution of each normal mode was computed by scaling the frequencies by 0.95, per NIST recommendations. 37 Corrections to the free energy were also computed by treating the two lowest frequency normal modes as hindered rotors. 38,39

9



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

Third, m-xylene was alchemically transformed to o- or p-xylene. The m- to o- transformation will be described in detail; the m- to p- is analogous. Each simulation contained a single dual-topology molecule consisting of 23 atoms as shown in Fig. S1. The potential energy functions, bonded and nonbonded terms, are modified using a parameter λ that defines the alchemical stage. At intermediate stages, 0 < λ < 1, carbons 2 and 3 are each bonded to both a hydrogen atom and a methyl group. As λ → 0, one hydrogen and one methyl group are transformed into an ideal gas state, freely traversing the simulation cell with no interactions, and the remaining atoms form a molecule of m-xylene. As λ → 1, the complementary hydrogen and methyl become ideal gas atoms, leaving a molecule of o-xylene. Molecular dynamics trajectories were computed using NAMD 40 and the OPLS force field. 41 The timestep was 1 fs and configurations were saved every 100 fs over the 5 ps trajectory. The temperature was maintained by applying a Langevin thermostat to the carbon atoms with a damping coefficient of 5.0 / ps. Electrostatic interactions for alchemical atoms were turned off over the range 0.0 ≤ λ ≤ 0.5, while van der Waals and bonded interactions were turned off over the range 0.0 ≤ λ ≤ 0.7. Additionally, the quantum energies of the 50 thousand saved configurations at 6 λ stages were computed using Q-CHEM 42 at a B3LYP level of theory with 6-31G* basis set. For λ = 0.0, 0.00001, and 0.00002, only the coordinates corresponding to a molecule of m-xylene were used. For λ = 0.99998, 0.99999, and 1.0, only the coordinates corresponding to a molecule of o-xylene were used. The free energy difference between m- and o-xylene was calculated using 3-stage MBAR. 11,43

3.2

Monte Carlo Simulations

All reaction ensemble Monte Carlo (RxMC) simulations were run using Cassandra 29 and the TraPPE force field. 44 The potential energy was calculated using periodic boundary conditions, Lennard-Jones interactions truncated at 14.0 Å, and tail corrections. Four series of simulations were run: bulk liquid and vapor phases at 1 bar and varying temperatures, bulk supercritical fluids at 700 K and varying densities, a single molecule inside nanotubes 10


Page 10 of 31

Page 11 of 31

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


of varying diameters, and fluids inside nanotubes of varying diameters in equilibrium with a bulk vapor at 300 K and 0.01 bar. Bulk liquid and vapor simulations were run at 1 bar and each of eight temperatures, ranging from 300 to 1000 K in 100 K increments. Each was initiated from 300 m-xylene molecules placed randomly in a cubic simulation cell. At temperatures below 412 K, the normal boiling point of a xylene mixture, the simulation cell was initially 42 Å on a side. At higher temperatures, the initial simulation volume was determined by the ideal gas law. Each MC sweep consisted of, on average, 300 molecule translations, 300 rotations, 1 volume change, and 10 reaction moves. The maximum displacement, rotation, and volume change was adjusted during equilibration to achieve 50% acceptance. During a reaction move, the carbon ring of the inserted xylene was aligned with that of the deleted xylene using carbons 1, 2, and 6. Unless specified otherwise, the ideal gas free energy of reaction from tabulated thermodynamic properties was used. Each simulation was equilibrated for 50 thousand sweeps, and then averages were taken over an additional 50 thousand sweeps with compositions sampled every 100 sweeps. Bulk supercritical fluid simulations were run at 700 K and each of nine densities: 1, 10, and 100 to 700 kg/m3 in 100 kg/m3 increments. An initial configuration of 300 m-xylene molecules was equilibrated for 10 thousand sweeps at each density with no reaction moves. Then reaction moves were turned on (10 per sweep) and acceptance rates were computed over an additional 10 thousand sweeps. Three simulations were run at each thermodynamic state point to test different insertion methods: random position and orientation, aligned position and random orientation, aligned position and orientation. Single-walled zigzag carbon nanotubes were generated using VMD 45 with the diameters shown in Table. 1. The carbon nanotube atoms were held fixed in place, so wall-wall interactions were ignored but LJ interactions with xylene molecules were computed using Lorentz-Berthelot mixing rules and wall parameters ǫ = 0.1547 kJ/mol and σ = 3.602 Å. 46 Nanotubes were centered on the z-axis and x- and y-dimensions of the box were large enough

11



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

that periodic images in those directions could be ignored. Interior xylene molecules were translated and rotated with maximum displacements of 0.25 Å and 20◦ , respectively. Table 1: Dimensions, in Å, of the single-walled zigzag carbon nanotubes.

b

(m,n) indices (12,0) (13,0) (14,0) (15,0) (16,0) (17,0) (18,0) (24,0) a diameter 9.4 10.2 11.0 11.7 12.5 13.3 14.1 18.8 b inner diameter 5.8 6.6 7.4 8.1 8.9 9.7 10.5 15.2 length 782.7 608.3 489.2 399.9 331.8 280.8 242.5 119.1 a center-to-center distance of opposite carbon atoms in the xy-plane estimated by subtracting the LJ diameter of a carbon atom, 3.6 Å, from the nanotube diameter Simulations with a single xylene molecule were run in nanotubes 34.0 Å long. One m-

xylene molecule was initially placed in the nanotube and equilibrated for 10 thousand sweeps. Each MC sweep consisted of, on average, 10 translations, 10 rotations and 1 reaction move. Composition was sampled every 10 sweeps over a 990 thousand sweep production run. Two-box simulations were used to compute the equilibrium between a bulk vapor at 300 K and 0.01 bar and the interior of the carbon nanotube. As in the Gibbs ensemble, the chemical potentials of each isomer are unknown but equilibrated between phases by swapping molecules from the vapor to the nanotube interior and vice versa. Volume moves are only attempted for the vapor box. Reaction moves are attempted in both boxes. The vapor pressure was set to 0.01 bar to facilitate swaps between phases while still resulting in liquid-like densities inside the nanotubes. The initial vapor configuration of 300 xylene molecules was equilibrated for 50 thousand sweeps using thermal, volume, and reaction moves. The nanotube was initially empty. The length of each nanotube, shown in Table 1, was chosen such that the interior volume was approximately 100 times the molecular volume of liquid xylene at the same temperature, 207 Å3 . Each MC sweep consisted of, on average, 300 translations, 300 rotations, 1 volume move for the vapor box, 10 reaction moves and 1,000 molecule swap moves. The system was equilibrated for 50 thousand sweeps and then compositions were sampled every 100 sweeps over a 50 thousand sweep production run.

12


Page 12 of 31

Page 13 of 31

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


4

Results and Discussion

4.1

Ideal Gas Free Energies

The ideal gas change in free energy due to the reaction as a function of temperature is shown in Fig. 2. Computed free energies will be distinguished using the subscripts m ⇀ ↽ p and m ⇀ ↽o to specify the reaction, and superscripts to specify the method: tab for tabulated values, HO/RR for harmonic oscillator and rigid rotor approximations, and alch for alchemical transformation. The free energy ∆Gtab m⇀ ↽p is approximately linear since the enthalpy and entropy of reaction have little temperature dependence. ∆Gtab m⇀ ↽o , on the other hand, is non-linear, reflecting the differences in heat capacity between o- and m-xylene. The higher heat capacity of o-xylene tab causes ∆Gtab m⇀ ↽o to drop below ∆Gm⇀ ↽p at temperatures above 600 K. Uncertainty in the ideal

gas free energies of reaction at 298 K due to experimental uncertainty is 1.4 kJ/mol. The harmonic oscillator and rigid rotor approximations yield free energies that are linear HO/RR

with temperature and qualitatively show that ∆Gm⇀ ↽o

HO/RR

> ∆Gm⇀ ↽p

HO/RR

> 0. While ∆Gm⇀ ↽p

HO/RR

underpredicts ∆Gtab ↽o m⇀ ↽p and is nearly constant with temperature, ∆Gm⇀

matches ∆Gtab m⇀ ↽o

at 300 K but overpredicts the free energy at higher temperatures. Using the HO/RR method, contributions to the free energy from the difference in potential energy of the optimized configurations, and translational, rotational, and vibrational degrees of freedom are easily separated (see Tables S3-S6). There is no translational contribution to either free energy difference since all isomers have the same mass. The rotational contribution depends on the symmetry of each isomer, which in turn depends on the relative orientation of the exo-ring methyl groups. The optimized m- and p-xylenes have staggered methyl groups and therefore C2 point group symmetry. The optimized o-xylene has antiplanar methyl groups, with one hydrogen atom of each methyl group in the plane of the aromatic ring and pointing away from the other methyl, consistent with prior results. 47 This configuration has C2V point group HO/RR

symmetry. The rotational contributions to ∆Gm⇀ ↽o 13

HO/RR

and ∆Gm⇀ ↽p


are less than 0.5 kJ/mol


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 14 of 31

HO/RR

at all temperatures. Other contributions to ∆Gm⇀ ↽p —differences in potential energy and vibrational—are also less than 0.25 kJ/mol in magnitude. The difference in potential energy of the optimized o- and m- configurations is -1.9 kJ/mol. The largest contribution HO/RR

to ∆Gm⇀ ↽o

is from the vibrational degrees of freedom, 4.6 kJ/mol at 300 K increasing to

9.5 kJ/mol at 600 K. Only the two lowest frequency vibrational modes make non-negligible contributions to the free energy. From visual inspection these two modes are symmetric and asymmetric rotations of the exo-ring methyl groups. Modeling these two modes as hinHO/RR

dered rotors instead of harmonic oscillators makes negligible corrections to ∆Gm⇀ ↽p HO/RR

significantly lowers ∆Gm⇀ ↽o

but

such that it becomes less than ∆Gtab m⇀ ↽o at all temperatures.

Sensitivities run with the same 6-311+G(2d,2p) basis set but at an B3LYP level of theory show that the MP2 results are more accurate for both reactions at all temperatures (see Table S7). alch The free energies computed using alchemical transformation, ∆Galch m⇀ ↽p and ∆Gm⇀ ↽o , do

not exhibit a linear or even monotonic temperature dependance. Free energy differences between individual λ stages are shown in Tables S8 and S9. The maximum error |∆Galch rxn − ⇀ ⇀ ∆Gtab rxn | is 2.2 kJ/mol and 6.0 kJ/mol for the m ↽ p and m ↽ o reactions, respectively. Since these errors are the same magnitude as ∆Gtab rxn , QM/MM alchemy should only be used to provide order of magnitude estimates of reaction free energies. Sensitivities run for the m ⇀ ↽o reaction with the same 6-31G* basis set but at an MP2 level of theory exhibit comparable accuracy and also lack a monotonic temperature dependence (see Table S10).

4.2

Bulk Fluids

The computed and experimental 20 equilibrium compositions of a mixture of xylene isomers at 1 bar and a range of temperatures is shown in Fig. 3. m-Xylene is the dominant isomer at all conditions and the temperature dependence of the equilibrium mole fraction is well captured by our simulations. The population of p-xylene shows little temperature dependence, either experimentally or in simulation. The mole fraction of o-xylene increases with temperature, 14


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


Free Energy of Rxn, kJ/mol

Page 15 of 31

10 8 6 m→o

4

m→p

2 0 300

400 500 Temperature, K

600

Figure 2: Ideal gas change in free energy of the m ⇀ ↽ p (black) and m ⇀ ↽ o (red) reactions computed using tabulated thermodynamic properties (line), using the harmonic oscillator and rigid rotor approximations (filled symbols) and treating the two lowest frequency modes as hindered rotors (red x’s), and using QM/MM alchemy (open symbols). Error bars depict standard error of the mean estimated via bootstrapping. overtaking p-xylene between 600 and 700 K. Sensitivities run at 1 bar, 300-600 K, and using ∆Galch and ∆GHO/RR are shown in Fig. S2. Constant density RxMC simulations were run at a supercritical temperature, 700 K, in order to test the efficiency of the semi-grand reaction move as a function of density. A supercritical temperature was chosen to avoid phase instabilities at intermediate densities. As shown in Fig. 4 all methods achieve 90%+ acceptance at 1 kg/m3 . With increasing density, the acceptance of aligned position and orientation moves decreases approximately linearly to 26% at 700 kg/m3 , while the acceptance of random position and orientation reaction moves decays more sharply, to less than 1% at 500 kg/m3 and higher densities. Reaction moves that align position but use random orientation result in acceptance rates bracketed by the other two methods, with 3% acceptance at 700 kg/m3 .

4.3

Under Confinement

Confining a single xylene molecule to the interior of a carbon nanotube has little effect on the equilibrium composition for diameters d > 12 Å (see Fig. 5a). Below 12 Å, the fraction of m-xylene decreases monotonically to 99.9% over the same interval. The population of o-xylene remains 15



0.6

meta

Mole Fraction

0.5 0.4 para

0.3 0.2

ortho

0.1 liquid

vapor

0 300

500 700 Temperature, K

900

Figure 3: Equilibrium mole fraction of a mixture of xylene isomers as a function of temperature at a pressure of 1 bar and using ∆Gtab . Simulation averages are shown as empty circles; experimental data as filled squares with a line to guide the eye. Error bars on the former are smaller than the symbol size. The normal boiling point (412 K) is indicated by a dashed line.

1 Acceptance Ratio

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.8 0.6 0.4 0.2 0 0

100 200 300 400 500 600 700 Density, kg/m3

Figure 4: The fraction of accepted reaction moves as a function of density at 700 K. Reaction moves performed by inserting a molecule at a random location and orientation (red diamonds), with aligned positions and random orientation (green circles), and with aligned positions and orientations (blue squares). Lines are a guide to the eye.

16


Page 16 of 31

Page 17 of 31

the smallest at all diameters despite a slight increase at d = 11.7 Å before decreasing to

Reaction ensemble Monte Carlo simulation of xylene isomerization in

Recommend Documents