Confined Ru‐catalysts in a Two‐phase Heptane/Ionic Liquid Solution: Modeling Aspects

A modeling approach for atomic‐resolution studies of sup‐ ported ionic liquid phase (SILP) catalytic systems in silica mesoporous confinement with surface hydroxyl and functional groups is proposed. First, a force field for the Ru‐based catalyst is developed. Second, its solvation behavior within a bulk two‐phase system of heptane and an IL is studied. Third, static and dynamic properties of the confined system are investigated. Using classical molecular dynamics simulations, experimentally inaccessible properties can thus be studied that are important for an optimization of a SILP system for performing a ring‐closing metathesis reaction.


Introduction
Recent advances in designing mesoporous silica materials with surface-anchored functional groups have shown a promising potential for applications with respect to the activity and selectivity in molecular heterogeneous catalysis. [1][2][3] For example, by immobilizing a Rh-complex on SBA-15 particles, a rate enhancement for 1-octene hydroformylation was observed compared to the homogeneous analogue. This finding can be explained by the suppression of the formation of inactive forms of the catalyst inside the pores. [4] Other examples include the improved selectivity for olefin or cyclooctane metathesis reactions carried out in mesoporous silica materials [5,6] or synergetic effects between multifunctional groups anchored in close proximity in confined spaces. [7] Selectivity and conversion of reactions can be affected due to controlling the factors that influence the spatial distribution and the mobility of reactants and products including molecular size, shape, configuration, degree of confinement, pore topology, strength of adsorption on pore walls, and the possibility of hydrogen bonding between reactants or/and products.
Instrumental for the application and rational design of functionalized mesoporous materials is the existence of a model for predicting how pore dimension and surface functionalization influence the properties of the fluid under confinement. [8] Fluids confined in mesoporous materials present a particular challenge to empirical tuning, as it is difficult to experimentally probe properties such as solvent composition and phase behavior within the mesopores. Mass transport of molecular compounds through porous solids is a decisive step in molecular heterogeneous catalysis. It is a multi-scale, hierarchical phenomenon: the effective diffusion through the macromesoporous material is influenced by parameters such as grain boundaries and particle packing on the macropore scale (> μm), as well as by factors such as particle size and connectivity of pores on the mesopore scale (> 10 nm, < μm). More importantly, meso-scale diffusion and macro-scale diffusion are first and foremost determined directly by processes on the molecular scale (< 10 nm), which depend on numerous factors like pore-size, pore-size distribution, interactions of the diffusing species with the solid surfaces and with the solvent. [9] For reacting systems, the transport of the reactants to and products from the catalytically active complex, as well as the interaction of the various species with the solvent and the functional groups on the surface is decisive. Therefore, there is a growing interest in computational modeling studies of relevant systems in order to provide a fundamental molecular-level description of confinement effects. [10][11][12][13][14][15][16][17][18] The molecular insight gained by such studies was used, for example, to obtain a unified description of diffusion inside mesoporous and microporous structures [20] or to rationalize the effect of surface chemistry on olefin metathesis in confined geometries. [21] An important group of industrial catalysts are supported liquid phase (SLP) catalysts. [22][23][24][25][26][27][28][29][30][31] In this context, MD simulations revealed for example fundamental insight into reactive ionic liquid films for supported ionic liquid phase (SILP) catalysis, relevant for electrochemical systems or the water-gas shift reaction. [32][33][34] Atomistic MD simulations may also be used to generate parameters for coarse-grained models that span significantly larger length and time scales than models with atomistic resolution. [35][36][37] It is therefore desirable to further widen the scope of such simulations by increasing the complexity of the studied systems. Within the context outlined above, we study a linker-free incorporation of a Ru-based catalyst inside a functionalized mesoporous environment filled with a two-phase fluid phase consisting of heptane and an ionic liquid.

Simulation Details
In this work, we study the microstructuring and the properties of IL/heptane mixtures in confinement by placing the solution within a nanometer-sized pore. To this end, we perform atomistic Molecular Dynamics (MD) simulations of different model systems. We begin with a two-phase system made of an ionic liquid and heptane. For the former, we consider 1-Butyl-3methylimidazolium trifluoromethanesulfonate [BMIM] + [Otf] À . We then place a Ru-catalyst within this two-phase bulk solution and monitor its diffusion. At a next step, we consider only a bulk IL solution with the Ru-catalyst in order to investigate configurational aspects of the catalysts in the solution. Having gained an insight from the behavior in bulk, we place the Rucatalysts, as well as substrates within a nanometer-sized pore and follow the dynamics of the system. The surface of the pore is covered by silanol groups, resembling a fully hydroxylated surface. [38] The inner walls of the pore have been functionalized with imidazolium (Im) molecules, in order to tune the polarity of the pore and its interaction with the IL solution. In the solution within the pore, chloride ions were added to play the role of counterions and neutralize the surface functional groups inside of the pore. The corresponding structures of all molecules considered in our work are shown in Figure 1.
For our atomistic MD simulations with the GROMACS 5.1.3 software package, [39][40][41][42] we have used the OPLS/AA force field for the ions, heptane, and the substrate. [43][44][45][46] We discuss below our own developed force field for the Ru-catalyst. We employed the reduced charged model (× 0.8) for the IL. Therefore, the partial charges of other ionic species including the surface functional group and chloride atoms are also uniformly reduced by the factor of 0.8. [16,45] This approach is commonly used to study bulk IL solutions. However, alternative approaches do exist (refer to [47,48] ). In order to resemble typical experimental conditions [19] all simulations were conducted at T = 353 K. The temperature was controlled by an improved velocity-rescaling thermostat [49] using a coupling time constant of 0.1 ps. Electrostatic interactions were treated through the particle mesh Ewald (PME) method, [50,51] where a real-space cutoff of 1.0 nm and a grid spacing of 0.16 nm with a fourth-order interpolation scheme were used. The Lennard-Jones (LJ) interactions were truncated at 1.0 nm and shifted to zero. The equations of motion were integrated using the leapfrog algorithm [52] with a time step of 2 fs. All bonds with hydrogen atoms were constrained by the LINCS algorithm. [53] An energy minimization was first performed using a conjugate-gradient method, followed by an equilibration time of 20 ns under the (NVT) ensemble. For the confined system, the simulation was extended for 100 ns under the (NVT) ensemble. The simulations of the bulk solutions, were performed under the (NPT) ensemble with the pressure kept constant at p = 1 bar by an isotropic Parrinello-Rahman barostat [54] with a coupling time constant of 2 ps and a compressibility of 4.5 × 10 À 5 bar À 1 . The total time of the NPT simulation was 50 ns. We have used computational boxes, which are periodic in the x, y, z-directions. In these boxes, apart from the fixed pore geometry, all other molecular species were randomly placed. The exact details and geometries of the computational boxes are given in Table 2 as well as in Figure 3. For the random insertion of molecules, the software package PACKMOL [55] for the bulk simulations was used. The corresponding numbers of molecules and ions are shown in Table 1 for each type of simulations.

Force-field for the Ru-catalyst
The OPLS/AA force field used for the IL, the heptane, and the substrates does not include a parameterization for the Rucatalyst. We have thus developed the respective parameters for the catalyst based on quantum-mechanical calculations within the density functional theory as implemented in the ORCA quantum chemistry code. [56] Specifically, the parameters for the catalyst were generated using the B3LYP [57] exchange-correla-  [21] and (e) the Rucatalyst. [19] tion functional with the D3BJ dispersion correction [58] and the def2-TZVP basis set. [59] For the self-consistent field (SCF) calculation the Pulay method [60] was utilized. The partial charges on the atoms were calculated through the CHELPG program. [61,62] Under these methods, we have relaxed the geometry of the catalyst and have used the equilibrium structural properties for the force-field parameterization. All the partial charges on each atom were in the end multiplied by the rescaling factor of 0.8 to match the reduced charge force field for ILs. We employed a semi-rigid model in which all bonds and angles are constrained by the same harmonic potential respectively. The optimized structure with the labels of each atom in the Ru-catalyst, as well as exact parameters of bond length, angles, partial charges are given in the supplementary information.

Bulk solutions
The purpose of the bulk phase simulations is to demonstrate that the catalyst is solvated in the IL phase and that heptane and the IL form a stable two-phase system with the catalyst residing in the IL phase. The system compositions are summarized in Table 1. Both the catalytic molecules and the IL ions are randomly inserted into the 3D-periodic box. A snapshot of the simulations for the two types of bulk solutions are depicted in Figure 2.

Confined catalyst within a Nanopore
The pore model was created utilizing the PoreMS Python package [63] Version 0.2.0. [64] The cavity pore system consists of two capsule-shaped cavities at a distance of 2 nm carved out of a β-cristobalite block. The cavities were constructed by combining a cylinder with half a sphere of the same diameter of 5 nm. Two bulk reservoirs were attached on each of the exterior surfaces, representing the (111) face of β-cristobalite. The interior pore surfaces have been functionalized with Im molecules with a coverage of 1.66 μmol m À 2 and a residual hydroxilation of 7.75 μmol m À 2 . The exterior surface is left unchanged with a hydroxilation degree of 8.79 μmol m À 2 . Figure 3 depicts the constructed pore system, while the structural details of this system are listed in Table 2.

Catalysts in the bulk solutions
For the analysis of the results, we will begin with the bulk solutions and then move on to the confined (solution, catalysts, substrates) system in the nanopore. For the bulk solutions, we first discuss the results related to the catalysts in the two-phase IL/heptane solution and then the catalysts in the bulk IL solution. For the former, our results underline that the heptane and the IL remain in two separated phases and do not mix. The catalysts remain throughout the simulation time within the IL part of the two-phase solution. In order to quantify these results, we calculate the mass density of all the species, heptane, IL cation, IL anion, catalyst in the two-phase bulk solution. For the calculation of the mass density, the simulation box is divided into slices and the average density of each slice is calculated. The results are depicted in Figure 4(a) and are shown along the z-axis, which is perpendicular to the interface between the IL and the heptane phase. It can be clearly observed that the heptane and IL do not mix and form two  phases. Some small fluctuations in the mass density can only be observed at the interface. This is expected, as there is no clear line that indicates the phase separation. Furthermore, the peak in the mass density of the catalysts lie well within the IL solution. Note, that the interface between the two phases is around 7.5 nm. Accordingly, we can infer that the catalysts are dissolved within the IL and not within the heptane phase. This fact is expected as the catalysts are divalent and strongly interact through electrostatic forces with the ionic species of the IL showing an attraction towards the anions of the IL solution.
In order to investigate the diffusion properties of all species in the bulk solutions, we have monitored the time evolution of the mean squared displacement (MSD) of all species. This quantity provides a very good estimation of the mobility of each species. The MSD for the two-phase bulk system can be seen in Figure 4(b). Inspection of this figure, clearly shows that the heptane molecules diffuse faster than the other ionic species (IL cations, IL anions, catalysts). This is physically intuitive due to the fact that there is less electrostatic interaction between the neutral species such as the heptane molecules. The IL cations and IL anions show a very similar diffusion, with former moving a bit faster within the IL phase. This observation is in agreement with previous studies, investigating the various factors that control diffusion in ionic liquids. [65][66][67] On the other hand, the Ru-catalysts diffuse slower than the IL species due to their stronger electrostatic interaction with the surrounding IL ions. Note, that the not so smooth MSD curve for the catalysts arises from the small number of the catalytic molecules in the solution and the corresponding poor statistics in gathering the respective data, as compared to the other species.
In order to study the solvation of the catalysts within the IL phase and address issues of clustering and aggregation, we have performed the bulk IL simulations in which we have placed the catalysts. For this solution, the center of mass radial distribution function (RDF) is calculated between the catalyst and all the species (IL cation, IL anion, catalyst). As revealed from Figure 5(a), the anions form the first solvation shell around the catalyst due to the strong electrostatic interaction between these two types of molecular species in the solution. The cations of the IL solution follow in the second solvation shell. The catalysts accumulate around the second solvation shell, which is mostly occupied by the cations. These trends can be explained by the strong electrostatic repulsion between the divalent catalysts and the less strong repulsion between the catalysts and the IL cations. In this way, the order of the peaks in the RDF for the reference catalyst closely follows -as expected -the type of interactions between the different species: the first shell corresponds to the anions (attraction to catalysts), the second shell to the cations (repulsion to catalysts, but attraction to IL anions), and the third shell to the catalysts (strong repulsion to catalysts, smaller repulsion to IL cations). Note, that the two first solvation shells formed by the IL anions and the IL cations, respectively screen the charge of the reference catalysts, so that its repulsion with the catalysts at the third solvation shell is negligible.
The MSD dynamics of the IL and the catalytic species in Figure 5(b) follow the trends of the two-phase system in Figure 4(b). A small decrease in the MSD can be observed in the case of the IL species as compared to the two-phase solution. The reason could be assigned to the fact that the IL species have smaller electrostatic interactions with the heptane molecules at the interface than with their own IL species in the bulk. This results in a slightly higher diffusion in the two-phase system close to the interface. The same decrease in the MSD and the diffusion can be seen also in the case of the catalysts. Again, this arises from the fact that catalytic molecules close to the interface in the two-phase solution interact less strongly with the heptane molecules at the interface, thus moving faster close to the interface. Note, that due to the increased number of catalytic molecules, the statistics in gathering the data for the MSD is better, thus the curve shows the expected linear behavior at longer times.

Confined catalysts in the Nanopore
We further move to the main system of interest, namely the confined catalysts in the pore. In this system, substrates are also randomly placed within the heptane solution. We calculate the mass density profile of all species along the axis perpendicular to the pore surface. Initially, the pore was filled only with the IL phase and the heptane phase covered only the part in front of the pore opening. During the simulation, the heptane molecules do enter the pore and move towards the bottom of the pore close to the catalyst as the black curve in Figure 6(a) underlines. However, the heptane phase does not cover the catalysts, which still remain within the IL phase. Nevertheless, the partial filling of the pore with the heptane solution does allow for the substrates to diffuse into the pore and reach the catalysts in order to eventually interact with the catalytic centers. The results show that the substrates accumulate further away from the catalyst, though in the pore. Note, that these results are time-averaged. In that respect, at many instances the substrates approach and interact with the catalytic centers, but overall, they can be found within the heptane solution in the pore. The IL species (cations and anions) accumulate at the inner surface of the pore, as these interact strongly with the positively charged functional groups on the pore surface. A close inspection of the cation and anion curves shows a larger accumulation of the cations closer to the positively functionalized surface, which at first seems counter intuitive. However,   this process is similar as in the case of electrical double layers (EDL). Even for a positive charge on the surface the bulkier (and fewer) cations accumulate closer to the surface, followed by the anions. The accumulation is especially evident at the bottom half of the pore. Note, that the density results are shown with respect to the axis perpendicular to the pore. In that respect, this perspective accumulates the 3D results on a 2D representation. Accordingly, the IL phase and the heptane do not coexist as falsely assumed by this representation, but still show a phase separation in the pore. Note that, compared to the bulk phase simulation, the pressure was slightly elevated in the pore phase simulation. Therefore, the density of heptane outside the pore shows a slightly higher value. The chloride ions can also be found mainly close to the pore walls in order to play their role as counterions in that region. The catalysts (blue curve) remain inside the pore and the IL phase during the whole simulation, which indicate that these molecules are well dissolved into IL phase inside of the pore and are immobilized by the IL solution.
Compared to the bulk simulations discussed above, in confinement the dynamics of the IL cations, IL anions, and catalysts are significantly slowed down. The MSD for these species in Figure 6(b) clearly show a 10-fold decrease compared to the bulk simulations in Figure 5(b). This clearly underlines the role of confinement within the nanometer-sized pore. From all species in the solution, the catalysts show the slowest diffusion revealing the fact that these are more or less immobilized in the pore. This immobilization is made possible by the presence of the IL solution. (Note again that the poor statistics of the few catalysts in the solution are mirrored on the MSD curve for large times). The chloride ions and the IL anions have a similar diffusion followed by a slightly larger diffusion of the IL cations. On the other hand, the heptane molecules show a diffusion of the same order as in the two-phase solution in Figure 4, which allows for the heptane solution to diffuse into the pore. At the same time, the heptane allows the simultaneous flow of the substrates, having a diffusion of a similar order, into the pore. Accordingly, the substrates can diffuse through the heptane phase and reach the catalyst inside the pore. This is a very important step for chemical reactions, i. e. the catalysis, to take place inside the pore. In short, our simulations have indicated that though initially separated at the pore entrance, the heptane diffuses together with the substrates into the pore, without mixing with the IL. The substrates can access the catalytic center, whose diffusion has been strongly hindered by the confinement and the presence of the IL ions. This is an expected flow of the catalytic reaction within the framework of the SILP technology. Note, that in this work, the focus and aim of our work is to capture the essential conditions prior to the catalytic reaction in order to optimize the environment for the confined heterocatalysis.

Conclusions
In this work, we have performed MD simulations of a two-phase bulk solution, a bulk IL solution, and a confined solution of heptane, IL, catalysts, and substrates within a nanometer sized pore. Our aim was to set up a model for a two-phase SILP system and understand the mechanism of dissolution of catalysts in the IL phase under confinement. The two-phase simulation of heptane and IL did show a stable phase separation of heptane and IL, as well as the dissolution of the catalysts in IL the phase. This is an essential part for the success of the SILP technology, namely the immobilization of catalysts by the IL solution. Furthermore, the catalysts within the bulk IL solution do not show any aggregation, also essential for the reaction to take place homogeneously on within the pore. With these two steps related to the bulk solutions, we have checked whether the conditions for a catalysis in the IL solution can be reached. Accordingly, we moved to the confined system with the pore filled with the solution and the catalysts. We found that the mobility of the catalysts is significantly slowed down due to their electrostatic interaction with the IL species and the functionalized pore wall. The respective results could clearly show the dissolved and immobilized catalysts in the IL phase within the pore. The IL ions are slowed down significantly in the confined space due to the electrostatic interaction between the surface silanol groups and the positively charged functional groups. The heptane phase shows a considerably larger mobility and enters the pore together with the substrate molecules. The latter can thus reach the catalysts in the pore. Accordingly, the reaction between the substrates and the catalysts, i. e. the catalysis, can take place. In the end, the pore is filled with a two-phase solution of IL and heptane, with the catalysts remaining in the IL phase and the substrates mostly diffusing in the heptane phase, but occasionally moving closer to the catalytic centers. The simulations thus provide access to local substrate concentrations that are important for an overall analysis of the catalytic cycle.
In this work, the electrostatic interaction between the surface of the pore and also between the ionic species have played an important role for the dissolution of the catalysts into the IL phase, as well as their mobility. Our aim was not to model the catalytic process itself, rather study the conditions that are important for the catalysis eventually to take place. We have set up a model for this aim, that paves the way for more in-depth investigations and analysis. These will allow to selectively set the conditions in the solution and the pore in order to accelerate the catalytic process. These conditions could refer to the density of the surface groups, the pore size, and type of IL solution, the chemistry of the catalytic centers, etc. The addition of the product molecules i. e. macrocyclic molecules is the next step that will provide insight into local concentrations around the catalytic center under reaction conditions. Furthermore, the refinement of the force-field of the Ru-based catalyst or a further development of atomistic force fields for other catalysts will be challenging tasks that can improve the validity of our simulation results and extend those. In summary, we have made a first step towards setting up the model for a two-phase SILP system and the understanding of the behavior of solution and catalysts in confinement. A number of studies, both theoretical and experimental are further needed to address these issues in order to be able to tune and optimize confined catalytic systems.