Hydrophobicity and Hydrophilicity Balance Determines Shape

Beijing National Laboratory for Molecular Sciences (BNLMS), Institute of Chemistry, ... of the shape selectivity of a reactant has been proposed and v...
2 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCC

Hydrophobicity and Hydrophilicity Balance Determines Shape Selectivity of Suzuki Coupling Reactions Inside Pd@meso-SiO2 Nanoreactor Yuejie Ai,*,† Xin Li,*,‡ Yongfei Ji,‡ Wei-Guo Song,§ and Yi Luo‡,∥ †

School of Environment and Chemical Engineering, North China Electric Power University, Beijing 102206, P.R. China Division of Theoretical Chemistry and Biology, School of Biotechnology, KTH Royal Institute of Technology, SE-10691 Stockholm, Sweden § Beijing National Laboratory for Molecular Sciences (BNLMS), Institute of Chemistry, Chinese Academy of Science, Beijing 100190, P. R. China ∥ Hefei National Laboratory for Physical Science at the Microscale, University of Science and Technology of China, Hefei, 230026 Anhui. P. R. China ‡

S Supporting Information *

ABSTRACT: Molecular sorting and catalysis directed by shape selectivity have been extensively applied in porous extended frameworks for a low-carbon, predictable, renewable component of modern industry. A comprehensive understanding of the underlying recognition mechanism toward different shapes is unfortunately still missing, owing to the lack of structural and dynamic information under operating conditions. We demonstrate here that such difficulties can be overcome by state-of-the-art molecular dynamics simulations which provide atomistic details that are not accessible experimentally, as exemplified by our interpretation for the experimentally observed aggregationinduced shape selectivity for Suzuki C−C coupling reaction catalyzed by Pd particles in mesoporous silica. It is found that both aggregation ability and aggregating pattern of the reactants play the decisive role in controlling the shape selectivity, which are in turn determined by the balance between the hydrophobicity and hydrophilicity of the reactants, or in other words, by the balance between the noncovalent hydrogen bonding interaction and van der Waals forces. A general rule that allows prediction of the shape selectivity of a reactant has been proposed and verified against experiments. We show that molecular modeling is a powerful tool for rational design of new mesoporous systems and for the control of catalytic reactions that are important for the petrochemical industry.



INTRODUCTION Molecular shape selectivity, by which a catalyst selectively chooses reaction species with certain shapes, has been the crown jewel of chemistry. Nature itself exhibits an art of molecular shape selectivity, which has been evolutionarily refined over millions of years. Owing to their exquisite molecular recognition properties, the enzymes induce shapespecific substrates to fit snugly into their enzymatic cavities. Learning from Nature, the shape selectivity has been applied in chemical industry.1,2 For instance, in petroleum processing, molecular shape-based selective sorting and catalysis employ extended porous frameworks of zeolites,1 hybrid metal−organic frameworks (MOFs),3 and even certain organic molecules.4 The most successful application of shape selectivity in industry lies in zeolites, which are microporous mineral materials.1,2 The shape of the internal pore structure of a zeolite can strongly affect the selectivity with which the shape of the reactants, intermediate species, and products are critical.1,5 However, such zeolite-based shape selectivity is hitherto confined to a single © XXXX American Chemical Society

small molecular behavior in the zeolite channels (less than 1 nm diameter) and gas-phase environment. However, in organic synthesis, a large number of catalytic reactions involve much larger complex intermediates and products. These reactions usually take place in liquid medium, and shape selectivity has been an elusive goal. Mesoporous materials opened up new possibilities for shape selectivity in liquid reactions. Because of their large specific pore volume and narrow pore size distribution, mesopores (between 2 and 50 nm) have an advantage over microporous materials for the adsorption and transformation of large organic molecules.6 In a recent study, Song et al.7 observed impressive shape selectivity over a well-designed Pd@meso-SiO2 nanoreactor which was composed of mesoporous silica hollow spheres and Pd nanoparticles inside the spheres (Scheme 1). Received: January 19, 2016 Revised: April 2, 2016

A

DOI: 10.1021/acs.jpcc.6b00616 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

qi and qj are the atomic charges of both atoms, and Aij, Bij, and Cij are constants specified for each type of pair of atoms i and j. These parameters were adopted from ref 10. A 24-6 LennardJones potential was further added to prevent two atoms from getting very close to each other where the original BKS potential cannot provide realistic repulsion. This modification has been shown to improve the reliability of the model without affecting its accuracy.11

Scheme 1. Suzuki Reaction of Substituted Phenylboronic Acid (PBA) and Iodobenzene Catalyzed by Pd@meso-SiO2 Nanoreactor in ref 7a

V (rij) =

a

qiqj rij

⎡⎛ ⎞24 ⎛ ⎞6 ⎤ σij σij + Aij exp(− Bij rij) − 6 + 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r rij ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠ (1) Cij

The geometries of the reactants, namely iodobenzene and arylboronic acids, were optimized at the B3LYP/6-31G(d) level of theory, as implemented in the Gaussian09 program package.12 Based on the optimized geometries, the electrostatic potentials (ESP) of the molecules were then calculated at the HF/6-31G(d) level of theory according the Merz−Singh− Kollman scheme.13,14 On the basis of the restrained electrostatic potential (RESP), partial atomic charges were derived and employed in subsequent MD simulations.15,16 Here, the General Amber Force Field (GAFF)16 was employed to model the organic reactants. For boronic acid, the LennardJones parameters of boron atoms were adopted from ref 17, while the bonded parameters were derived from the Hess2FF scheme.18,19 Classical MD Simulations. The initial structure for the meso-SiO2 pore was built from an α-quartz phase20 whose unit cell is hexagonal. In order to make an orthorhombic supercell, the unit cell was first transformed to an orthorhombic one by taking the new lattice vectors as a′ = 5(a + b), b′ = 8(−a + b), and c′ = c. Several silicon and oxygen atoms were removed manually in a circular direction with a diameter about 2 nm at the center of the new unit cell in order to build a mesoporous pore in the nanoreactor. This large unit cell was then replicated along the c direction 10 times to yield a mesoporous silica system of size 42.5 Å × 39.3 Å × 54.1 Å with 5970 atoms (1990 silicon atoms and 3980 oxygen atoms). The classical MD simulations consist of two parts. In the first part, we aim to investigate the diffusion behavior of the reactants in the meso-SiO2 nanoreactor filled with ethanol. Here, the simulation box was enlarged by 5 nm along the z-axis, after which the reactant molecules (one iodobenzene molecule and 10 phenylboronic acid (PBA) (or 3-BIPBA) molecules, or one iodobenzene molecule and five n-(C4H9)-PBA (or t(C4H9)-PBA) molecules) were randomly inserted. The system was then filled with ethanol molecules using the genbox tool of the Gromacs program package21 and then subjected to energy minimization and constant-NpT simulation with a time step of 1 fs. The temperature of the system was maintained at 350 K by means of the Nosé−Hoover thermostat,22,23 and all dimensions of the simulation pore were coupled independently (anisotropic scaling) to reference pressures of 1 bar using the Parrinello− Rahman scheme.24,25 During the simulations, the LennardJones potential was truncated at 1.0 nm, while the Coulomb interaction was treated by the particle-mesh Ewald (PME)26 with a real-space cutoff of 1.0 nm. The simulation was conducted for 100 ns with coordinates saved ever 5 ps. Diffusion constants of organic reactants were computed from the mean square deviation of the ensemble of atoms according to Einstein’s relation, using the g_msd tool of the Gromacs Program package.21

Yields of the catalyzed reactions are given in parentheses.

On this catalyst, the Suzuki coupling reactions preceded with sharp shape selectivity; i.e., n-butyl-substituted phenylboronic acid led to 86% yield, while tert-butyl-substituted phenylboronic acid led to almost zero yield. It has been proposed that such shape selectivity of the nanoreactor arises from a collective diffusion barrier induced by the preferential adsorption of bulky phenylboronic acid onto mesoporous silica. 7 Such an explanation is convenient and empirical. Unfortunately, it does not provide much insight about the origin of the shape selectivity at the molecular level. In this contribution, we employ molecular dynamics (MD) simulations to explore the above-mentioned catalytic system. In particular, two aspects have been addressed. One is the confinement effect of the mesopores. Compared to the extensively studied silica surface,8 the disordered atomic structures and the presence of nanopores provide large matrix/pore interfacial area and resistance to energy transport through the nanostructure.9 The other is focused on the aggregation among reacting species which may serve as a way to create local obstruction and to affect the shape selectivity inside the nanoreactor. Through quantum chemical calculations and MD simulations, we propose a generally applicable volcano curve showing a linear relation between the effective size of clustering and the hydrophobic parameter. The theoretical simulations, allied to recent experimental studies, pioneer the rational design of a wide range of new shape selective reactions that the chemical industrial processes demand.



COMPUTATIONAL DETAILS Force Field and Parameters. We employed the BKS potential10 to model the mesoporous silica nanoreactor, as shown in eq 1. The first three terms compose the original BKS potential, where rij is the distance separating two atoms i and j, B

DOI: 10.1021/acs.jpcc.6b00616 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

while the hydroxyl group was treated as hydrophilic contribution. A sizable hydrophobic contribution to ΔSASA, arising from the hydrophobic groups being encapsulated inside the aggregation, indicates a strong hydrophobicity of the molecule that drives the formation of aggregation.

In the second part, we aim to investigate the aggregation abilities of different reactant molecules in ethanol solution. For each arylboronic acid reactant, 20 molecules were solvated in a cubic box with dimensions of approximately 50 Å by ethanol molecules and then subjected to energy minimization and constant-NpT simulation with the same simulation parameters as mentioned in the first part, except that the isotropic pressure coupling scheme was employed. The simulations were conducted for 50 ns for each arylboronic acid, and snapshots were extracted from the last 10 ns of trajectory. To assess the shape and size of the aggregation formed by each arylboronic acid, which is essential in determining their ability of entering the mesopore of silica nanreactor, we first computed the radius of gyration tensor defined in analogy to the moments of inertia ⎡ x2 ⎢∑ i 1⎢ S = ⎢ ∑ xiyi N⎢ ⎢⎣∑ xizi



RESULTS AND DISCUSSION Comparison of Different Aggregation Ability: Transport of Iodobenzene Together with PBA and 3-BIPBA in Mesoporous Pore and Ethanol Solution Box. There are several factors that may influence the reaction efficiency of Suzuki reaction in silica mesopores. First, the confinement within the silica mesopores sets up a confined environment for an efficient relative orientation of the reagents and the catalyst and will enhance the concentration and activity of those accessible Pd surfaces. Second, as illustrated in ref 7, since reactions are confined within the mesopores, it must diffuse through the mesopores. Thus, another important factor is the transport of the reactants from one end to the other where the catalyst Pd particles reside. Considering the relatively low concentration of iodobenzene in the mesoporous silica filled with solution molecules and boric acid compounds, the mobility of the iodobenzene molecules will be the main important factor that determines the rate of reactions. As a starting point, we address the different shapes of compounds phenylboronic acid (PBA) and 3-biphenylboronic acid (3BIPBA) by investigating their diffusion in mesopore and adsorption onto silica surface using classical MD simulations. From the calculated self-diffusion constants for iodobenzene as shown in Table 1, one may conclude that its self-diffusivity in

∑ xiyi ∑ xizi ⎤⎥ ⎥

∑ yi2 ∑ yzi i ⎥ ⎥

∑ yzi i ∑ zi2 ⎥⎦

(2)

where S is the radius of gyration tensor, N is the number of atoms in the maximal cluster, and {xi, yi, zi} is the position vector of each atom with respect to the center of mass of the cluster. Diagonalization of S gives the three eigenvalues (λ1, λ2, and λ3) as well as the effective radius along each principal axis (a, b, and c), from which the triaxial ellipsoid shape parameter G is obtained.27 a=

3λ1

b=

3λ 2

c=

3λ3

1 ⎡ a 2 + b2 + c 2 ⎤ ⎥ G= ⎢ 5 ⎣ (abc)2/3 ⎦

Table 1. Self-Diffusion Constant D of the Iodobenzene in Meso-SiO2 Nanoreactor with Different Reactantsa

(3)

(4)

Then we introduce the effective size of the maximal cluster, defined as n neff = max (5) G

a

reactant

D(iodobenzene) (10−7 cm2 s−1)

PBA 3-BIPBA n-(C4H9)- PBA t-(C4H9)- PBA

11.27 (0.50) 2.69 (0.80) 0.73(0.03) 0.69 (0.26)

Standard errors of mean are shown in parentheses.

PBA packed mesopore (11.27 × 10−7 cm2 s−1) is much higher than that of 3-BIPBA (2.69 × 10−7 cm2 s−1). This is further confirmed by the relative displacements of iodobenzene molecule in mesoporous silica in which the moving range for iodobenzene in PBA is covered from 0 to 5 nm along the whole mesoporous pore, whereas the transport of iodobenzene in 3BIPBA is confined within a range of around 2 nm, as shown in the Supporting Information (SI, Figure S1). In order to identify the fundamental interactions and processes that control this sharp shape selectivity, we then give detailed structural characteristics inside the mesopores, as depicted in Figure 1. Adsorption of ethanol molecules and boric acid compounds onto the mesoporous silica surface is observed to occur with the help of hydrogen bonds between the hydroxyl group of boronic acid and the oxygen atoms of silica surface. MD simulation reveals that the 3-BIPBA molecule is packed tightly inside the mesopore channels and blocks the transport channel for the iodobenzene molecule, as magnified in Figure 1b. On one hand, the larger molecular volume of 3-BIPBA molecule makes it more likely to block the iodobenzene. On the other hand, the 3-BIPBA molecules are able to form tighter dimers than PBA molecules. We optimized the geometries for

where nmax is the number of non-hydrogen atoms in the maximal cluster and G is the triaxial ellipsoid shape parameter. The value of parameter G is 0.6 for ideal sphere and increase as the shape of cluster deforms. Thus, in eq 5, nmax reflects the actual size of the maximal cluster, while the parameter G takes into account the flexible shape of cluster which provides the possibility of deforming and entering a pore with a smaller diameter. Solvent accessible surface area (SASA) is another indicator that characterizes the degree of aggregation. The SASA of a group of molecules shrinks upon formation of aggregation since part of the molecular surface becomes buried inside the cluster. The shrinkage of SASA, or ΔSASA, defined as the difference between the SASA of the whole aggregation and the sum of SASA of each constituent, reflects the degree of aggregation. For each arylboronic acid, we used the g_sas tool of the Gromacs program package21 to compute the hydrophilic and hydrophobic contributions to ΔSASA. The “g_sas” computes hydrophobic, hydrophilic, and total solvent accessible surface area. It is based on the method from Eisenhaber et al.28 In our work, the aromatic rings were treated as hydrophobic groups, C

DOI: 10.1021/acs.jpcc.6b00616 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

molecules. Then two phenyl rings of 3-BIPBA are expected to interact with the surrounding 3-BIPBA molecules and are likely to exhibit high hydrophobicity inside the mesopores. In contrast, the PBA molecules consisting of only one phenyl ring form weaker π−π stacking (3 kcal/mol weaker than that in 3-BIPBA dimer). As shown in Figure 2, there is a horizontal displacement of around 1.52 Å between the two corresponding carbons of benzenes. The average horizontal displacement of heavy atoms between two PBA molecules is 1.78 Å (see Table S2 in SI). Therefore, the relatively loose stacking structure of PBA molecules gives more free space for the iodobenzene to go through, as shown in Figure 1c.The final outcome of the above factors is the aggregation ability of different boric acid compounds. In order to further investigate the different aggregation abilities between these two boric acids, we studied their aggregation properties in ethanol solution. The calculated dynamic properties of geometric clustering including the size of the maximal cluster, the number of the clusters in total, as well as the computed hydrophobic, hydrophilic, and total solvent accessible surface area are depicted in Figure 3. The black line with a larger maximal cluster than the red one signifies the stronger aggregation ability of the 3-BIPBA molecule in ethanol solution. Correspondingly, the number of clusters of 3-BIPBA in ethanol is smaller than that of PBA. In adition, despite of the similar hydrophilic contributions to SASA (about 15 nm2) for both molecules, for 3-BIPBA the hydrophobic contribution (around 50 nm2) accounts for a dominant proportion of the total SASA (around 65 nm2). Only a small portion of hydrophobic SASA of about 30 nm2 was observed for the PBA clusters. As the aggregation of the molecules is strongly correlated to the size of the maximal cluster and also to the hydrophobic SASA, we propose that the observed shape selectivity between PBA and 3-BIPBA molecules is not only because of the molecular size but more importantly caused by their different aggregation ability induced by hydrophobic properties. We also argue that despite the lack of full consideration of real mesoporous silica materials the approach of MD simulations in both simple mesopore and solution serves as a useful starting point that can help us arrive at a further description of mesoporous catalysis in solution. Here, in mesoporous silica, the effect of aggregation ability gives a new interpretation of the traditional concept of shape selectivity, which predicts that the 3-BIPBA molecule will not react with iodobenzene because their aggregation can effectively block the mesoporous pores. Still, the hydrophobicity-based shape selectivity encounters difficulties in interpreting the yield of n-(C4H9)-PBA and t(C4H9)-PBA in mesoporous silica reactor, as the traditional shape selectivity concept does. In the next section, we will address the aggregation pattern of these two molecules in mesoporous silica which gives rise to their distinct reaction yield, despite their similar structures and aggregation abilities. As the separation of branched hydrocarbons from linear hydrocarbons in the crude oil is of considerable interest in petroleum industry, this comparative example displays how silica can help to control reaction activities by favoring the formation of particular pattern of aggregation. Comparison of Different Aggregation Patterns: Transport of Iodobenzene Together with n-(C4H9)-PBA and t-(C4H9)-PBA in Mesoporous Pore and Ethanol Solution Box. The aggregation properties of n-(C4H9)-PBA and t-(C4H9)-PBA that related to the size of the maximal

Figure 1. Side (a) and magnified (b) views for the aggregation of 3BIPBA that blocked the passage of iodobenzene taken from one snapshot of the MD trajectory. The 3-BIPBA is shown in balls and sticks, the ethanol molecules are shown in lines, while the iodobenzene is drawn with VDW sphere. The snapshot was visualized using VMD.29 (c) The magnified views for the aggregation of PBA in mesopore.

dimers of PBA and 3-BIPBA at the MP2/cc-pVTZ level of theory without any geometric constraints using the Gaussian09 program.12 The calculated binding energy between two neighboring 3-BIPBA molecules (11.48 kcal/mol) is larger than that of PBA (8.33 kcal/mol). The corrected interaction energy considering basis set superposition error (BSSE) is 7.31 and 5.93 kcal/mol for the 3-BIPBA dimer and PBA dimer, respectively (see Table S1). Different configurations with various dimer orientations result in interaction energies of about 1.34 to 5.93 kcal/mol for PBA and 7.31 to 13.71 kcal/ mol for 3-BIPBA dimers, which can be found in Figure S2 and Table S1. The benzene rings in the two neighboring 3-BIPBA molecules are found to stay close to the optimal distances (3.498 and 3.615 Å) of π−π stacking, as shown in Figure 2. Moreover, these weak interactions are further complemented by hydrophobic interactions in the presence of polar solvent

Figure 2. Structures for PBA, 3-BIPBA, and their dimers optimized at the MP2/cc-pVTZ level of theory in combination with a polarizable continuum model (IEF-PCM) of ethanol. All all the geometrical parameters are optimized without any constraints. D

DOI: 10.1021/acs.jpcc.6b00616 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 3. Dynamical properties analyses (max cluster size, number of clusters, and solvent-accessible surface area (SASA, in nm2)) for the PBA and 3-BIPBA in the ethanol solution boxes.

cluster and also the solvent accessible surface are modeled within the ethanol solution. The corresponding dynamic properties are presented in Figure 4. From the pictures, we

Figure 5. Magnified views for the aggregation of n-(C4H9)-PBA and t(C4H9)-PBA molecules in the ethanol solution box (left column) and mesoporous pore (right column).

molecules is mainly driven by the hydrophobic properties of tertiary butyl groups. In mesoporous silica, the Suzuki coupling reaction shows selectivity toward the linear n-(C4H9)-PBA over the branched t(C4H9)-PBA molecules, as shown in Scheme 1. One may expect from the possible industrial applications that the branched hydrocarbons should be able to be separated from the linear ones. We then studied the diffusion of iodobenzene in mesoporous pore filled with n-(C4H9)-PBA and t-(C4H9)PBA molecules. The self-diffusion constants of iodobenzene are quite close; see Table 1. However, the displacements of iodobenzene molecule in Figure S1 show that it moved past the pore over 2 nm together with linear n-(C4H9)-PBA while it is confined within a range of around 1.5 nm inside the mesopore that was filled with the branched t-(C4H9)-PBA. It thus suggests that the iodobenzene has a higher self-diffusivity in the presence of the linear n-(C4H9)-PBA as compared to t-(C4H9)-PBA in mesopore. This behavior is in agreement with the observed yields of Suzuki coupling reactions. A careful inspection into

Figure 4. Dynamical properties analyses (max cluster size, number of clusters, and solvent-accessible surface area (SASA, in nm2)) for the n(C4H9)-PBA and t-(C4H9)-PBA in the ethanol solution boxes.

concluded that the aggregation abilities of these two molecules in ethanol solution are quite similar. On one hand, the black and red lines of max cluster size almost overlap, suggesting the similar aggregation abilities for both kinds of molecules. On the other hand, in addition to the similar total SASA value which is around 60 nm2, the hydrophilic or hydrophobic contributions are also comparable for both n-(C4H9)-PBA and t-(C4H9)-PBA clusters. However, when we take a close look at the structures, different aggregation patterns were found between the linear and branched molecules. As illustrated in the left column of Figure 5, the aggregation of n-(C4H9)-PBA molecules in ethanol solution took place via the hydrogen bonding between the hydroxyl groups, whereas the aggregation of t-(C4H9)-PBA E

DOI: 10.1021/acs.jpcc.6b00616 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

hydrophobic groups are substituted on the PBA molecule. Conversely, a negative HHB corresponds to the presence of hydrophilic substituents and the dominant role of hydrophilic effects. An imbalance between the hydrophobic and hydrophilic effects results in large absolute value of HHB and in turn generates a strong driving force toward aggregation. As shown in Figure 6 with blue columns, the PBA and the 3-BIPBA have

the mesopore suggests that the iodobenzene was blocked by tertiary butyl groups of two t-(C4H9)-PBA molecules, as illustrated in Figure 5. The aggregated rigid tertiary butyls form a wall that the iodobenzene molecule gets rebounded back from when it comes close; see detailed snapshots in Figure S3. On the contrary, as shown in Figure 5, the intramolecular hydrogen bonds between n-(C4H9)-PBA molecules are broken due to formation of the intermolecular hydrogen bonds between n(C4H9)-PBA molecules and the silica surface. This is further verified by the quantitative dynamic analysis for the number of hydrogen bonds shown in Figure S4. As most of the n-(C4H9)PBA molecules form hydrogen bonds with the mesoporous surface, the linear tail are highly flexible and may orient themselves when an iodobenzene molecule approaches (Figure S3). Thus, the iodobenzene may come through the linear tails and reach the end of mesopore to meet the Pd catalyst. As shown in Figure S5, the larger fluctuant distances between the butyl groups of n-(C4H9)-PBA indicate that compared to a rigid tertiary butyl of t-(C4H9)-PBA, the linear n-butyl is more flexible and adaptable to local changes of its surroundings, which serves as a favorable factor for the transport. Therefore, we propose that in addition to the aggregation ability the aggregation pattern plays an important role in affecting the diffusive behavior inside the mesoporous silica. Here again, we reaffirmed the importance of the specific structure and property of mesoporous silica which supplies not only the spatial construction but also some functional groups for intermolecular interactions. On the basis of the discussions above, one can conclude that the shape selectivity observed in Pd@meso-SiO2-catalyzed Suzuki coupling reaction is mainly induced by the aggregation of the reactants with different stabilities and patterns. In addition, at times, it is expected that the specific characteristic of mesoporous silica and the molecular reorientations during diffusion that caused by molecular flexibility are both important for the transport process. General Rule for the Aggregation and Shape Selectivity. With comprehensive modeling, the origins of the shape selectivity for two sets of molecules have been identified, which allow us to explain the corresponding experimental results. However, it would be more desirable if we can find a more general rule that can correctly predict the reaction yield without carrying out all those detailed and expensive computational modeling. As we have found, the aggregation behavior of the molecules in mesopore is closely related to two inter-related important factors: the size of the maximal cluster and the hydrophobicity of molecules in solution. It can be seen that the maximal cluster size is a direct evidence of the aggregation ability, which is mainly driven by the hydrophobicity of the molecules in solution. In other words, the large maximal cluster size is the result of high hydrophobicity. However, both hydrophobic and hydrophilic groups are present in a molecule and compete with each other in solution. To some extent, the hydrophobicity will be offset by the hydrophilic groups. Therefore, in order to provide a realistic reflection of the hydrophobic ability, we here define a parameter, named as the hydrophobicity−hydrophilicity balance (HHB), i.e., ΔSASA (hydrophobic) − ΔSASA (hydrophilic), which is the difference between the variations of hydrophobic and hydrophilic solvent accessible surfaces area (ΔSASA) upon aggregation. If taking the PBA molecule as the standard reference molecule that has no substituents with a HHB of 0.26, a more positive value of HHB indicates the

Figure 6. Max cluster size and deviations of hydrophobic and hydrophilic solvent accessible surfaces area (ΔSASA) for PBA and 3BIPBA in the ethanol solution boxes, averaged over the last 10 ns of MD trajectories.

average sizes of maximal clusters of 2.49 and 3.40, respectively. More important, a larger HHB value for 3-BIPBA indicates its stronger aggregation ability than PBA. The maximal cluster size and HHB are then considered for other seven reactants (p-, m-, and o-carboxyphenylboronic acid (p-, m-, o-CAPBA), 2,4,6-tri-tert-butylphenylboronic acid (tri-t(C4H9)-PBA), 1-naphthaleneboronic acid (Naph-PBA), and n(C4H9)-PBA, t-(C4H9)-PBA)) that have been studied experimentally.7 All of the dynamical analyses are summarized in Figure S6. On the basis of the original data of maximal cluster size, considering the flexible shape of molecule itself, we introduce the effective size of cluster neff (see definition of eq 5 in the computational details). As a complement to the maximal cluster size, the parameter neff gives a more realistic evaluation of the actual molecular volume of the molecule assembly entering a pore with a smaller diameter. As shown in Figure 7, the HHB value demonstrates a good correlation with the effective size of the cluster for all molecules under investigation. We have also plotted in Figure 7 the relationship between the HHB value and the yield of nine reactants measured in ref 7, represented as a volcano curve. It simply implies that the yield of the reactant can be characterized by its absolute HHB value. A larger HHB value leads to better aggregation and in turn leads to a low reaction yield. It is noted that a positive HHB value represents the hydrophobic substituents like phenyl and alkyl groups, while the negative ones correspond to the hydrophilic substituent which then carries carboxylic acid groups. The limit case is represented by tri-t-(C4H9)-PBA, which is at the bottom-right of yield volcano and has the largest value of both two parameters. It is interesting to observe that the more negative HHB value means higher potential for molecular aggregation through hydrogen bonds between hydrophilic substituents, whereas the strong molecular aggregation via hydrophobic groups in solution is mainly dependent on the van der Waals forces F

DOI: 10.1021/acs.jpcc.6b00616 J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C



Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b00616. Dynamic properties analyses for the nine reactants in the ethanol solution box, etc. (PDF)



AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.



Figure 7. Two-dimensional volcano plot (according to the Sabatier principle that relates the optimal concentrations of catalysts and substrates30,31) of HHB as a function of effective cluster size neff (the dotted curve) and corresponding production yield of nine reactants (the solid curve) for the Suzuki reaction in mesoporous pore from ref 7. The color column represents nine different reactants. The slash area represents the active zone. The data are achieved by averaging over the last 10 ns of MD trajectories.

ACKNOWLEDGMENTS This work is supported by the National Natural Science Foundation of China (21403064, 20925311, 21161160557) and the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine. The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC and NSC.



between them. Essentially, HHB is a measure of the competition between these noncovalent interactions. We can further mark out an active zone of the reactants for the Suzuki reaction in mesoporous pore. That is to say, once the computed values of two parameters (neff and HHB) for reactants inside the active zone, the reaction will likely take place. The closer to the center of the active zone, the greater the product yield. This clearly gives us a great power to predict the reactivity through simple theoretical simulations.



REFERENCES

(1) Smit, B.; Maesen, T. L. M. Towards a Molecular Understanding of Shape Selectivity. Nature 2008, 451, 671−678. (2) Vermeiren, W.; Gilson, J. P. Impact of Zeolites on the Petroleum and Petrochemical Industry. Top. Catal. 2009, 52, 1131−1161. (3) Li, J.-R.; Sculley, J.; Zhou, H.-C. Metal-Organic Frameworks for Separations. Chem. Rev. 2012, 112, 869−932. (4) Mitra, T.; Jelfs, K. E.; Schmidtmann, M.; Ahmed, A.; Chong, S. Y.; Adams, D. J.; Cooper, A. I. Molecular Shape Sorting Using Molecular Organic Cages. Nat. Chem. 2013, 5, 276−281. (5) Zaera, F. The New Materials Science of Catalysis: Toward Controlling Selectivity by Designing the Structure of the Active Site. J. Phys. Chem. Lett. 2010, 1, 621−627. (6) Hartmann, M. Ordered Mesoporous Materials for Bioadsorption and Biocatalysis. Chem. Mater. 2005, 17, 4577−4593. (7) Chen, Z.; Cui, Z.-M.; Li, P.; Cao, C.-Y.; Hong, Y.-L.; Wu, Z.-y.; Song, W.-G. Diffusion Induced Reactant Shape Selectivity Inside Mesoporous Pores of Pd@meso-SiO2 Nanoreactor in Suzuki Coupling Reactions. J. Phys. Chem. C 2012, 116, 14986−14991. (8) Rarivomanantsoa, M.; Jund, P.; Jullien, R. Classical Molecular Dynamics simulations of Amorphous Silica Surfaces. J. Phys.: Condens. Matter 2001, 13, 6707−6718. (9) Coquil, T.; Fang, J.; Pilon, L. Molecular Dynamics Study of the Thermal Conductivity of Amorphous Nanoporous silica. Int. J. Heat Mass Transfer 2011, 54, 4540−4548. (10) Vanbeest, B. W. H.; Kramer, G. J.; Vansanten, R. A. Force-fields for Silicas and Aluminophosphate based on Abinitio Calculations. Phys. Rev. Lett. 1990, 64, 1955−1958. (11) Guissani, Y.; Guillot, B. A Numerical Investigation of the Liquidvapor Coexistence Curve of Silica. J. Chem. Phys. 1996, 104, 7633. (12) Frisch, M. J., et al. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009. (13) Singh, U. C.; Kollman, P. A. An approach to Computing Electrostatic Charges for Molecules. J. Comput. Chem. 1984, 5, 129− 145. (14) Besler, B. H.; Merz, K. M.; Kollman, P. A. Atomic Charges Derived from Semiempirical Methods. J. Comput. Chem. 1990, 11, 431−439. (15) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A wellbehaved Electrostatic Potential based Method using Charge Restraints for Deriving Atomic Charges - The RESP model. J. Phys. Chem. 1993, 97, 10269−10280.

CONCLUSIONS

Rational and predictive engineering of highly selective catalytic reactions constitutes a key step in future development of synthesis and characterization of controllable nanomaterials and has triggered considerable efforts in industrial and laboratoryorientated pursuing toward energy’s sustainability. The shape selectivity plays a significant role in molecular shape sorting and catalysis of chemical industry for producing petrochemical product in a controlled, predictable way. In this work, we have assessed and predicted the aggregation-induced shape selectivity of the Suzuki C−C coupling reaction taking place in a mesoporous silica nanoreactor by means of theoretical simulations for various reactants. We also suggest that the shape selectivity is promoted by the aggregation of molecules via their hydrophobic properties in solution, which is further supported by a volcano curve connecting the effective size of cluster and hydrophobic parameters. On the basis of this volcano curve derived from the results of MD simulations, mechanistic interpretations were provided for the observed different reaction yields of each arylboronic acid. We believe that the approach described here opens up a promising way to predict shape selectivity of synthesis reactions in mesoporous nanostructures and will likely lead to rapid developments in selective catalytic synthesis and separation in experimental and industrial conditions as well as further theoretical simulations in elucidating the underlying fingerprints and mechanisms. G

DOI: 10.1021/acs.jpcc.6b00616 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (16) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force field. J. Comput. Chem. 2004, 25, 1157−1174. (17) Tafi, A.; Agamennone, M.; Tortorella, P.; Alcaro, S.; Gallina, C.; Botta, M. AMBER Force field Implementation of the Boronate Function to Simulate the Inhibition of Beta-lactamases by Alkyl and Aryl Boronic acids. Eur. J. Med. Chem. 2005, 40, 1134−1142. (18) Seminario, J. M. Calculation of Intramolecular Force fields from Second-derivative Tensors. Int. J. Quantum Chem. 1996, 60, 1271− 1277. (19) Nilsson, K.; Lecerof, D.; Sigfridsson, E.; Ryde, U. An Automatic Method to Generate Force-field Parameters for Hetero-compounds. Acta Crystallogr., Sect. D: Biol. Crystallogr. 2003, 59, 274−289. (20) Levien, L.; Prewitt, C. T.; Weidner, D. J. Structure and Elastic Properties of Quartz at Pressure. Am. Mineral. 1980, 65, 920−930. (21) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. Gromacs 4: Algorithms for Highly Efficient, Load-balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (22) Nose, S. A Molecular-dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (23) Hoover, W. G. Canonical Dynamics- Equilibrium Phase- Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (24) Parrinello, M.; Rahman, A. Polymorphic Transitions in Singlecrystals - A New Molecular-Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (25) Nose, S.; Klein, M. L. Constant Pressure Molecular-dynamics for Molecular-systems. Mol. Phys. 1983, 50, 1055−1076. (26) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (27) Harding, S. E. A General-method for Modeling Macromolecular Shape in Solution - A Graphical (II-G) Intersection Procedure for Triaxial Ellipsoids. Biophys. J. 1987, 51, 673−680. (28) Eisenhaber, F.; Lijnzaad, P.; Argos, P.; Sander, C.; Scharf, M. The Double Cubic Lattice Method: Efficient Approaches to Numerical Integration of Surface Area and Volume and to Dot Surface Contouring of Molecular Assemblies. J. Comput. Chem. 1995, 16, 273−284. (29) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. J. Mol. Graphics 1996, 14, 33−38. (30) Rothenberg, G. Catalysis: Concepts and Green Applications; Wiley-VCH, 2008. (31) Knözinger, H.; Kochloefl, K. Heterogeneous Catalysis and Solid Catalysts; Wiley−VCH, 2005.

H

DOI: 10.1021/acs.jpcc.6b00616 J. Phys. Chem. C XXXX, XXX, XXX−XXX