Elucidating the Crystal Face- and Hydration ... - ACS Publications

ARTICLE pubs.acs.org/JPCC

Elucidating the Crystal Face- and Hydration-Dependent Catalytic Activity of Hydrotalcites in Biodiesel Production Kuang Yu and J. R. Schmidt* Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin, Madison, Wisconsin 53706, United States ABSTRACT:

Hydrotalcites (HTs) are an important class of heterogeneous catalysts for a wide variety of traditionally base-catalyzed organic transformations and have been considered as an alternative to homogeneous catalysts in the transesterification of fats/oils in biodiesel production. Previous experimental observations show that HT-like materials possess a crystal face-dependent activity that varies between related organic transformations, as well as a profound influence for HT hydration on its catalytic activity. We use a combination of molecular dynamics simulations and periodic plane-wave density functional theory to elucidate both the crystal faceand hydration-dependent activity of Mg-Al HT. We develop a new force field for HT, based on the existing ClayFF force field, and utilize this model to examine the interaction of a model ester reactant with HT in the presence of differing solvents, representative of conditions for common transesterification and hydrolysis reactions. Our results suggest that the variations in the activity of the HT crystal faces can be explained in terms of dramatically different adsorption free energies of the ester onto the various faces, changing with the nature of the solvent. We further conclude that the observed hydration dependence of the HT activity is due to the deprotonation of interfacial acidic Mg(Al)-O-H hydroxyl hydrogen by interlayer hydroxide counterions in the absence of interlayer hydration, leading to a loss of Brønsted basicity.

’ INTRODUCTION Biodiesel has the potential to provide an important source of renewable carbon-neutral energy, yielding substantially more energy that that required for its production.1 Biodiesel consists of a mixture of fatty-acid alkyl esters, sharing many of the chemical and physical properties of conventional petroleumbased diesel fuel and thus maintaining compatibility with existing diesel engines. Due to its increased oxygen content and the lack of sulfur, it also exhibits a typically lower emission profile than petroleum-based diesel. Yet biodiesel is not without limitations, which include a slightly lower energy density (approximately 90% that of conventional diesel) and poor flow behavior at lower temperature.2 These limitations can be mitigated by blending biodiesel with petroleum-based diesel in blends that range from 5% (B5) to 100% (B100) biodiesel. r 2011 American Chemical Society

The principal costs associated with biodiesel production include those of the feedstock materials (typically plant-derived oils or waste animal fats) as well as the costs associated with processing the feedstock into usable fuel.1 The latter conversion is necessary because the native oils/fats consist of triglycerides with a much higher viscosity and lower volatility than petroleumbased diesel, both of which adversely effect engine performance. Conversion of triglycerides to biodiesel is conventionally accomplished by mixing the precursor oils with methanol in the presence of a homogeneous basic catalyst, resulting in a transesterification reaction yielding 3 equiv of fatty-acid methyl esters Received: June 24, 2010 Revised: October 21, 2010 Published: January 11, 2011 1887

dx.doi.org/10.1021/jp105849d | J. Phys. Chem. C 2011, 115, 1887–1898

The Journal of Physical Chemistry C

Figure 1. Schematic illustration of the basic Mg-Al hydrotalcite (HT) structure. HT consists of layers of Mg2þ (Al3þ) cations octahedrally coordinated by structural hydroxyls in the ab crystallographic plane and stacking in the c direction.

and a glycerol byproduct.3 This homogeneously catalyzed route is inefficient as it requires the use of batch reactor setups, neutralization of the caustic catalyst, and disposal of an aqueous salt waste stream.4 In contrast, a heterogeneous catalyst would enable trivial separation and reuse of the catalyst material, as well as facilitate efficient continuous-flow reaction processes. A variety of heterogeneously catalyzed pathways have been explored in the literature, involving both solid acids and solid bases, and many of these have been addressed in recent reviews.5-10 Among these methods, hydrotalcites have received considerable attention as heterogeneous catalysts for a wide variety of traditionally base-catalyzed organic transformations, including the transesterification reactions utilized in biodiesel production.11-35 Hydrotalcites are a family of structural analogues to the mineral hydrotalcite, Mg6Al2(OH)16CO3 3 4H2O, a layered double hydroxide (LDH)36 consisting of layers of divalent Mg2þ ions coordinated in an octahedral fashion by structural hydroxyl groups. These metal-hydroxyl complexes form layers via edge sharing in the ab crystallographic plane, which then stack in the c plane to form three-dimensional structures; the structure is shown schematically in Figure 1. HT is obtained by substituting a fraction of the divalent Mg2þ cations with a trivalent Al3þ, yielding layers of positive charge. This charge is balanced by the interlayer region, consisting of varying amounts of water and a charge-balancing number of anions. While the mineral hydrotalcite contains carbonate counterions, in almost all catalytic applications HT is “activated” by converting these anions to hydroxide either as a part of the synthesis, by ion exchange, or via calcination/rehydration. In the latter case, calcination results in a mixed Mg/Al oxide, and subsequent exposure to water re-forms a HT-like material with the original LDH structure where the interlayer counterions have been replaced by hydroxide. Unless otherwise noted, all subsequent references to HT refer specially to the “active” hydrotalcite-like LDH. For clarity, we will refer to all structural metal-coordinated OH groups as “hydroxyl” and the interlayer OH- anions as “hydroxide”. Abundant evidence demonstrates that HT acts as a true heterogeneous catalyst rather than simply leaching hydroxide into solution,21 although the interlayer charge-balancing hydroxide anions are nonetheless thought to play a crucial role in HT catalysis, acting as Brønsted basic sites.11,27,37 These interlayer counterions are accessible to the bulk solution only at the edge of the HT crystallites, and thus it has been suggested that these faces serve as the active sites for HT catalysis.27,31,32,34 Roeffaers et al.33 carried out in situ fluorescence microscopy experiments on a related Li-Al LDH and were able to spatially resolve the crystalface-dependent activity in both hydrolysis and transesterification

ARTICLE

reactions on a small fluorescent ester. Like its Mg-Al analogue, the Li-Al LDH has a layered structure consisting of Al3þ cations octahedrally coordinated by structural hydroxyl groups. However, it should be noted that the Li-Al LDH contains a fraction of octahedral vacancies due to the differing stoichiometry of the system. Nonetheless, the Li-Al LDH can be grown with large crystals with easily identifiable faces, a requirement for the experiment. They find that hydrolysis (reaction of the fluorescent ester with a water solvent) occurs predominantly on the edge of the crystallite, as might be expected. However, for the corresponding transesterification reaction (where the solvent/reactant is changed from water to 1-butanol) almost no selectivity between crystal faces was observed, with reactions occurring on both the basal plane (parallel to the ab crystallographic axes) and the crystal edge. This result contradicts the assumed inertness of the basal plane, which was deduced from indirect correlations between structural and reactive properties.27,31,32,34,35 Several other authors have also noted that the amount of “edge” surface area does not correlate well with observations and have suggested a possible role for the structural hydroxyl groups in catalysis via an acid-base mechanism, or that crystallite morphology may play an important role in modulating the observed activity.17,28 Additionally, Greenwell et al.35 presented evidence for the activity of surface-bound hydroxide anions from density functional theory (DFT) calculations. These studies suggest that HT is acting as more than a simple hydroxide carrier and that the metal cations, structural hydroxyl groups, and/or basal surfacebound counterions may participate directly in the catalytic events. The hydration state of HT is also an important variable in controlling its activity. Water occurs naturally in HT in at least two forms: surface-bound water molecules, physisorbed to the external face of the crystallites; and interlayer water molecules, confined between the laminar HT planes. Xi and Davis21 examined the activity of HT as a function of the interlayer water content. They found that the loss of surface water did not appreciably influence the activity of the HT but that the activity steadily decreased toward zero with the loss of interlayer water molecules. They hypothesize that this loss of activity may be due to the interaction of the interlayer hydroxide anion with the brucite-like layers as the interlayer water content decreases. In order to probe the catalytically active site/face of HT and to explain the observed crystal face- and hydration-dependent activity, we use a combination of plane-wave DFT and molecular dynamics (MD) simulations to model the structure of the HT interface under representative reactive conditions. We develop a validated force field for the relevant HT/solvent/ester interfaces and utilize this force field to model the structure of the interface under conditions representative of transesterification and hydrolysis reactions. These simulations, in conjunction with additional DFT calculations, elucidate the dominant ester/HT interaction motifs in the presence of solvent and allow us to explain the observed crystal face- and hydration-dependent activity of HT.

’ COMPUTATIONAL METHODS All plane-wave DFT calculations were carried out with version 4.6 of the Vienna ab initio simulation package (VASP),38-43 within the context of the non-spin-polarized generalized gradient approximation and the Perdew-Burke-Ernzerhof (PBE)44,45 exchange-correlation functional. The wave function for the valence electrons was expanded into a plane-wave basis set with 1888

dx.doi.org/10.1021/jp105849d |J. Phys. Chem. C 2011, 115, 1887–1898

The Journal of Physical Chemistry C an energy cutoff of 500 eV, while the interaction of the core and valence is described by the projector-augmented wave approximation. Our HT supercell and surface slab models consist of between 136 and 272 atoms, with Mg/Al in a 2:1 ratio. Given the large supercell and surface slab models employed in the DFT calculations, the Brilloiun zone was sampled only at the γ point. Molecular dynamics (MD) simulations were conducted with a modified version of GROMACS 4.0.5.46 As described in detail below, a significantly reparametrized version of the ClayFF force field47 was used to model the HT and interlayer anions. Water molecules were described by the simple point charge (SPC) water model,48 chosen due to its known compatibility with the ClayFF force field. Organic solutes and solvents were described with the OPLS/AA force field,49 as this force field has been optimized to fit many experimental properties of the liquid. Interactions between organic solutes/solvents and the HT substrate were explicitly parametrized, as described below, on the basis of DFT calculations, while other interactions were modeled by use of the standard Lorentz-Bertelot mixing rule as prescribed by the ClayFF force field. The particle mesh Ewald (PME) summation was utilized to calculate electrostatic interactions in conjunction with a 10 Å real-space cutoff. Van der Waals interactions were calculated by the force switch method to ensure that the potential and its derivative vanish smoothly between 8 and 9 Å.50 All simulations are carried out in an NPT ensemble with the exception of replica exchange and a subset of the umbrella sampling simulations, which are done in an NVT ensemble. Nose-Hoover temperature coupling in conjunction with Parrinello-Rahman pressure coupling are used to enforce the temperature and pressure constraints; unless otherwise noted, all simulations are at 300 K and 1 bar. All MD simulations are based on a 3 3 1 supercell containing three HT layers with Mg and Al atoms in a 2:1 ratio (1215 total atoms), corresponding to the most commonly employed experimental HT stoichiometry. An appropriate number of interlayer water molecules are also added to yield the desired HT hydration state. Surface models for the lowest energy HT crystal faces are generated by cleaving the supercell and then solvating with methanol, 1-butanol, or water. In all cases, the initial geometry for HT was obtained by from the DFToptimized supercell or surface model. In several cases, replica exchange molecular dynamics (REMD) was employed in order to ensure efficient and complete sampling. Briefly, REMD uses an ensemble of systems, each at a different temperature, in order to efficiently sample phase space and surmount kinetic barriers separating disparate regions of configuration space. Members of the ensemble are exchanged between temperatures in such a way that the appropriate equilibrium average at the desired (low) temperature is maintained; full details on the method can be found elsewhere.50,51 In particular, we found that REMD was necessary where an adsorbate approaches and binds to the HT surface, where the strong interactions between adsorbate and HT would otherwise result in poor sampling. In these cases, between 24 and 40 replicas were used, with the temperatures of each replica ranging between 300 and 500 K (700 K) with 24 (40) replicas; an exchange between systems at adjacent temperatures was attempted every 1 ps.

’ MODEL PARAMETRIZATION Assessing the ClayFF Force Field. We initially tested a model for HT based on the existing ClayFF force field. ClayFF

ARTICLE

was originally parametrized on the basis of a combination of DFT and experimental structural data for a diverse set of clay systems and has been widely employed for modeling bulk HT. 47 The primary advantage of ClayFF is that it contains few explicit bonded terms, treating the cation/anion interactions as purely ionic, thus allowing (in principle) significant freedom to describe not only bulk HT but also HT interfaces. However, our initial simulations of bulk HT using this force field suggest that ClayFF may contain a hidden thermodynamic instability, resulting in gradual loss of crystalline over the course of a simulation; this and other instabilities have been observed in previous studies of HT using ClayFF.52,53 Although in the case of bulk HT we did not observe any bulk “melting”, even for long simulations, cleaving the HT perpendicular to the [110] lattice vector to generate a HT surface leads rapidly to a disordered phase, suggesting that a true thermodynamic order-disorder transition is occurring at moderate to low temperatures. Note that, experimentally, HT does not actually melt but rather spontaneously dehydroxylates into an Mg/Al mixed oxide at temperatures above approximately 473 K. 21 (Since the ClayFF force field treats the structural hydroxyl group as a single unit, dehydroxylation is not possible within this simulation model.) Thus the observed “melting” is not a physically meaningful phenomenon but rather an artificial disordered phase that becomes thermodynamically favored over the (realistic) ordered phase at moderate temperatures. As this order-disorder transition occurs in a temperature regime similar to those used in transesterification or hydrolysis reactions, we require (at a minimum) a model that is thermodynamically stable at catalytically relevant temperatures of up to ∼373 K. We calculate the melting point for our HT model, choosing for simplicity the charge-balancing interlayer anions to be chloride (HT-Cl), as has been commonly studied in previous MD simulations of HT.52,54-57 We estimate the thermodynamic melting point of the ClayFF model of HT by measuring the melting velocity of an interface between a crystalline HT-Cl [110] surface and an adjacent disordered phase. The system is held at fixed temperatures varying between 360 and 480 K, and the melting velocity is quantified via an order parameter, ξ, defined by 3+ * 2 X ξ ¼ exp4i k3B r j =a5 B j

Here the summation extends over all metal ion positions, Br j , within a single HT layer, a = 3.21 Å is the norm of the B a lattice vector (using an average from our simulation), and B k is the corresponding reciprocal lattice vector. Note that it is necessary to compute this average separately for each of the three HT layers, since they slide with respect to one another throughout the course of the simulation. The average is taken over nine HT layers from three independent simulations. We define the melting velocity, v, as proportional to the initial linear decay of ξ(t), and plot the resulting temperature dependence in Figure 2. The thermodynamic melting point is defined to be the temperature where the melting velocity goes to zero. (Since the crystallization of HT is extremely slow, even in the presence of an existing interface, it is not possible to make use of 1889



ARTICLE

Table 1. Comparison of Original47 and Modified ClayFF Force Field Parameters and Associated Properties for HTa original ClayFF

modified ClayFF

DFT

Charges (e) H

0.425

0.464

O Al

-0.871 1.575

-1.113 2.018

1.050

1.438

Mg

k (kJ 3 mol-1 3 rad-2) Mg-O-H

250

53 -1

-2

[001]

3.81

3Å ) 3.84

[110]

1.56

3.77

4.24

[100]

3.38

6.26

6.35

Surface Energy (kJ 3 mol

Figure 2. Calculated melting velocity, v, from the unmodified ClayFF force field. Simulation data () are shown alongside the physically motivated fit (solid line); the melting point is determined by extrapolating the melting velocity to zero.

-7

Diffusivity (10

2

2.40

-1

cm 3 s ) 6.87 ( 0.4 0.21 ( 0.02

crystallization data.) We fit the melting curve to the physically motivated functional form58 - ðT - Tm ÞΔS v∼exp - EA =kb T 1 - exp kb T

a

where the liquid (disordered) phase diffusion activation energy is EA, Tm is the thermodynamic melting point, and ΔS is the corresponding entropy change upon melting. EA is obtained from the temperature dependence of the diffusion coefficient in the disordered phase (measurements for the various different ions yield nearly identical results), while the other parameters are obtained from the fit. The extrapolated melting temperature of HT-Cl is 360 K, far below the experimentally observed decomposition temperature of HT and within the range of catalytically relevant temperatures. The artificially low melting point suggests that the cohesive energy predicted by ClayFF is too small. We investigate this hypothesis by calculating the unrelaxed surface energy for HT-Cl along several of the known lowest energy lattice planes.59 We cleave the HT-Cl crystal (using the DFT-optimized bulk structure) normal to the given lattice vector and introduce a 10 Å vacuum gap between the two crystal faces thus generated. We then calculate the unrelaxed surface energy as the difference between the bulk and cleaved energies, using both DFT and the ClayFF force field. The results are summarized in Table 1. Note that the ClayFF-predicted [110] and [100] surface energies are far too low, consistent with the fact that ClayFF underestimates the cohesive energy. We presume that since the vast majority of previous MD simulations examined only bulk HT or the HT basal surface (which does not require breaking any covalent bonds or disrupting the metal coordination sphere), the thermodynamic instability that we highlight here was likely not observed due to the kinetic stability of bulk HT. However since the “edges” of HT have been implicated in its catalytic activity, we require a new HT simulation model that is capable of describing both the HT basal plane and the orthogonal edges. Reparametrizing ClayFF. To resolve the instabilities with the existing ClayFF force field, we reparametrize the force field specifically for HT. We obtain a new point charge model for HT based on fitting the electrostatic potential (ESP) obtained from DFT calculations of several HT interfaces. While this procedure is well-established for molecular systems, the present case is more

complicated since the absolute ESP of a bulk solid is undefined (since for an infinite periodic solid the ESP will not vanish at infinity). Although one possible solution to this dilemma was recently proposed,60 we instead utilize plane-wave DFT calculations of the three lowest energy surfaces of HT-Cl ([001], [110], and [100]),59 which would be expected to dominate the surface area of crystalline HT under equilibrium conditions. In each case, we cleave the HT-Cl crystal orthogonal to the given surface normal and calculate the ESP within the vacuum gap, assigning a value of zero to the plateau region that occurs between the two surfaces. The ESP is obtained directly from VASP on a dense grid of points, which is consistent with the potential arising from the nuclei and the diffuse electron density of the system. We use an analogous procedure to calculate the ESP arising from an empirical point-charge model, in this this case calculating the ESP via the classical Ewald summation. Sampling the ESP on planes between 1.5 and 2.1 Å from the three interfaces (the distance is chosen on the basis of energetic criteria and to ensure similar variations in the ESP among the three surfaces), we find the point-charge model that best reproduces the calculated ESP. We constrain the charge of the Cl- anion to be -1, and the remaining Mg, Al, O, and H charges are optimized to best reproduce the DFT ESP at a variety of symmetry-unique points along the planes. We subsequently verified that the ESP is also well reproduced over a wide variety of distances out of plane. The DFT and fitted ESP are shown for the two thermodynamically most favorable interfaces in Figure 3, and the final fitted charges are given in Table 1. The ESP for the [100] surface is somewhat less accurate due to the presence of exposed Cl- anions and the associated constraint on the Cl charge in the fitting procedure. However, the surface energy of this interface is considerably higher than the others and is unlikely to make a significant contribution to the surface area of HT crystallite; we do not consider it further. It is also interesting to note that even the original ClayFF force field accurately reproduced the ESP in the [001] direction (not shown), motivating the success of previous simulations using ClayFF

H2O Cl-

6.30 ( 0.6 0.36 ( 0.12

DFT calculated values are also given for comparison, where appropriate.

1890



ARTICLE

Figure 3. Comparison of DFT- and modified ClayFF-calculated ESP for (A, B) [001] and (C, D) [110] surfaces. Positions are given in direct lattice coordinates for the supercell, and all ESP values are given in volts.

Figure 4. Comparison of surface phonon spectra for the [110] interface as calculated by (A) DFT, (B) ClayFF force field, and (C) our modified ClayFF force field, including Mg-O-H bond bending terms.

for examining bulk HT or the HT basal surface. Yet when the original charges are compared to the new ones, the new point charges are far more ionic, suggesting that interactions with other HT surfaces may differ considerably between the original and revised force field. Applying the refitted charges in conjunction with the original Lennard-Jones (LJ) parameters, we calculate a lattice constant of 3.07 Å, in excellent agreement with the experimental lattice constant (3.05 Å)61 and an improvement from the original ClayFF (3.20 Å). Thus we do not make further adjustments to the LJ parameters. We also recalculate the HT melting point, which is now 690 K and thus higher than the observed decomposition temperature, ensuring the stability of our model at all relevant temperatures. Interestingly, we note that the calculated diffusion constants for the interlayer water and Cl- counterions do not change substantially from the original force field (see Table 1). This is almost certainly because both the original and reparametrized force fields yield similar ESP values in this region, both in agreement with that calculated via DFT. Yet we would expect dramatic differences in properties involving HT edge interfaces. The resulting MD model yields an excellent description of the bulk structural properties of HT. However, the structure of the

Figure 5. (A) Calculated HT swelling, (B) hydration energy, and (C) interlayer water excess chemical potential as a function of interlayer water content, given as a percent by mass. Simulation data are given as points, while the dashed line gives (A) experimental HT basal spacing, (B) heat of vaporization, and (C) chemical potential of bulk SPC water.

[110] interface is more disordered than predicted on the basis of DFT calculations, the former predicting significant structural disorder of the interfacial structural hydroxyl anions at room temperature. We quantify the interfacial structure and dynamics by calculating a surface-projected phonon spectrum, shown in Figure 4; for simplicity, we evaluate this term using pure brucite, Mg(OH)2, whose structure is identical to that of HT. While the surface phonon spectrum from the MD model is denser due to the larger supercell employed, the primary difference when 1891



ARTICLE

compared to the DFT model is the lack of surface phonons in the regions from 500 to 1000 cm-1, corresponding to in-plane torsions of interface structural hydroxide groups. We remedy this deficiency by adding a harmonic Mg-O-H (Al-O-H) angle bending term of the form Ebend = 1/2k(θ - θ0)2 where θ0 is chosen to be the DFT equilibrium bond angle in the interfacial region, 132. The force constant k is roughly evaluated by comparing the surface phonon spectrum and is given in Table 1; the resulting phonon spectrum is in excellent qualitative agreement with that obtained via DFT. Note the original ClayFF force field also contains an (optional, almost universally neglected) angle bending term, although the force constant is much larger than present. Finally, we require a model for the interlayer hydroxide anions, for which the Cl- counterions are exchanged in “activated” HT. We employ a hydroxide model based on SPC water where one of the hydrogen atoms is removed, and the associated charge is added to the O atom,62 yielding charges quite consistent with those calculated by fitting to the ESP. Since the prior hydroxide model was parametrized only for aqueous solution, we also adjusted slightly the effective LJ diameter of the hydroxide oxygen to σ = 3.535 Å for all hydroxide/nonwater interactions based on DFT calculations of the hydroxide/HT interaction energy.

’ RESULTS AND DISCUSSION Equilibrium Hydrotalcite Hydration. In the majority of cases HT is activated by calcination and subsequent rehydration by exposure to liquid water. In this case, the interlayer water will exist in equilibrium with the bulk. While this equilibrium water content can be estimated on the basis of the weight loss via thermogravimetric analysis (TGA),21,23,31 the relatively low decomposition temperature of HT makes it difficult to tell whether all or only a fraction of the interlayer water molecules have been driven off prior to decomposition. In addition the peaks in the differential TGA signal corresponding to the loss of physisorbed surface water and the loss of interlayer water overlap, making a definitive estimate of the interlayer water content hard to reach.21 Since it has been demonstrated that the HT hydration state plays a crucial role in modulating its activity,21 we will first utilize the reparametrized ClayFF force field to study the hydration of activated HT and to determine the water content of HT under conditions of equilibrium with bulk water. As HT is hydrated, water molecules progressively occupy the interlayer region and the clay material swells as the distance between successive HT layers increases. We use our simulation model to examine the swelling with respect to the HT water content, as shown in Figure 5. The basal spacing between two adjacent HT layers shows three-stage swelling. Initially, the interlayer molecules form a two-dimensional hydrogen-bonding network incorporating the interlayer hydroxide anions, where both the water and counterions lie in plane; a representative example of this structure is shown in Figure 6. This result lies in contrast to previous MD simulations of HT-Cl, which showed distinctly out-of-plane hydrogen bonding, with water coordinating the interlayer Cl- from both above and below.63 Here the smaller hydroxide anion results in decreased interlayer distance and thus a more confined space that does not admit out-of-plane hydrogen bonding. At higher water content the interlayer structure transitions to a double-layer hydration structure and the basal spacing rapidly increases to accommodate the

Figure 6. Depiction of the equilibrium HT hydration state (top), showing a two-dimensional interlayer hydration structure (bottom). In the latter case, the HT layers have been removed to clearly show the 2D hydrogen-bonding arrangement.

additional water. This double layer is formed initially only in a fraction of the supercell, resulting in buckling of the basal surface. Continued hydration results in a complete double layer and reforms a flat structure. We validate our MD simulations against DFT calculations, optimizing the interlayer spacing for a small HT/water cell in several hydration states. Although the DFT-optimized systems utilize smaller unit cells and are effectively at 0 K, they yield qualitatively identical hydrogen-bonding arrangements and reproduce the planarity of the interlayer region. However the DFT results give systematically slightly larger basal spacings; this is due, in part, to the fact that the buckling deformations observed in the MD simulation are not possible in this smaller unit cell, forcing instead an overall increase in the interlayer distance. Using the experimentally measured HT interlayer spacing of 7.7 Å for a sample of rehydrated HT and comparing with our calculated swelling curve,27 we would estimate a HT interlayer 1892


The Journal of Physical Chemistry C hydration of 18% by weight. However, this estimate contains large uncertainties due to the weak dependence of the HT interlayer distance, making this estimate extremely sensitive to small errors in the predicted MD interlayer spacing. As an alternative estimate of the hydration, we calculate the incremental hydration energy of HT by hEðNi þ 1 Þi - hEðNi Þi Ehyd Ni þ 1=2 ¼ ΔN where ÆE(Ni)æ is the average potential energy from a simulation containing Ni water molecules and ΔN = Niþ1 - Ni is the change in the number of water molecules between the two different hydration states. The resulting plot is shown in Figure 5 and also provides evidence for three-stage hydration. Initially the HT shows strong affinity for water due to the strong interactions with the charged HT layers. Additional water can be accommodated only by swelling and/or buckling the HT layers, which is energetically unfavorable since mutual layer/counterion/layer electrostatic attractions must be overcome.54 After the layers have swelled sufficiently, additional waters can be accommodated without swelling, leading once again to very energetically favorable hydration. Neglecting entropic effects, we can estimate the extent of HT hydration by equating this hydration energy to the heat of vaporization of bulk SPC water, yielding an estimate of approximately 16% by mass. However, due to the highly confined nature of the interlayer region, it is likely that the entropy of the interlayer water is significantly lower than that of bulk water, and thus this estimate is an upper bound on the actual hydration. We account for entropic effects by calculating the excess chemical potential of the HT interlayer water as a function of hydration. We use the Bennett estimator for the excess chemical potential64,65

! F - β ΔE - μ0 d βμex ¼ ln þ βμ0ex F β ΔE - μ0 i FðxÞ ¼ 1=½1 þ expðxÞ where ΔE is the potential energy difference upon particle insertion or deletion, μ0ex is an initial estimate for the excess chemical potential, and the subscripts d and i denote averages over trial particle deletions and insertions, respectively. The Bennett method is a more sophisticated implementation of the standard Widom insertion method,66 using a combination of trial insertions and deletions, along an initial estimate of the chemical potential, to reduce the variance of the estimator and achieve far more accurate results in the case of a dense liquid. The procedure is iterated to self-consistency between μ0ex and μex, with the resulting excess chemical potential shown in Figure 5. Equating the excess chemical potential of the interlayer water to that of bulk SPC water suggests that HT at equilibrium with the bulk contains approximately 14% water by mass, slightly lower than would be predicted on the basis of the hydration energy due to the latter’s neglect of entropic effects. Experimental TGA data yields a slightly higher estimate for the interlayer water content, approximately 14-16% by weight; however, the overlap of peaks for surface and interlayer waters makes it hard to distinguish between these two contributions (see, for example, Figures 1 and 2 of ref 21). Entropic effects can be estimated (neglecting PV terms) from the calculated hydration energy and chemical potential, along with the corresponding values from

ARTICLE

bulk SPC water. On the basis of these values, we estimate that the partial molar entropy of HT (between 11 and 15% water by weight) is approximately 150 J/(mol 3 K) lower than that of bulk SPC water. Thus the partial molar entropy is much lower in the interlayer region, leading to an entropic penalty for HT hydration. Such a large difference is not surprising as the interlayer region is strongly confined, leading to loss of significant of translational and rotational freedom. As discussed above, the structure of the interlayer water at equilibrium hydration is unique, consisting of a highly confined 2D hydrogen-bonding network incorporating both water molecules and hydroxide anions. This arrangement is facilitated due to the strong “external pressure” exerted by the mutual attraction of the adjoining HT layers, mediated by the hydroxide anion. For convenience we use HT-OH 3 2H2O (two interlayer water molecules per hydroxide anion), with the same equilibrium structure, as our model system for all subsequent simulations. Role of Hydration in Hydrotalcite Reactivity. Recently Xi and Davis21 found that the catalytic activity of HT was influenced significantly by the extent of interlayer hydration (which they controlled by progressively heating the HT to increasingly higher temperatures), with a dramatic loss in activity as interlayer water was removed. In order to examine this phenomenon, we created models of the HT [110] interface at two different hydration states corresponding to (i) a complete lack of interlayer water and (ii) the equilibrium hydration state determined above. We then relaxed and optimized these two interfacial structures by DFT. Examining the dehydrated structure, we find that all of the interlayer hydroxide anions in the interfacial region deprotonate their neighboring structural hydroxyl groups. This striking result appears to be a strictly interfacial phenomenon, as corresponding DFT calculations on bulk dehydrated HT were stable, exhibiting no deprotonation. This suggest that the interfacial structural hydroxyl groups are more acidic than their bulk counterparts, most likely because the undercoordinated Mg2þ/Al3þ cations are better able to stabilize the resulting conjugate base anion. In contrast, the hydrated HT sample was stable with respect to spontaneous deprotonation during the course of the optimization. Here the presence of interlayer water molecules solvates and stabilizes the interlayer hydroxide anions, decreasing their basicity and preventing the deprotonation of the structural hydroxyl groups. Deprotonation of the structural hydroxyl groups by the interlayer hydroxide anions corresponds to a significant loss of Brønsted basicity for the catalyst. Inasmuch as these anions are intimately involved in the HT catalysis, their neutralization would lead to a significant loss in the activity of HT. As this acid-base neutralization is hindered by the presence of interlayer solvent, we would thus suggest that the decrease in HT activity observed with the loss of interlayer solvent may be related to loss of solvent stabilization, resulting in interfacial hydroxide neutralization and leading to a reduction of the available Brønsted basicity. This conclusion supports the general hypothesis of Xi and Davis21 that the decreased activity arises from the direct interaction of the hydroxide anion with the brucite-like layers, although the present work yields a much more specific picture of the nature of those interactions. Structure of Hydrotalcite-Solvent Interface. The vast majority of HT-catalyzed transesterification reactions for biodiesel production utilize solutions of methanol and triglyceride (esters), producing 3 equiv of methyl esters (biodiesel). We thus 1893



ARTICLE

Table 2. Comparison of DFT- and Force Field-Calculated Adsorption Energies for Methanol and Methyl Acetatea force field

DFT

Methanol [001]

-104.8 (-104.8)

-104.6

[110]

-63.9 (-63.9)

-64.5

Methyl Acetate [001]

-70.6 (-59.9)

-51.9

[110]

-51.6 (-90.3)

-83.8

a

All values are given in kilojoules per mole. Values in parentheses are for the optimized force field, after reparametrization.

Table 3. Summary of Optimized Force Field Parameters Governing the Interaction of HT Interfaces with Methyl Acetatea σAB (Å)

A

B

carbonyl O

hydroxyl O

carbonyl O

εAB (kJ/mol)

[001] Surface 3.551

0.2325

hydroxyl O

2.590

0.7558

carbonyl O

hydroxyl H

1.949

0.8956

carbonyl O

Al

4.416

0.002 211

[110] Surface

carbonyl O chargeb a

-0.56

These parameters are used only for the specified interactions. b This charge is used only to calculate HT-ester interaction. The compensating charge is distributed equally among the remaining atoms in the molecule.

study the structure of the two lowest energy HT-methanol interfaces in the presence of two different model esters to explain the role of HT-ester and HT-solvent interactions in determining the spatial selectivity of HT catalysis. We utilize the OPLS/AA force field to describe the methanol solvent, as well as the methyl acetate and methyl butyrate, our model esters.49 Comparison across the two esters allows for an assessment of the role of the alkyl chain length in the adsorption process; furthermore, methyl butyrate serves as a proxy for the (even larger) glyceryl tributyrate often utilized as an experimental model system in catalysis studies.9,67 Since OPLS and ClayFF were parametrized by different methodologies, we use DFT calculations to benchmark the combined force field by calculating the binding energy of both methanol and methyl acetate on the HT [001] and [110] surfaces and comparing to the corresponding calculation (at the same DFT geometry) using our empirical force field; results are shown in Table 2. It is clear that the results for methanol are excellent, while significant errors exist for methyl acetate. We compensate by adjusting the cross terms between the two force fields governing the interaction between the ester and HT. We select one binding site for the ester on the [001] surface (the carbonyl oxygen binding with the structural hydroxyl hydrogen), and three binding sites on the [110] surface (the carbonyl oxygen binding with the Mg, Al, and structural hydroxyl hydrogen). These structures are optimized at the DFT level, and the molecule is scanned perpendicular to the surface. These energies are used to optimize the force field cross terms, whose optimized values are listed in Table 3. While some residual errors remain, the optimized adsorption

Figure 7. Comparison of the number density of methanol oxygen atoms in the vicinity of Mg2þ (dashed red line) and Al3þ (solid black line) cations.

energies are dramatically improved and qualitatively correct. Given the similar nature of the interactions, the force field governing the methyl butyrate interactions with HT was derived from that of the corresponding methyl acetate interactions, using the prescribed OPLS interaction parameters and charges for the additional akyl chain atoms; the small residual charge (∼0.1e) resulting from this procedure was distributed equally among the approximately 30 atoms in the ester to enforce charge neutrality, resulting in a negligible change. Using our validated force field, we examine the structure of the HT-methanol [001] and [110] interfaces. The interface is created by cleaving the bulk HT along the plane corresponding , which is to the desired surface and creating a gap of at least 20 Å then filled with an appropriate number of solvent molecules consistent with the bulk density. We pin several Mg2þ cations far from the interface of one of the HT layers to prevent translation of the bulk HT; however, since only one interface is pinned, the layers are still free to slide with respect to one another (although for the low equilibrium interlayer water concentrations studied here, no sliding was observed). The [001] interface is wellordered, with all hydroxide counterions strongly hydrogenbonded to the structural hydroxyl groups and oriented largely perpendicular to the HT surface. The dynamics of the surfacebound counterions consist of primarily activated hopping between various surface binding sites. Our simulation model also contains adsorbed surface water molecules (whose presence under experimental conditions will depend sensitively on the details of the catalyst preparation) that bridge two counterions via hydrogen bonding. The interfacial methanol molecules interact with the surface primarily by hydrogen-bonding to the structural hydoxyl groups, orienting the hydrophobic methyl group toward the bulk liquid phase and forming a well-organized solvation layer near the surface. The structure of the [110] interface is more complex and exposes many coordinatively unsaturated metal cations and oxygen atoms, giving rise to strong interactions with solvent. These unsaturated cations originate due to the formation of the interface; the bridging hydroxyl groups that were formerly charged between cations across the cleavage plane must now associate with one side of the interface or the other. For symmetry and to minimize the dipole moment of the system, we have assigned half of the hydroxyl groups to each side of the 1894


The Journal of Physical Chemistry C interface. Due to the natural disorder associated with this particular HT facet, this interface does not contain a wellorganized solvation shell. While some of the methanol donates hydrogen bonds to the structural hydroxyl groups at the edge of the HT crystallite, the majority of the methanol molecules interact via a Lewis acid-base mechanism involving their oxygen atom and an exposed Mg2þ or Al3þ cation. Surprisingly, the large majority of the methanol molecules interact with the Mg2þ cations rather than the more Lewis-acidic Al3þ; Figure 7 compares the number density of methanol oxygen atoms in the vicinity of Mg2þ/Al3þ cations. This phenomenon arises due to the surface reconstruction of the HT, whereby the Al3þ cations tend to bury themselves slightly into the interior of the crystal and become inaccessible to the surrounding solvent. This phenomenon appears to be significant in altering the surface properties and reactivity of HT and is reproduced in the corresponding DFT optimizations of the HT interface in vacuo. We further examine the interaction of the model esters with the [001] and [110] HT surfaces. As shown in Figure 8, the predominant interaction motif on the [001] surface involves the donation of a surface hydroxyl hydrogen bond to the ester carbonyl oxygen, with the carbonyl group nearly always oriented perpendicular to the HT surface. Along the [110] interface the ester interacts primarily with the exposed Mg2þ cations, forming a strong Lewis acid-base adduct with the ester carbonyl oxygen. The resulting adduct is extremely strong, remaining stable over several nanoseconds of simulation in the absence of enhanced sampling methods. We calculate the adsorption free energy of the model ester onto each of the HT surfaces, using umbrella sampling in conjunction with REMD to enhance the sampling of binding configurations as the ester approaches the surface.50 Briefly, umbrella sampling introduces a biasing potential to steer the system toward a particular region of phase space. By adjusting the location of the position of the biasing potential, a particular coordinate of interest can be systematically explored. As the magnitude of the biasing potential is known, its effect can be compensated to calculate equilibrium free energy differences via a potential of mean force (PMF). We define a one-dimensional reaction coordinate for adsorption via the normal distance between the ester carbonyl oxygen and the interface, given by the average position of the interfacial metal cations. We combine the sampling data from the various umbrella regions using the weighted histogram analysis method, and the resulting PMFs are shown in Figure 9. Adsorption on the [001] surface is unfavorable with respect to bulk solution, as adsorption of the ester necessarily displaces adsorbed methanol solvent, which adsorbs more strongly on this surface. Nonetheless, the PMF clearly displays two local minima. The first is located approximately 4 Å from the interface and corresponds to a true adsorbed species stabilized due to the interaction of the carbonyl oxygen and the surface hydroxyl hydrogen. The second minimum occurs further from the surface, at approximately 6.5 Å, and is due to interaction of the ester with the adsorbed first solvation layer. In this case, the ester interacts with the adsorbed methanol and/or surface-bound counterions, all oriented so as to create a surface dipole that stabilizes the adsorbed ester. In contrast, adsorption on the [110] surface is strongly favorable, once again consistent with the trend in the gas-phase DFT adsorption energy. Here, as well, the PMF contains two local minima. The first minimum corresponds to the strongly surface-bound Lewis acid-base adduct, while the

ARTICLE

Figure 8. Representative illustration of the interaction of a model ester, methyl butyrate, with (top) [001] surface and (bottom) [110] surface. The HT is shown as balls-and-sticks, the methanol solvent as lines, and the ester as a space-filling model. Consistent with the periodic boundary conditions, a second HT interface occurs opposite the one shown here.

second minimum is due to hydrogen bonding between the carbonyl oxygen and the structural HT hydroxyl hydrogen atoms. Comparing the adsorption of methyl acetate and methyl butyrate, we find that the resulting PMFs are qualitatively similar but with significant shifts in the absolute free energy as the esters approach the HT surface. In particular, the free energy at the (local) minimum nearest the surface, corresponding to the adsorbed ester, is almost 10 kJ/mol higher in energy for methyl butyrate than for methyl acetate. We attribute this decrease in stability to the influence of the long alkyl tail of the former, which has both energetic and entropic consequences. As the ester and associated alkyl tail approach the surface, its steric bulk displaces adsorbed surface solvent, resulting in an energetic penalty (as the alkyl tail cannot interact with the unsaturated cations via strong Lewis acid/base interactions). Similarly, the presence of the HT interface results in a loss of configurational entropy for longchain tails. As the alkyl tail is further increased in length to the C16-C20 esters commonly utilized in biodiesel production, we expect these trends to continue, and thus it is likely that ester adsorption, at least on the [001] basal face and possibly even on the [110] edge, is a highly activated process Solvent Effects in Hydrotalcite Spatial Selectivity. To connect with the recent in situ fluorescence microscopy experiments 1895



ARTICLE

Figure 9. (A, B) Comparison of the PMF for methyl acetate adsorption onto (A) [001] and (B) [110] surfaces in water (solid black line), methanol (dashed red line), and 1-butanol (dot-dashed blue line) solvents. (C, D) Adsorption of methyl acetate (black) and methyl butyrate (red) in methanol solvent on (C) HT [001] and (D) [110] surfaces.

of Roeffaers et al.,33 who observed different crystal face selectivity of HT depending on the solvent used, we also model and contrast the HT-butanol and HT-water interfaces. In the case of the HT-butanol interface, we use an OPLS/AA model for 1-butanol, with parameters modified in a manner consistent with methanol. For the HT-water interface, we utilize the SPC water model to maintain consistency with the treatment of the interlayer water; note that in this case we have not specifically parametrized the HT-water interactions at the interface, and thus these results are likely not quantitatively accurate. However, preliminary DFT calculations on HT surface models verify that the model captures the qualitative trends in the HT-alcohol and HT-water interactions. (A reparametrization akin to that of methanol is not possible in the present case, since the same SPC water model describes both the interlayer and bulk water models.) The corresponding adsorption PMFs are shown in Figure 9. The substantial variations in the ester adsorption free energy with solvent present one possible explanation for the spatial selectivity of HT. Since HT acts as a true heterogeneous catalyst, the first step in the reaction is almost certainly the adsorption of the ester reactant onto the surface of the HT. Our simulations show that, for all solvents considered, adsorption on the [001] surface is endergonic (increasing free energy) and must overcome a substantial barrier, while adsorption on the [110] surface is exergonic (decreasing free energy) and barrierless. Assuming both surfaces are potentially active, one would a priori expect the [110] surface to exhibit greater reactivity, as the adsorption free energies differ by 30-40 kJ/mol. Furthermore, adsorption on the [001] surface is approximately 10 kJ/mol less endergonic in aqueous solution (hydration reactions) than for 1-butanol (transesterification reactions), suggesting that the reactivity on this surface would be greater in 1-butanol than in water. In the absence of a full kinetic analysis, including all steps (adsorption, reaction, and desorption), it is difficult to discern the impact of

the adsorption step on the overall rate. However, both observations are largely consistent with the experimental results of Roeffaers et al.,33 who found that HT-catalyzed ester hydrolysis took place predominantly along the [110] surface, while ester transesterification in 1-butanol displayed little selectivity between the [001] and [110] surfaces. In any case, the profound differences in the adsorption free energy make a compelling argument that this selective adsorption may play an important role in explaining the experimentally observed crystal face selectivity.

’ CONCLUSIONS Using a combination of molecular simulations and density functional calculations, we propose explanations for the hydration- and crystal face-dependent activity of hydrotalcite (HT). In the former case, we find that HT which lacks interlayer water molecules undergoes a loss of Brønsted basicity due to neutralization of the interlayer hydroxide anions by acidic interfacial structural hydroxide groups. We further find that this is a purely interfacial phenomenon and that structural hydroxyl groups in the bulk interior are not substantially acidic and thus not deprotonated. However, for hydrated HT, the presence of interlayer water molecules stabilizes the interlayer hydroxide anion via hydrogen bonds with the surrounding solvent. In this case the interlayer hydroxides do not deprotonate the structural hydroxyl groups, preserving these active species. This picture provides a molecular-level explanation for the observations of Xi and Davis21 that the activity of HT steadily decreases with the loss of interlayer water. We also propose that the crystal face-dependent activity of HT in transesterification and hydrolysis reactions can be explained in terms of the adsorption free energies of the ester onto the various crystal faces, modulated by the effect of solvent. Consistent with in situ fluorescence microscopy experiments of Roeffaers et al.,33 we find that adsorption on the [110] surface is much 1896


The Journal of Physical Chemistry C (approximately 30 kJ/mol) more favorable as compared to the [001] surface. Inasmuch as the adsorption step influences the overall catalyst turnover rate, one would expect a much higher rate on the former surface. The large adsorption energy on the [110] surface is attributed to the formation of a strong Lewis acid-base adduct involving the ester carbonyl oxygen and interfacial coordinatively unsaturated Mg2þ cations. As the solvent is changed from water to methanol to 1-butanol, ester adsorption on the [001] surface becomes relatively less endergonic, consistent with the experimental observation that reaction on the [001] surface is much more prevalent in transesterification reactions (1-butanol solvent) than in hydrolysis reactions (aqueous solvent). Although we have not examined the molecular details of subsequent steps in the HT-catalyzed transesterification reaction, and thus cannot quantitate the influence of the adsorption rate in the overall reaction, we will examine these steps in a future publication.

’ AUTHOR INFORMATION Corresponding Author

*E-mail [email protected].

’ ACKNOWLEDGMENT We are grateful for support and funding from the University of Wisconsin via the Wisconsin Alumni Research Foundation. We also thank Arun Yethiraj and Juan de Pablo for helpful discussion. ’ REFERENCES (1) Di Serio, M.; Tesser, R.; Pengmei, L.; Santacesaria, E. Energy Fuels 2007, 22, 207. (2) Lotero, E.; Liu, Y.; Lopez, D. E.; Suwannakarn, K.; Bruce, D. A.; Goodwin, J. G. Ind. Eng. Chem. Res. 2005, 44, 5353. (3) Demirbas, A. Energy Policy 2007, 35, 4661. (4) Di Serio, M.; Ledda, M.; Cozzolino, M.; Minutillo, G.; Tesser, R.; Santacesaria, E. Ind. Eng. Chem. Res. 2006, 45, 3009. (5) Smith, B.; Greenwell, H. C.; Whiting, A. Energy Environ. Sci. 2009, 2, 262. (6) Sivasamy, A.; Cheah, K. Y.; Fornasiero, P.; Kemausuor, F.; Zinoviev, S.; Miertus, S. ChemSusChem 2009, 2, 278. (7) Melero, J. A.; Iglesias, J.; Morales, G. Green Chem. 2009, 11, 1285. (8) Lestari, S.; M€aki-Arvela, P.; Beltramini, J.; Lu, G.; Murzin, D. ChemSusChem 2009, 2, 1109. (9) Jothiramalingam, R.; Wang, M. K. Ind. Eng. Chem. Res. 2009, 48, 6162. (10) Marchetti, J. M.; Miguel, V. U.; Errazu, A. F. Renewable Sustainable Energy Rev. 2007, 11, 1300. (11) Constantino, V. R. L.; Pinnavaia, T. J. Inorg. Chem. 1995, 34, 883. (12) Choudary, B. M.; Kantam, M. L.; Reddy, C. V.; Aranganathan, S.; Santhi, P. L.; Figueras, F. J. Mol. Catal. A: Chem. 2000, 159, 411. (13) Corma, A.; Hamid, S. B. A.; Iborra, S.; Velty, A. J. Catal. 2005, 234, 340. (14) Climent, M. J.; Corma, A.; Hamid, S. B. A.; Iborra, S.; Mifsud, M. Green Chem. 2006, 8, 524. (15) Figueras, F.; Kantam, M. L.; Choudary, B. M. Curr. Org. Chem. 2006, 10, 1627. (16) Bulbule, V. J.; Borate, H. B.; Munot, Y. S.; Deshpande, V. H.; Sawargave, S. P.; Gaikwad, A. G. J. Mol. Catal A: Chem. 2007, 276, 158. (17) Lei, X.; Zhang, F.; Yang, L.; Guo, X.; Tian, Y.; Fu, S.; Li, F.; Evans, D. G.; Duan, X. AIChE J. 2007, 53, 932. (18) Liu, Y.; Lotero, E.; Goodwin, J. G.; Mo, X. Appl. Catal., A 2007, 331, 138.

ARTICLE

(19) Gao, L. J.; Xu, B.; Xiao, G. M.; Lv, J. H. Energy Fuels 2008, 22, 3531. (20) Macala, G. S.; Robertson, A. W.; Johnson, C. L.; Day, Z. B.; Lewis, R. S.; White, M. G.; Iretskii, A. V.; Ford, P. C. Catal. Lett. 2008, 122, 205. (21) Xi, Y.; Davis, R. J. J. Catal. 2008, 254, 190. (22) Chuayplod, P.; Trakarnpruk, W. Ind. Eng. Chem. Res. 2009, 48, 4177. (23) Xi, Y. Z.; Davis, R. J. J. Catal. 2009, 268, 307. (24) Climent, M. J.; Corma, A.; De Frutos, P.; Iborra, S.; Noy, M.; Velty, A.; Concepcion, P. J. Catal. 2010, 269, 140. (25) Gao, L. J.; Teng, G. Y.; Lv, J. H.; Xiao, G. M. Energy Fuels 2010, 24, 646. (26) Liu, H. H.; Xu, W. J.; Liu, X. H.; Guo, Y.; Guo, Y. L.; Lu, G. Z.; Wang, Y. Q. Kinet. Catal. 2010, 51, 75. (27) Winter, F.; Xia, X. Y.; Hereijgers, B. P. C.; Bitter, J. H.; van Dillen, A. J.; Muhler, M.; de Jong, K. P. J. Phys. Chem. B 2006, 110, 9211. (28) Greenwell, H. C.; Holliman, P. J.; Jones, W.; Velasco, B. V. Catal. Today 2006, 114, 397. (29) Chimentao, R. J.; Abello, S.; Medina, F.; Llorca, J.; Sueiras, J. E.; Cesteros, Y.; Salagre, P. J. Catal. 2007, 252, 249. (30) Debecker, D. P.; Gaigneaux, E. M.; Busca, G. Chemistry-a European Journal 2009, 15, 3920. (31) Abello, S.; Medina, F.; Tichit, D.; Perez-Ramirez, J.; Groen, J. C.; Sueiras, J. E.; Salagre, P.; Cesteros, Y. Chem.;Eur. J. 2005, 11, 728. (32) Roelofs, J.; van Dillen, A. J.; de Jong, K. P. Catal. Today 2000, 60, 297. (33) Roeffaers, M. B. J.; Sels, B. F.; Uji-i, H.; De Schryver, F. C.; Jacobs, P. A.; De Vos, D. E.; Hofkens, J. Nature 2006, 439, 572. (34) Roelofs, J.; Lensveld, D. J.; van Dillen, A. J.; de Jong, K. P. J. Catal. 2001, 203, 184. (35) Greenwell, H. C.; Stackhouse, S.; Coveney, P. V.; Jones, W. J. Phys. Chem. B 2003, 107, 3476. (36) Evans, D. G.; Slade, R. C. T. Structural aspects of layered double hydroxides. In Layered Double Hydroxides; Springer-Verlag: Berlin, 2006; Vol. 119, p 1. (37) Trave, A.; Selloni, A.; Goursot, A.; Tichit, D.; Weber, J. J. Phys. Chem. B 2002, 106, 12291. (38) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 59, 1758. (39) Kresse, G.; Hafner, J. Phys. Rev. B 1994, 49, 14251. (40) Kresse, G.; Hafner, J. Phys. Rev. B 1993, 47, 558. (41) Kresse, G.; Furthm€uller, J. Phys. Rev. B 1996, 54, 11169. (42) Kresse, G.; Furthm€uller, J. Comput. Mater. Sci. 1996, 6, 15. (43) Bl€ochl, P. E. Phys. Rev. B 1994, 50, 17953. (44) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865. (45) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1997, 78, 1396. (46) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J. Chem. Theory Comput. 2008, 4, 435. (47) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255. (48) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981. (49) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225. (50) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, CA, 2002. (51) Hukushima, K.; Nemoto, K. J. Phys. Soc. Jpn. 1996, 65, 1604. (52) Thyveetil, M. A.; Coveney, P. V.; Suter, J. L.; Greenwell, H. C. Chem. Mater. 2007, 19, 5510. (53) Kim, N.; Kim, Y.; Tsotsis, T. T.; Sahimi, M. J. Chem. Phys. 2005, 122, No. 214713. (54) Wang, J. W.; Kalinichev, A. G.; Kirkpatrick, R. J.; Hou, X. Q. Chem. Mater. 2001, 13, 145. (55) Wang, J. W.; Kalinichev, A. G.; Amonette, J. E.; Kirkpatrick, R. J. Am. Mineral. 2003, 88, 398. 1897



ARTICLE

(56) Kim, N.; Harale, A.; Tsotsis, T. T.; Sahimi, M. J. Chem. Phys. 2007, 127, No. 224701. (57) Greenwell, H. C.; Jones, W.; Coveney, P. V.; Stackhouse, S. J. Mater. Chem. 2006, 16, 708. (58) Broughton, J. Q.; Gilmer, G. H.; Jackson, K. A. Phys. Rev. Lett. 1982, 49, 1496. (59) Churakov, S. V.; Iannuzzi, M.; Parrinello, M. J. Phys. Chem. B 2004, 108, 11567. (60) Campana, C.; Mussard, B.; Woo, T. K. J. Chem. Theory Comput. 2009, 5, 2866. (61) Allmann, R.; Jepsen, H. P. Neues Jahrb. Mineral., Monatsh. 1969, 544. (62) Zangi, R.; Engberts, J. B. F. N. J. Am. Chem. Soc. 2005, 127, 2272. (63) Kirkpatrick, R. J.; Kalinichev, A. G.; Wang, J. Mineral. Mag. 2005, 69, 289. (64) Bennett, C. H. J. Comput. Phys. 1976, 22, 245. (65) Delgado-Buscalioni, R.; De Fabritiis, G.; Coveney, P. V. J. Chem. Phys. 2005, 123, No. 054105. (66) Widom, B. J. Chem. Phys. 1963, 39, 2808. (67) Cantrell, D. G.; Gillie, L. J.; Lee, A. F.; Wilson, K. Appl. Catal., A 2005, 287, 183.

1898


Elucidating the Crystal Face- and Hydration ... - ACS Publications

Recommend Documents