Pore diameter dependence of catalytic activity: p-nitrobenzaldehyde conversion to an aldol product in amine-functionalized mesoporous silica

There has long been heuristic recognition that the occurrence of single-file diffusion (SFD) inside catalytically-active nanopores can lead to a strong reduction in reactivity for catalytic conversion.^1–5 For SFD, molecular species cannot pass each other due to steric constraints,^6,7 and its occurrence has been demonstrated for diffusion of various molecular species in zeolites.^8,9 For irreversible conversion, reactants enter the pore and convert to products at catalytic sites, but SFD impairs product removal and the entire pore center becomes exclusively populated with products.¹⁰ Consequently, reactant population, and thus reactivity, is limited to near pore openings. For reversible conversion, just the excess reactant concentration above the equilibrium value is localized near pore openings. However, the same conclusion applies regarding limited reactivity. This behavior can be verified by coarse-grained stochastic modeling and Kinetic Monte Carlo (KMC) simulation^1–5,10 and also by a generalized hydrodynamic theory which captures the subtle interplay between SFD and fluctuations in adsorption-desorption processes at pore openings.¹¹ 

SFD applies for sufficiently narrow nanopores, but this constraint is clearly relaxed as the pore diameter is increased. Consequently, penetration of reactant into the pore, and thus reactivity, should increase with increasing pore diameter.¹⁰ This increased reactivity is controlled by the enhanced propensity for reactants and products to pass each other with increasing pore diameter. However, systematic experimental exploration of associated behavior has been lacking, as has theoretical investigation of this key passing propensity, $P$⁠, for system-specific molecular models of catalytic systems. One exception is a prescient analysis of $P$ based on transition state theory for the passing propensity of methane and ethane in zeolites.¹² 

Mesoporous silica nanoparticles (MSN) provide a natural platform within which the above phenomena can be investigated.^13–15 MSN can be readily synthesized with a range of nominal pore diameters from 2 to 10 nm. Also, MSN have been successfully functionalized to provide effective catalysts for a variety of solution-phase reactions.^16–18 It is certainly the case that the solvent phase, and the mixed solvent + reactant + product phase, inside the pore can have distinct properties than in the external fluid. For example, formation of ordered layers of water against the mesoporous silica pore walls is possible.¹⁹ However, it is also expected that for pore diameters in the range considered here at least the pore center is close to a bulk liquid-like state.¹⁹ In fact, other studies have demonstrated that the freezing temperature of water can be strongly suppressed in mesoporous silica, i.e., confinement mesopores favor a fluid state over a solid state.²⁰ In addition, there have been numerous studies quantifying solution-phase like diffusion of various molecular species inside silica mesopores, including functionalized pores.^21–23 

Of particular relevance for this study, there have been multiple studies of aldol and nitroaldol condensation reactions in functionalized MSN catalysts.^24–26 To provide a benchmark for our multiscale modeling study, we draw upon experiments for the conversion of p-nitrobenzaldehyde (PNB) to an aldol product, 4-hydroxy-(4-nitrophenyl)-2-butanone, by reaction with acetone in amine-functionalized MSN. Both acetone and hexane serve as solvents. In these studies, the MSN have diameters of 100-200 nm and are traversed by a hexagonal array of parallel linear nanopores with no connections between pores, with a range of nominal pore diameters from 2.8 nm up to about 5 nm. Reaction yield was observed to be high for the largest pore diameters but dropped dramatically for a smaller diameter of 2.8 nm. However, the longest dimensions of the reactant and product species are about 1.0 nm for PNB and 1.4 nm for the aldol product (and 0.6 nm for acetone) well below the smallest diameter. Furthermore, the shortest molecular dimensions in an orthogonal direction are significantly smaller. Thus, one might anticipate that even in the narrowest pores, there is plenty of room for passing. However, actually the effective pore diameter, d = d_p(eff), of the narrowest pore is reduced below 1.5 nm during the reaction by attachment of reactant species at the pore walls.¹⁸ We will propose that this feature allows a transition at least close to the SFD regime for the smallest d, thereby explaining the dramatic variation in reaction yield.

Theory, modeling, and simulation have the potential to provide detailed insight into the above behavior. Appropriate analysis of energetics from ab initio based formulations can potentially assess key reaction barriers and mechanisms.²⁷ Molecular Dynamics (MD) simulations are beginning to achieve sufficient fidelity to unravel specific features of cooperative reaction processes at solid-solution interfaces including assessment of the fluxional dynamics of surface functional groups.²⁸ Indeed, both these approaches have been successful in providing detailed molecular-level insight for the catalytic processes of interest here. However, the goal and nature of the current theoretical effort is somewhat different being focused on elucidating the experimentally measured reaction kinetics. Furthermore, our goal is to go beyond traditional mean-field treatments of kinetics²⁹ in order to reliably treat the subtle interplay between inhibited transport and reaction in determining catalytic yield. To address these challenges, detailed MD analysis is not viable as it cannot assess the experimental time scales for this large complex many-species system. Thus, our approach instead will utilize a novel multiscale modeling framework with distinct non-trivial components. Specifically, we will incorporate a Langevin-based assessment of molecular passing propensity in narrow pores into coarse-grained stochastic modeling of the overall reaction diffusion process.

We should emphasize several complexities associated with the key initial challenge of assessing the variation of the passing propensity, $P$ = $P$(d), with pore diameter d for the above reaction system. Both PNB reactant and aldol type product species can be regarded as elongated molecules with their longest dimensions significantly exceeding those in the orthogonal directions. Thus, for the narrowest range of pores where passing is still possible but $P$ is small, the alignment of these species with each other and also with pore axis will clearly be required for passing. As a result, both translational and rotational diffusion and any coupling between them will play a key role and must be described effectively for the accurate prediction of $P$⁠. Ideally, again, MD including a single reactant-product pair, as well as an explicit treatment of solvents and a detailed treatment of the functionalized pore wall structure, would provide the most reliable assessment of $P$⁠. Unfortunately, even here, it is difficult to use MD to assess rare passing events for narrow pores which we shall see require O(10⁵) simulation trials to precisely determine $P$⁠. We should however note that analysis of the rare dynamical events is a long-standing and still developing field^30–32 and that the associated approaches could potentially be applied to facilitate more computationally efficient analysis of rare passing events.

However, rather than performing such MD simulations with an explicit solvent, a key component of our analysis of $P$ will utilize two-particle overdamped Langevin simulations³³ with rigid reactant and product species and an implicit treatment of the solvent to track the solvent-mediated evolution of the six degrees of freedom for each of these species. This approach allows significantly longer time steps in numerical simulation of the dynamics in order to more efficiently determine $P$⁠. An illustration of the basic features of these simulations is shown in Fig. 1. This analysis will be augmented exploiting analytic results motivated by the equivalent Fokker-Planck formulation for this diffusive passing problem and associated first-passage problems.³⁴ We also utilize analytic results which provide insight into the form of $P$(d) for both large and small d where Langevin simulation is inefficient. For example, for small d, these analytic results reveal a non-trivial scaling of the form³⁴ $P$ ∼ A(d – d_c)^σ for d just above the critical value, d_c, of the pore diameter for the onset of SFD due to steric blocking.

However, as indicated above, connection with the experimentally observed reaction yield also requires incorporation of results for $P$ into a suitable multiscale modeling framework which can describe the complex heterogeneous many-particle catalytic reaction-diffusion system on the appropriate time scales and length-scales. We will apply a spatially coarse-grained stochastic model¹⁰ wherein the pores are treated as linear arrays of cells. Reactant and product species adsorb and desorb from the end cells at the pore openings, hop to adjacent unoccupied cells, and exchange between adjacent occupied cells within the pore, at appropriate rates. The rate of exchange of the PNB reactant and aldol product type species reflects $P$⁠. Within this multiscale approach, an additional issue is the appropriate or consistent determination of $P$ in the context of “parameter passing” from two-particle Langevin analysis to the coarse-grained modeling framework.

In Sec. II, we provide a brief description of the experiments and an overview of several different components of the multiscale modeling. The latter includes a description of the following: the development of the appropriate geometric models for reactant and product species; formulation of strongly-damped Langevin (i.e., Brownian) dynamics for translational and rotational degrees of freedom; determination of the 6 × 6 translation-rotation diffusion tensors required as input for the Langevin simulations; appropriate formulation of the passing propensity; and the coarse-grained modeling formulation. Section III presents detailed results from theory and simulations for the 6 × 6 diffusion tensors, for the passing propensity versus the pore diameter obtained from Langevin simulations, and for the catalytic activity versus the pore diameter obtained from the coarse-grained simulations. Conclusions are provided in Sec. IV.

We first comment further on the synthesis of the functionalized MCM-41 type MSN and on the pore wall structure during reaction. The synthesis procedure utilized a pore expander agent to achieve a range of pore diameters and also involved co-condensation wherein amine groups are incorporated during synthesis rather than by post-synthesis grafting. See Ref. 25 for further details. For the experiments on which we focus, pore diameter amine-functionalized MSN or AP-MSN have the measured values d_p(AP-MSN) = 2.8, 3.3, and 3.6 nm. These dimensions were estimated by measuring nitrogen sorption isotherms and applying Brunauer-Emmett-Teller and Barrett-Joyner-Halenda methods of analysis. However, experimental analysis indicates that the effective diameter of the pore is significantly further reduced from the above values as a result of reaction. It was proposed, and confirmed by NMR studies,²⁵ that PNB reacts with the aminopropyl groups on the functionalized pore walls to form a Schiff base. To obtain a rough assessment of the effective pore diameter during reaction, we have estimated the additional distance, δL(Schiff-amine) ≈ 0.75 nm (see the supplementary material), that the larger Schiff base protrudes toward the pore center relative to the smaller amine groups. Given that the Schiff base is attached all around the MSN pore diameter protruding further into the pore center than the amine groups, one obtains an estimate of the effective pore diameter of d = d_p(eff) ≈ d_p(AP-MSN) − 2δL(Schiff-amine) ≈ 1.3, 1.8, and 2.1 nm corresponding to the values of d_p(AP-MSN) listed above.

With this reduced effective diameter, it is reasonable to anticipate that effects of inhibited passing of reactant and product species could significantly reduce reactivity for the narrowest pore with d = 1.3 nm. It should be noted that experimental estimation of the diameter of the smallest pore after reaction by nitrogen sorption yielded a crude higher estimate around 2 nm. However, the instrument cannot assess diameters below 2 nm. Also, the configuration of the surface pendant groups on dried samples analyzed in these studies could be quite different than in the solution, and also adsorption studies may not probe the effective diameter relevant for assessing SFD. Finally, we note that formation of the Schiff base also inhibits the access of PNB to amine catalytic sites, so this effect can also reduce the local reactivation rate.

The conditions for the catalytic reaction study on which we focus were as follows: 3 mol. % catalyst, 0.39 mmol PNB, 1.5 mL acetone, and 1.5 mL hexane. Acetone is present in excess as a reactant and also as a solvent, and hexane is present as a solvent. There is limited expansion of the acetone + hexane solution upon dissolution of PNB implying that there are about X_b = 0.08 PNB molecules per nm³ in the bulk solution external to the MSN. However, PNB is not soluble in hexane, but only in acetone, which implies a preference for PNB to be located inside the MSN pores rather than in the bulk solution. Inside the pore, proximity to the pore walls can reduce the interaction of PNB with hexane. In addition, there is an affinity of PNB for the pore surface due to interaction between the carbonyl (and nitro) and the surface silanol groups. Thus, one expects a substantially higher concentration, X_p, of PNB inside the MSN pores relative to the external bulk solution. This feature will prove an important factor in determining the reaction kinetics.

The reaction was run at 60 °C for 2 h. The yield after 2h was measured in terms of percentage or fractional conversion, α, of the PNB reactant to the 4-(4-nitrophenyl)-2-butanone aldol-type product through reaction with excess acetone nearby catalytic sites inside the pores. As an aside, we note that conversion (the fraction of the starting reactant which has been lost during reaction) will not exactly match the yield (how much of a particular product is formed) since some of the PNB reactant is lost in forming the Schiff base and also since there is a small amount of secondary product (although aldol dominates under our conditions). However here we ignore this difference. Conversion, α, as a function of the effective pore diameter, d = d_p(eff), is shown in Table I. The central question is whether this dramatic, roughly 20×, decrease in yield from larger to smaller d is consistent with a transition to a (near-)SFD regime for the smallest pore diameter. As an aside, a similar variation in yield has also been found for higher catalyst concentrations, 11.3 mol. % (versus 3 mol. % used above). Again, the yield is about 2% for d = d_p(eff) = 1.3, and it increases to 65% for d = d_p(eff) = 3.2 nm. Studies for other conditions (e.g., a 3 h reaction run versus 2 h) also lead to consistent results, although they suggest some uncertainty in the estimate of the yield for the smallest d_p. One final observation is that studies of this reaction for a different class of SBA-15 type MSN (rather than MCM-41 MSN) with larger pore widths achieved conversions of around 90% after 2 h. This indicates that while the reaction is not strictly irreversible, equilibrium is strongly biased toward the product under these reaction conditions. See Sec. S1 of the supplementary material.

Our Langevin simulations treating the solvent implicitly as a continuous viscous fluid (see Sec. II C), determining the passing propensity, P, of a single reactant-product pair (see Sec. II D), will require a reasonable description of these species and of the pore geometry. To this end, each of these species is treated as a rigid molecule with steric non-overlapping constraints between the molecules themselves and also between the molecules and the pore wall. The pore is described by a cylindrical geometry. Clearly, it is important for the molecular models of the reactant and product to accurately capture both their size and shape, as these determine the minimum or critical pore diameter, d = d_c, above which passing is possible. A reliable description of the shape is also essential to capture the ease of passing for slightly larger pores where the alignment of molecules with themselves and with the pore axis is necessary. A detailed model could represent all 16 atoms in PNB (3 O, 1 N, 7 C, 5 H) as overlapping spheres with the appropriate van der Waals radii and similarly could represent all 26 atoms in the aldol product (4 O, 1 N, 10 C, 11 H). See Figs. 2(a) and 2(c) for schematics of such detailed models. Sec. S2 of the supplementary material provides a complete listing of atom positions and sizes obtained from the NCBI PubChem database. The origins in Figs. 2(a) and 2(c) are the center of mass, and the axes correspond to the listed coordinates. However, at each step of the simulation with this detailed model, it would be necessary to check the overlap between 16 × 26 = 416 possible pairs of spheres with atom one in each of the species.

To reduce this computational cost, we instead utilize somewhat simplified “steric” models with fewer spheres but which still reliably capture the overall size and shape of the reactant and product species. For PNB, two C and one N are in the center of the molecule, and removal of the associated spheres does not significantly affect the overall molecular size or shape. Also, the five remaining peripheral C—H pairs can be replaced by a single sphere, which is shifted slightly outward from the actual C position to better capture the outer boundary of the overall PNB molecule. However, we retain the actual radius for C in order to capture the “thickness” of the PNB. The resulting planar 8-sphere molecule reasonably captures the size and shape of the 16-atom PNB species. See Fig. 2(b). For the aldol product, one end of this molecule including 2 O, 1 N, 6 C, and 4 H is the same as PNB, so again the central N and two C are removed, and each of the four additional peripheral C—H pairs are replaced by a single sphere. The remaining significantly out-of-plane 13 atoms (2 O, 4 C, 7 H) are represented by two larger suitably positioned spheres. The resulting 8-sphere model for the aldol product also reasonably captures its overall size and shape. See Fig. 2(d) and the supplementary material for more details. The origin and axes in Figs. 2(b) and 2(d) are preserved from Figs. 2(a) and 2(c).

The 8-sphere reactant and 8-sphere product steric models will be utilized in our Langevin simulations to determine the passing propensity. These require overlap checks for only 8 × 8 = 64 (versus 416) pairs of spheres after each move. Note, however, that these models have either low (planar) or no symmetry. This latter feature will be reflected in our analysis of diffusivity, which will also utilize further slightly simplified 8-sphere models in order to exploit an available hydrodynamic framework to determine diffusivity for general molecular shapes.³⁵ Finally, we note that a simpler and more computationally efficient approach to describe the molecular shape in the Langevin simulations might be to use general triaxial ellipsoids with sizes and shapes selected, attempting to best match those of the reactant and product.^36–38 A theory of diffusivity is available which shows that the diffusivity for such general ellipsoids can match that of the most (but not all) general 3D molecular shapes.³⁹ See Sec. S3 of the supplementary material. However, the higher fidelity of our 8-sphere models in matching molecular shape should justify the additional computational expense.

Finally, we remark that it will be useful to consider inscribing spheres surrounding the reactant (R) and product (P) species with centers at the center of mass and with radii R_R and R_P, respectively. Thus, their longest dimensions of the reactant and product are given by 2R_R ≈ 1.03 nm and 2R_P ≈ 1.39 nm, respectively, as noted in Sec. I. In addition, the shortest dimensions of the reactant and product species, which are in a direction orthogonal to the longest dimension, will be particularly significant for our analysis. These determine the critical diameter, d_c, of the cylindrical pore below which SFD is imposed due to steric blocking of passing of PNB and aldol product species. For our geometric models of these species, we estimate that d_c = d_c(eff) = 0.93 nm.

Non-trivial theoretical analysis is required for estimation of the complete 6 × 6 tensor, D, which describes both translational and rotational diffusion, and the coupling between them, for general shaped species immersed in a solvent.^40–42 However, it is instructive to note that some insight into the overall translational (and rotational) diffusion coefficients, D_t (and D_r), can be provided through approximate Stokes-Einstein-Debye (SED) relations which follow from a hydrodynamic treatment wherein the solvent is taken as a continuum viscous fluid. These SED diffusion coefficients have the form

with units of (distance)²/time for D_t and (radian)²/time for D_r, and where η denotes the fluid viscosity, and R_eff the effective hydrodynamic radius of the particle. These treatments are designed for macroscale objects immersed in a fluid. They should be used with caution on the molecular scale, although sometimes they are surprisingly effective.^43,44 Given the challenge of obtaining the full 6 × 6 diffusion tensor, we will utilize a hydrodynamic formulation,³⁵ noting that for our application only the relative (and not the absolute) magnitudes of the diffusion coefficients are relevant.

The 6 × 6 diffusion tensor, D, described below is determined with respect to a prescribed tracking point (TP) embedded in the reactant or product species, and also with respect to a set of body-fixed axes, the origin for which is the tracking point. Any tracking point can be utilized (including the center of mass), but the entries in the diffusion tensor depend somewhat on both its choice and on the choice of body-fixed axes. The diffusion tensor, D, and a corresponding related 6 × 6 grand resistance (or friction) tensor, Ξ ≡ k_BTD⁻¹, are naturally partitioned into 3 × 3 blocks, which correspond to translation (tt), rotation (rr), and translation-rotation coupling (tr and rt) as³⁵ 

In general, D_tr and D_rt are not symmetric. However, one has that D_rt = D_tr^T and Ξ_rt = Ξ_tr^T where T denotes the transpose, thus imposing symmetry on the overall 6 × 6 matrix. The overall translational diffusion coefficient is recovered from the 6 × 6 tensor via D_t = 1/3Tr(D_tt) and the overall rotational diffusion coefficient via³⁵ D_r = 1/3Tr(D_rr).

For some molecular shapes with sufficient symmetry, and always for 2D systems,^40–42 there exists a point called the center of hydrodynamic stress (CoH) for which translation-rotation coupling vanishes, i.e., D_tr = D_rt = 0.^35,45 However, such a CoH will not exist for the 3D shapes of relevance in this study. There only exists a center of diffusion for which D_tr and D_rt, are individually symmetric.^35,45 Finally, we remark that according to hydrodynamic theory, the diffusion tensor and specifically the above centers are not fixed for a molecular species in the presence of confinement but can vary according to the position and orientation relative to the boundaries.⁴⁵ However, we just use a fixed tracking point for simplicity.

The formulation which we use to determine the diffusion tensors is described in a review by Carrasco and Garcia de la Torre.³⁵ It is based on models for molecules either composed of groups of equal sized overlapping spheres⁴⁶ or multiple such groups which do not overlap each other and where spheres can have different sizes in different groups.⁴⁷ The steric models for PNB and the aldol product used for the Langevin simulations thus require further slight modification to apply this formulation. For PNB, we make a minor adjustment to the 8-sphere Langevin model making all spheres equal in size. For the aldol product, positions of the two larger spheres representing the 13 atoms at one end of the molecule are adjusted not to overlap the smaller spheres representing the remainder of the molecule. See Sec. S2 of the supplementary material. Specific results for D for PNB and for the aldol product will be presented in Sec. III A.

The motion of an individual reactant or product species is described by Brownian dynamics, i.e., overdamped Langevin dynamics, through a viscous isotropic homogeneous fluid representing the implicit solvent.^45,48 In addition, there are non-overlapping constraints with the other species (i.e., the aldol product cannot overlap the PNB reactant and vice versa) and with the walls of the cylindrical nanopore. To describe the coupled translational and rotational Brownian motion for a single species, it is convenient to introduce a 6-component state vector, Q, for each molecule which includes 3 components for the position, r, of a specified tracking point in the molecule (see below) and 3 orientation angles, θ. Components will be labeled i = 1-3 for translation and i = 4-6 for rotation. Then, the velocity vector d/dt Q includes 3 components specifying translational velocity, $v$⁠, and 3 components for angular velocity, ω. Overdamped motion is determined by a balance between (i) drag forces and torques described by a 6-component vector F^drag, and (ii) random forces and torques described by another 6-component vector, F^rand. The relationship between (a) the translational and angular velocities and (b) the drag forces and torques has the form

where the symmetric Ξ denotes the 6 × 6 grand resistance (or friction) tensor mentioned in Sec. II B, which can be obtained from the 6 × 6 diffusion tensor via Ξ ≡ k_BTD⁻¹. Note that M = Ξ⁻¹ corresponds to a grand mobility tensor. The random force, F^rand, has zero mean and has an amplitude determined by the fluctuation-dissipation relations so that

Our analysis of Langevin dynamics will use both a simplified description of diffusivity relative to real systems and also simplified numerical algorithms. In general, there are finite interactions between the diffusing reactant and product species and between these species and the wall. Also the diffusivity can depend on the location of the molecule relative to any boundary such as the pore wall. However, consistent with our model, we adopt a simplified treatment including only steric interactions. In addition, we ignore confinement effects and adopt for simplicity a configuration-independent diffusivity (and thus mobility) in the body-fixed frame corresponding to that of the isolated molecule in an infinite fluid. Given this configuration-independent body-fixed D tensor, we first determine its “square root,” μ, and then directly evaluate incremental motion in the body-fixed frame according to

and W is a column vector of normally distributed random numbers with zero mean and unit variance. At each time step, after incrementally propagating positions and orientations for both reactant and product, we check for an overlap and reject the move should such a configuration occur. Finally, we transform updated coordinates back to the space-fixed reference frame to facilitate tracking of the configuration of both species within the pore. One caveat is that this simple and efficient numerical algorithm does not precisely account for stochastic drift terms,^45,48 but we anticipate that these terms and their treatment will not significantly affect resulting passing propensity values. Further technical details including the numerical integration protocol are given in Sec. S4 of the supplementary material.

Langevin dynamics simulations for a PNB reactant and an aldol product pair, incorporating diffusion tensors as prescribed in Sec. II C, are used to assess a suitably defined passing propensity, $P$⁠. We establish a space-fixed Cartesian axis system with the z-axis along the pore axis and the (x, y)-axes in the orthogonal directions. Our strategy is to formulate the definition and determination of $P$ in the Langevin simulations so as to be consistent with the spatially discrete coarse-grained (CG) model described in Sec. II F to which this parameter is “passed.” The cell width, w, in this spatially discrete CG model will be chosen based on the characteristic linear dimensions of the reactant and product species. Specifically, we assign w = R_R + R_P, where again R_R and R_P are the radii of the inscribing spheres for the reactant and product species, respectively. Thus, for the reactant and product located in adjacent cells, the difference, Δz, of z-coordinates of the centers is Δz = +w. Passing corresponding to exchange leads to Δz = −w and separation where either species hops to neighboring empty cells leads to Δz = +2w.

Thus, to match these scenarios in our Langevin simulations, we start with PNB on the left, say, and the aldol product on the right with the Δz = Δz₀ = +w. Thus, in this initial configuration, the reactant and product species do not interact with each other (i.e., they cannot overlap), and they are chosen so that their other coordinates (x, y, and orientations) are randomly distributed within the constraints of not overlapping the pore wall. The simulation is run until either Δz reaches +2w, which is designated as the separation of the reactant and product, or Δz reaches −w, which is designated as passing. See Fig. 3. The passing propensity is extracted from running typically O(10⁵) simulations and defining $P$ = P_pass as the fraction of these in which passing (rather than separation) is achieved. If P_sep is the fraction of outcomes corresponding to separation, then P_sep + P_pass = 1. This process is repeated for various pore diameters, d = d_p(eff), to obtain $P$ = P(d) as a function of pore diameter, d. We will present numerical results for P(d) versus d in Sec. III B.

By definition, $P$(d) → 0, as d → d_c (the critical diameter for SFD). The non-trivial scaling behavior for d close to d_c will also be discussed in Sec. III B. In the regime of large d, it follows that there is negligible interference between the diffusion of PNB and the aldol product. As a result, determination of the passing propensity reduces to a simpler process corresponding to a first-passage problem for a one-dimensional random walk (RW) with traps.⁴⁹ Specifically, the RW in Δz starts at Δz = Δz₀ = +w and evolves in the interval −w ≤ Δz ≤ +2w until reaching a “trap” either at Δz = +2w (corresponding to separation) or at Δz = −w (corresponding to passing). Since the change in Δz for separation is one half of that for passing, it follows from RW theory that P_sep:P_pass = 2:1. Consequently, one has that $P$(d) → 1/3, as d → ∞.

Many-particle MD analysis of the solution-phase reaction system with an explicit solvent, or even many-particle Langevin dynamics with an implicit treatment of solvents, cannot access the appropriate time scales for the overall catalytic reaction-diffusion process of interest. Thus, we implement a spatially discrete coarse-grained description of the process where each pore is divided into a linear array of cells.¹⁰ The width of the cells matches to the characteristic linear dimension of reactant and product species. Specifically, we assign w = R_R + R_P as noted in Sec. II E so that w ∼ 1 nm. The cell volume is thus somewhat above 1 nm³. Each cell can contain at most one PNB reactant (R) or aldol product (P) species, where other solvent species are treated implicitly. Again, these solvent species include both hexane and the reactant acetone which is present in excess. Diffusion of PNB and aldol product species within the pores is described by hopping to adjacent unoccupied cells at rates h = h_R and h_P, respectively. We emphasize that here unoccupied means that cells are not populated by either PNB or aldol product, not that cells are free of solvents. Diffusion also occurs by exchange of PNB and the aldol product species between adjacent populated cells at rate h_ex = ½(h_R + h_P)P_ex. Thus, P_ex serves as a measure of the passing propensity. Furthermore, P_ex = 0 corresponds to SFD since, without exchange, the prescribed hopping dynamics imposes a single-file constraint on PNB and the aldol product. We assign all cells as catalytically active and prescribe that the reaction of PNB to the aldol product occurs in each cell populated by PNB at a microscopic rate k per cell. One can readily refine the model to incorporate, e.g., a reversible reaction for any specified equilibrium constant. See Fig. 4 for a schematic of this spatially-discrete stochastic CG model.

The results for $P$ values from the appropriately crafted Langevin simulations in Sec. II E are used determine the parameter P_ex. To clarify the relationship between these quantities, consider the coarse-grained model with a pair of isolated PNB reactant and aldol product species in adjacent cells with PNB on the left, say, and the aldol product on the right so that Δz = +w in this configuration. The following distinct events can occur starting from this initial configuration: (i) PNB and the aldol product can exchange with rate ½(h_R + h_P)P_ex corresponding to passing and resulting in Δz = −w; (ii) the PNB can hop left at rate h_R, or the aldol product can hop right at rate h_P, either event corresponding to separation and resulting in Δz = +2w. Consequently, the passing propensity is given by

As a result, given $P$ obtained from the Langevin simulations, the corresponding P_ex can be determined from P_ex = 2$P$/(1 − $P$⁠). This P_ex can vary from 0 to 1.

It is also necessary to consider adsorption and desorption of PNB reactants and aldol products at the pore openings.¹⁰ Of course, at the onset of the reaction with negligible conversion, adsorption of product will also be negligible. The adsorption rate will control the concentration of PNB reactant and aldol product in the pore, and it will be selected in the simulations to achieve appropriate concentrations.

Precise analysis of this spatially discrete stochastic model is achieved by Kinetic Monte Carlo (KMC) simulation wherein processes are implemented with probabilities proportional to their physical rates.¹⁰ For any specified conversion, α, of reactant to product in the external fluid, the concentration distribution of the reactant and product within the pore quickly achieves a quasi-steady state from which the reactivity for that conversion can be determined. There are substantial fluctuations in concentration within a single pore. However, mean behavior of concentration distributions at a specific α, and thus of reactivity, is determined precisely by averaging over many simulations, analogous to the automatic averaging over many pores in MSN in the actual reaction system. Alternatively, one can average a single simulation over time for fixed α to obtain precise mean values. Such simulations reveal negligible penetration of reactant into the pore for SFD as noted in Sec. I,^2,5 (except for extremely small values of k/h). Reactant penetration naturally increases with increasing P_ex.¹⁰ Also, for this model,⁵⁰ one finds that reactivity decreases linearly to zero with increasing conversion, α. Thus, analysis of reactivity at the onset of reaction (i.e., for negligible conversion) can provide a complete picture of kinetics.

Diffusion tensors are determined using the formulation of Ref. 35 for the models of PNB and the aldol product described in Sec. II B. The body-fixed axes for PNB and aldol product are shown in Fig. 2. For PNB, (x, y)-axes are in the plane of the molecule with the x-axis in the direction of the longest dimension, and the z-axis is out-of-plane. For the aldol product, one end of which has the same structure as PNB, the axes are chosen analogously to PNB with the x-axis in the direction of the longest dimension. As noted in Sec. II E, any body-fixed location can be used as the tracking point in simulations. This includes the center of mass which however is not intrinsically related to diffusion behavior. Thus, we use a slightly shifted tracking point (TP) for the calculation of diffusion tensors which is closer to the center of diffusion (cf. Ref. 35). However, this slight change in TP makes little difference to even the dominant diagonal entries in the diffusion tensors. See Sec. S3 of the supplementary material. For the 6 × 6 diffusion tensors, D, components i = 1, 2, and 3 corresponding to translational degrees of freedom are in the direction of x-, y-, and z-axes, respectively. Components i = 4, 5, and 6 correspond to the rotational degrees of freedom, which specifically correspond to rotations about the x-, y-, and z-axes, respectively. Distances are measured in angstroms. In our analysis, the solvent fluid viscosity, which implicitly incorporates the time scale, is set to unity. Since our simulation analysis of and results for passing propensity depend only on the relative magnitude of diffusion coefficients, the absolute values are not required.

For PNB, results for the 3 × 3 blocks of D are

so that D_t = 1.61 × 10⁻² and D_r = 5.14 × 10⁻⁴ (in the units specified above). Given the near planar molecular structure, the coupling in D_tt between the x-component or y-component and the out-of-plane z-component is very weak. The coupling between the x- and y-coordinates is more than two orders of magnitude stronger, although D_tt is still near-diagonal. The magnitude of the diagonal entries reflects the cross section orthogonal to the axis, smaller cross sections corresponding to larger translational diffusion coefficients. Similarly, D_rr is also near-diagonal, with the largest entry corresponding to rotation about the x-axis which presents the smallest cross section. The translation-rotation coupling is quite weak.

For the aldol product, results for the corresponding 3 × 3 blocks of D are

so that D_t = 1.27 × 10⁻² and D_r = 2.68 × 10⁻⁴ (in the units specified above). Noting the significant out-of-plane molecular structure, the coupling in D_tt between y- and z-components is no longer weak (and is actually now the strongest), although D_tt is still near-diagonal. Similarly, D_rr is also near-diagonal, with the largest entry corresponding to rotation about the x-axis. The translation-rotation coupling is again quite weak. Results for the corresponding μ tensors for both PNB and the aldol product which are used in the Langevin simulations are given in Sec. S3 of the supplementary material.

Finally, one can estimate effective hydrodynamic radii within a Stokes-Einstein-Debye framework from the relation R_eff = [(3/4)D_t/D_r]^1/2 (in angstroms) following from (1). This yields values of R_eff = 4.8 Å for PNB and R_eff = 6.0 Å for the aldol product. These rough estimates are in good agreement with the geometric size of the molecules and appropriately reflect the larger dimensions of the aldol product. Also, this type of crude analysis does provide insight into the relative magnitude of diagonal entries in D_tt and D_rr.

Computationally demanding Langevin simulations were performed to precisely determine the passing propensity, $P$⁠, for four choices of effective pore diameters, d = d_p(eff) = 1.5, 2.0, 2.5, and 3.0 nm (taken as the diameter of the cylindrical pore in the simulations). These choices effectively cover the experimentally relevant range.

For narrower pores closer to the critical diameter of d_c = d_c(eff) = 0.93 nm for SFD, analytic theory will be used to extrapolate the simulation results and similarly for large d. We run a total number N_trial of Langevin simulation trials for a PNB and the aldol product pair with the difference Δz of z-coordinates of tracking points initially separated by Δz₀ = +w (as described in Sec. II F). In each trial, the other coordinates are assigned different initial values randomly sampling the phase space of allowed values where molecules do not overlap the pore walls. The simulation is run until either Δz reaches +2w (separation) or −w (passing). $P$ = P_pass is the fraction of these trials where passing is achieved. The number, N_trial, of trials required for accurate determination of $P$ depends on the number of passing events N_pass ≈ $P$ N_trial, with uncertainty ±√N_pass. Thus, e.g., for $P$ = 0.025, selecting N_trial = 10⁵ yields an uncertainty in $P$ of ±2%. We select N_trial = 2 × 10⁵ since $P$ for the smallest simulated d = 1.5 nm is close to this value. See Sec. S5 of the supplementary material for further discussion of uncertainties. Precise values for $P$ obtained from using such large N_trial will also facilitate effective integration of these simulation results with analytic theories, as described below.

Another important technical issue for precise determination of $P$ is the selection of a sufficiently small time step, Δt, for accurate numerical integration of the Langevin equations (5). Results presented in Table II indicate reasonable convergence when the time step is decreased to below Δt ≈ 0.1 (with units corresponding to unity viscosity) so that D_tΔt is below about 10⁻³ (cf. Ref. 34). Next, we provide a more comprehensive picture of the variation of $P$(d) versus d by combining the simulation results for Δt = 0.01 with analytic insights into behavior for d close to d_c and for large d.

For small values of gap size, g, defined as g = d − d_c, accurate determination of the small values of $P$ from Langevin simulation becomes computationally expensive. Extremely large N_trial is needed to obtain sufficiently large N_pass ≈ $P$ N_trial required to reduce uncertainty. An effective alternative is to exploit insights from previous analyses of passing processes for simpler shaped molecules which utilized simple and approximate Transition State Theory (TST)¹² and also an exact Fokker-Planck equation formalism.³⁴ TST determines the (entropic) free energy barrier for passing as δF ∼ −k_BTln(Ω_min), where Ω_min denotes the 10-dimensional phase space volume for the two molecular species with Δz fixed at the transition state (TS) for passing. Ω_min at this “crowded” TS is smaller than the phase space volume, Ω, for other Δz. Then, one has $P$_TST ∼ exp[−δF/(k_BT)] ∼ Ω_min. Analysis of the system geometry at the TS indicates that $P$_TST ∼ Ω_min ∼ $gσTST$⁠, as g → 0, where σ_TST reflects the number of degrees of freedom at the TS.²⁴ One finds that σ_TST = 2.5 for two spherical molecules, 4.5 for two spheroids, and 8.5 for two general shaped molecules. However, analysis of passing for simple systems based on the exact Fokker-Planck equations, which are equivalent to the Langevin equations, indicates that the exact passing probability satisfies

where σ is significantly below σ_TST (e.g., σ = 1.7 rather than σ_TST = 2.5 for two spheres).³⁴ Thus, we adopt the form (9) for the PNB + aldol product system. The parameters A = 0.272 and σ = 4 are selected to match both the value of $P$ and its rate of change with d, at d = 1.5 nm, as determined by the Langevin simulation results. This σ-value appears consistent with expectations from the above-mentioned comparison of TST and Fokker-Planck analyses.

In Sec. II F, we have noted that $P$(d) → 1/3, as d → ∞. For large d, the reduction in $P$ from its asymptotic value should scale like 1/d, so we adopt the approximate form $P$(d) ≈ 1/3 − B/(d − d^*) and select the parameters B = 0.137 and d^* = 1.5 to match both the value of $P$ and its rate of change with d, at d = 3.0 nm, as determined from the Langevin simulation results. An interpolation of simulation results, combined with the analytic forms for $P$(d) for larger and small d, is shown in Fig. 5.

In our stochastic modeling of the overall reaction-diffusion process as described in Sec. II B, we select a pore of 200 cells labeled n = 1–200 corresponding to a length of L ∼ 200 nm, compatible with the typical MSN diameter. The hop rates for reactant (and product) species, h = h_R (and h_P), to adjacent unoccupied cells are related to D_t = 1/3Tr(D_tt) via h = w⁻²D_t. Thus, based on the results in Sec. III A, we have h_R:h_P = 1.611:1.269. We will assign h_R = 1 and h_P = 0.788, which sets a time scale for the coarse-grained model. Actually, this time scale is not relevant as we consider only quasi-steady-state behavior. The key exchange probability, P_ex = P_ex(d), characterizing the propensity for passing of reactants and products for various effective pore diameters, d, is determined from $P$_ex(d) = 2P(d)/[1 − $P$(d)] with $P$(d) given in Sec. III B. The results presented below consider a reaction with microscopic reaction rate per cell, k, for conversion of reactant to product. To recover experimental behavior, we find that k must be selected to be well below the hop rates. We anticipate that the reaction rate per cell should reflect the total number of catalytic sites within that cell which in turn reflects the pore perimeter length. Consequently, we set k = k₀d with d in nm, and will adjust k₀, or equivalently k₀/h (with h = h_R = 1) to determine if the model can mimic experimental behavior. Finally, we must also specify the rate of adsorption and desorption at the pore openings. Rates for desorption from the end cell of the pore are selected to be the same as for diffusion within the pore. Our analysis focuses on the assessment of reactivity at the onset of reaction when conversion of reactant to product in the external fluid is negligible. The reactant adsorption rate is selected to ensure that the reactant concentration per cell inside the pore would equal a target concentration of X_p per cell, which as noted in Sec. II A should significantly exceed the reactant concentration in the bulk corresponding to X_b ≈ 0.1.

In this model, we let ⟨R_n⟩ denote the mean probability or “concentration” of PNB reactant in cell n and ⟨P_n⟩ the mean concentration of aldol type product. Then ⟨X_n⟩ = ⟨R_n⟩ + ⟨P_n⟩ is the total (reactant + product) concentration per cell. This concentration is not exactly uniform,¹⁰ and thus adsorption rates are chosen so that ⟨X_n⟩ is close to the target concentration X_p near the pore openings. These concentration profiles can be extracted precisely from KMC simulation. Then, the total reactivity, $R$⁠, is simply obtained from the total amount of reactant in the pore multiplied by the irreversible reaction rate, k, as $R$ = k∑_n ⟨$Rn$⟩.

A key goal is to assess whether the experimental observation of a dramatic increase in yield upon increasing of d from around 1.3 nm to 2.1 nm can be achieved in our model. The variation of reactivity with pore diameter (and other parameters) is non-trivial and not described by standard mean-field (MF) reaction-diffusion equations. However, a qualitative assessment can be provided by generalized hydrodynamic (GH) formulation of reaction-diffusion kinetics in these systems.¹¹ See also Sec. S7 of the supplementary material. This GH formulation indicates that the penetration depth, L_p, of the reactant into the pore scales like L_p(SFD) ∼ (h/k)^1/4 for SFD is in contrast to conventional MF reaction-diffusion behavior L_p(MF) ∼ (h/k)^1/2 which would apply for wide pores. Reactivity measured here at the onset of reaction scales like $R0$ ∼ 2kL_p, accounting for reactant penetration from both ends of the pore. Thus, the mean-field enhancement of initial reactivity relative to its SFD should scale like

Consequently, the desired large target values require sufficiently large h/k, i.e., sufficiently small k, given our choice h = O(1). One caveat is that the above results for L_p only apply for pore length L > 2 L_p since otherwise L_p saturates at ½L. As a result, $R0$ ∼ kL would be similar for SFD and for wide pores. Thus, while k must be selected to be sufficiently small to achieve large h/k and large $E0$(MF), it cannot be selected to be too small.

Results of KMC simulations are shown in Fig. 6 for the actual enhancement of initial reactivity relative to the SFD value, $E0$(d) = $R0$(d)/$R0$(SFD) ≥ 1, where $R0$(SFD) corresponds to the value of $R0$ for d ≤ 0.93 nm. Results are shown for selected pore concentrations, X_p, which are significantly above the external bulk fluid concentration of X_b ≈ 0.1 and for a broad range of k₀. It is clear that there are two requirements to achieve the experimentally observed strong increase in yield and thus reactivity, upon increasing d from 1.3 nm to 2.1 nm. The first is that the reaction rate should be selected close to an optimum value of around k₀ ≈ 0.001. For larger k₀ (and thus smaller h/k), the variation of reactant penetration depth, L_p, with d is insufficient to produce the large variation in reactivity. On the other hand, for smaller k₀, one has that L_p ∼ ½L, so again reactivity displays a weak dependence on d. The second requirement based on the results in Fig. 5 is that the reactant concentration interior to the pore must be sufficiently high. Results for other values of X_p shown in Sec. S7 of the supplementary material fit the trend illustrated in Fig. 5 with the normalized reactivity, $E0$(d), for k₀ = 0.001 increasing smoothly as a function of X_p. Behavior for X_p above about 0.4 appears consistent with the experiment. A comparison of the variation of concentration profiles for SFD (d ≈ 0.93 nm) and for d = 1.3 nm and 2.5 nm is shown in Fig. 7 for the case X_p = 0.40. There is a strong enhancement in the penetration of the reactant into the pore with increasing d.

Since only the reactant concentration in the external fluid outside the pore is readily determined, it is not straightforward to confirm our prediction regarding the internal concentration. However, we have provided two reasons in Sec. II A as why significant enhancement is expected. In principle, MD simulation of the complex pore + internal and external fluid system with realistic potentials could provide insight into this issue. With regard to our prediction of low k₀, we are not aware of previous estimates. However, conventional assessment would be based on traditional mean-field analysis of measured kinetics which is not appropriate for this system with inhibited passing. Our modeling framework in principle provides a more reliable way to extract reaction rates.

The above analysis considers only the reactivity, $R0$(d), at the onset of reaction versus d rather than the reaction yield measured in terms of conversion, α_f(d), after t_f = 2 h as measured in the experiment. To connect the above results to the reaction yield, we exploit the result for our model that the reactivity, $R$(α), for conversion α should have the form¹⁴ $R$(α) = $R0$(d) (1 − α). Thus, one has d/dt α = c $R$(α) which can be integrated to obtain the yield, α_f(d), at time t = t_f which is given by α_f(d) = 1 − exp[−K_f $E$(d)], where K_f = c $R0$(SFD) t_f. Thus, the enhancement in yield upon increasing d above the SFD regime satisfies

where $E0$(d) = $R$₀(d)/$R$₀(SFD) ≥ 1 is the enhancement in initial reactivity, as given above. The enhancement in yield is somewhat below $E0$(d) but still substantial, consistent with experiment. See Sec. S8 of the supplementary material for further discussion.

The above analysis successfully demonstrates that the substantial enhancement in yield for PNB conversion to aldol for a modest increase in pore diameters in amine-functionalized MSN can be reasonably tied to a strong variation in the passing propensity for reactant and product species over this range of diameters (the narrowest of which is close to the onset of SFD). Certainly, we cannot completely rule out other factors such as dependence on the pore diameter of the activation barrier for reaction or on different ordering in the solution-phase. However, we expect these dependencies to be relatively weak. We emphasize that this connection between yield and passing propensity can only be made by the utilization of an appropriate multiscale modeling framework.

Passing of reactant and product species within narrow pores requires orientational alignment. Consequently, accurate assessment of passing propensity, $P$⁠, requires the development of a Langevin formulation incorporating appropriate treatment of the associated translational and rotational diffusion. Precise simulation results for $P$ facilitate integration of insights from an equivalent Fokker-Planck formulation. This in turn provides a comprehensive picture of passing propensity including behavior in the regime of narrow pores close to the threshold for SFD. In addition, consistent formulation P-values and inputs into coarse-grained models are required to treat the overall catalytic reaction-diffusion process on the appropriate time scales and length-scales. This ultimately enables reliable assessment of the consequences of the increase in $P$ for the enhancement in experimentally measured reactivity and yield.

Certainly, our treatment has included simplifications regarding molecule and pore geometry and their rigidity, a hydrodynamic treatment of diffusivity, and maximally-coarse spatial coarse-graining (where the latter can be straightforwardly refined⁵¹). However, we believe that our basic insights are reliable. This work constitutes the first such effort integrating various non-trivial and diverse theoretical approaches to achieve system-specific multiscale modeling address basic issues for catalysis in nanoporous materials.

See supplementary material for more details of the following: experimental procedures (S1); molecular models (S2); theory of diffusivity (S3); Brownian dynamics (S4); Langevin analysis of passing propensity (S5); diffusive dynamics in the CG model (S6); generalized hydrodynamic formulations (S7); coarse-grained stochastic simulation results (S8).

We thank David Ackerman for discussions on Langevin simulations and Chi-Jen Wang for insights from Fokker-Plank equation analysis. This work was supported by the U.S. Department of Energy (USDOE), Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences through the Ames Laboratory Chemical Physics program (A.G. and J.W.E.), and the Ames Laboratory Homogeneous and Interfacial Catalysis program (I.I.S.). The work was performed at Ames Laboratory which is operated for the USDOE by Iowa State University under Contract No. DE-AC02-07CH11358.