Mapping of guest localization in mesoporous silica particles by solid-state NMR and Ab initio modeling: New insights into benzoic acid and p-fluorobenzoic acid embedded in MCM-41 via ball milling

Mapping of Guest Localization in Mesoporous Silica Particles by Solid-State NMR and Ab Initio Modeling: New Insights into Benzoic Acid and p‑Fluorobenzoic Acid Embedded in MCM-41 via Ball Milling Published as part of The Journal of Physical Chemistry virtual special issue “Advanced Characterization by Solid-State NMR and In Situ Technology”. Katarzyna Trzeciak, Sławomir Kazmierski, Kacper Druzbicki,*and Marek J. Potrzebowski* Cite This: J. Phys. Chem. C 2021, 125, 10096−10109 Read Online ACCESS Metrics & More Article Recommendations * sıSupporting Information ABSTRACT: We present a novel strategy for mapping the localization of guest molecules (GMs) in mesoporous materials by combining mechanochemistry with solid-state nuclear magnetic resonance (ssNMR) spectroscopy. To this end, we consider model guest−host systems of benzoic acid (BA) and para-fluorobenzoic acid (4-FBA) embedded in mesoporous MCM-41 material and examine the recently proposed loading (MeLo) procedure for efficient encapsulation of molecular species. Application of high-resolution NMR experiments (1H and 19F NMR) has allowed detection of a multimodal distribution of the spectral signals ascribed to embedded GMs. This peculiarity reflects the presence of distinct molecular ensembles subjected to intrinsically different local environments and exhibiting different dynamical behavior. Furthermore, a considerable fraction of an amorphous phase has been found as a byproduct of the ball-milling. The stability of the phase mixture was further checked by subjecting the samples to chemical and physical stimuli, and a detailed interpretation of the NMR data was corroborated by theoretical calculations. To this end, we have undertaken a challenge to predict the NMR spectra of the GMs@MCM-41 using advanced ab initio molecular dynamics (AIMD) simulations, providing an accurate and exhaustive analysis of the NMR spectra. On the basis of the ab initio modeling validated against the experimental results, we find that the multimodal signal distribution reflects the level of the pore filling, and can be ascribed to the presence of both interface and fluid molecular species trapped inside the pores. This has been confirmed for both BA@MCM-41 and 4-FBA@MCM-41 systems. The presence of the third, amorphous fraction can be linked to interstitial space between randomly ordered crystallites, pointing at the importance of the external surface in further stabilization of encapsulated materials. Altogether, a consistent experimental and theoretical methodology has been presented, paving the way for a more accurate analysis of complex nanoconfined systems. I. INTRODUCTION In recent years, mesoporous silica particles (MSPs) have found a plethora of applications as catalyst carriers and systems for drug transportation. MSPs may be prepared under a wide range of conditions allowing for pore-size tuning (easily adjustable from ca. 20 Å to about 100 Å) and modification of the host composition (e.g., to include other metallic oxides and sulfides as well as aluminophosphates). Such flexibility allows their targeting for many potential catalytic applications, starting from more traditional acid and redox processes, 1 through the gasification of biomass 2 enzyme immobilization, 3 or ending with their use in photocatalysis, 4 to name just a few. Following the seminal paper from Vallet-Regi et al., 5 MSPs have also received a great deal of attention as drug delivery systems (DDSs). 6−13 The attraction and utility of MSPs as DDSs are due to their unique geometrical features, such as high surface area and large pore size, that can be further adapted to meet specific needs. 14,15 Recently, MSPs have been approved as drug carriers by the Food and Drug Administration (FDA), being recognized as safe-to-use oral delivery ingredients in amounts up to 1500 mg per day. 16 Notwithstanding plenty of Received: February 24, 2021 Revised: April 25, 2021 Published: May 3, 2021 Articlepubs.acs.org/JPCC © 2021 The Authors. Published by American Chemical Society 10096 https://doi.org/10.1021/acs.jpcc.1c01675 J. Phys. Chem. C 2021, 125, 10096−10109 Downloaded via CSIC on January 25, 2022 at 11:39:17 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. commercial applications of MSPs, there is often little understanding of the structure and dynamical processes occurring at the level of the silica pores. These phenomena are of great importance to understand drug release kinetics, or a tailored design of novel catalysts, where small molecules can act as valuable probes addressing these questions. This has been illustrated by Dervin et al., 17 studying the structure and diffusion of benzene confined in a mesoporous-silica-based catalyst. The mesoporous material with larger pores compared to microporous zeolites have the potential for greater freedom of movement, where the hydrogenation was found to be internally limited by mass transport. Studying the structure and dynamics of the sorbed molecules in such systems is a key factor facilitating design of better catalysts. Because of the pharmaceutical and industrial relevance and both economic and ecological reasons, the crucial step in the successful application of MSPs is also developing a simple and efficient method for loading requested chemicals into a chosen carrier. 18,19 Among different loading strategies, the method based on the milling of solid components seems to be particularly attractive. Such an approach, called mechanochemistry, or more precisely mechanoloading (with acronym MeLo), is a solid-state method in which neat or liquid-assisted grinding is performed. Recent studies have indicated that the MeLo route can be more effective in preparing solid complexes compared to conventional solution methods. 20 The problem which has to be taken into consideration during the ballmilling processes is the resistance of the nanocarriers to mechanical stress. This issue was well-recognized for MCM-41 (Mobil Composition of Matter No. 41), which is a hexagonally ordered silica with an average pore size of ca. 30 Å. 21 Abu-Zied et al. have studied the effects of ball-milling on the structure, texture, and morphology properties of mesoporous MCM41. 22 The authors have proven that the surface area and the pore volume showed continuous decreases with an increase in milling time. MCM-41 maintains its initial crystallinity until 30 min of ball-milling; thereafter, a severe loss of crystallinity can be observed by grinding the material for a period of 1−2h. Having such experimental data and keeping in mind the unique features of MCM-41, we have chosen this matrix as a proper object for our MeLo studies. Among a plethora of spectroscopic techniques used for the characterization of nanoconfined species, NMR spectroscopy emerges as the leading method, providing means to elucidate the local structure and dynamics of model GMs at different spatial and temporal scales. Recent examples include studies of simple liquids (with the prime example of water, 23,24 small nonpolar species, 25−30 or polar aromatic guests, 31−37 to name just a few (for an extensive review, see ref 38). Furthermore, as illustrated by Pruski and co-workers, state-of-the-art solid-state NMR (ssNMR) experiments, including fast magic angle spinning (MAS), 39−42 dynamical nuclear polarization (DNP), and surface-enhanced NMR spectroscopy (SENS), 43−45 allow elucidating the structural information on the functionalization of the MSPs surface with an unprecedented level of detail. In common opinion, supported by numerous studies cited above, MCM-41 is considered as a very well-organized homogeneous functional material. The pores are very uniform causing narrow pore size distributions. Moreover, the pores are unidirectional and arranged in a honeycomb structure over micrometer length scales. This simple picture was recently discussed by Limbach and co-workers, 46 proving that the nature of MCM-41 is more complex. This complexity stems from the presence of at least two types of surfaces illustrated in Figure 1. The internal pore surface is very large compared to the external surface of the crystallites, which has usually been neglected. To date, only a little is known about the external volumes, which might depend on the way how the crystallites are compacted and subjected to different forces even at the stage of physicochemical sample characterization. For instance, MAS NMR spectroscopy is a primary analytical technique in this regard. However, the crystallites can be displaced toward the rotor walls by centrifugal forces affecting the interstitial space in a powder sample (see Figure 1). Our project aims to answer the question of whether such phenomena and different localization of the guest in the host matrix can occur during the mechanoloading of solid organic molecules into the pores of MCM-41. In the present work, we focus on the use of high-resolution 1H NMR spectroscopy to study benzoic acid (hereafter, BA; see Figure 2a). BA along with ibuprofen (hereafter, Ibu) are considered as model GMs for studies of encapsulation of hydrogen-bonded active pharmaceutical ingredients (API) in DDSs. In this regard, we highlight the importance of a particularly insightful series of papers by Azais and co-workers, Figure 1. Illustration of internal and outer-space accessible to GMs in the powder samples of MCM-41 silica according to TorresBarthelemy et al. 46 A single pore, with a diameter of ca. 3 nm and a depth much greater than ≫100 nm, give rise to a large internal volume (left). These pores are ordered in a regular honeycomb manner (middle), resulting in a non-negligible external surface (the gray area marked with the vertical arrow). The external surface give rise to interstitial space (right) formed between the randomly ordered crystallites (see the horizontal arrows). Figure 2. Molecular geometry of (a) BA and (b) 4-FBA selected as model GMs in the present study. The Journal of Physical Chemistry C pubs.acs.org/JPCC Article https://doi.org/10.1021/acs.jpcc.1c01675 J. Phys. Chem. C 2021, 125, 10096−10109 10097 providing a comprehensive NMR analysis of model GMs embedded in silica matrices with different sizes of pores. 47−50 To solidify our findings, we extend the analysis toward a fluorinated BA analog, i.e. 4-fluorobenzoic acid (hereafter, 4FBA; see Figure 2b). Much needed insights into the local structure and dynamical events occurring in the guest−host systems require the use of accurate modeling framework, 51−57 which is, however, limited by the construction of a realistic model accounting for their complexity. In this regard, only a few papers appeared to date, providing alternative realistic models of amorphous silica and MCM-41. Most importantly, Ugliengo et al. have introduced a realistic, large-scale and fully periodic model of MCM-41 (based on the solid-state B3LYP calculations), which successfully accounts for an average stoichiometry, infrared spectroscopy, and diffraction data. 58 Alternatively, the authors have also proposed a well-established, simplified slab model of the amorphous silica surface, validated against several physicochemical properties. 59 In our work, we utilize both structural models proposed by Ugliengo and co-workers to shed more light on the loading ability of the MCM-41 matrix, the local structure around GMS and the associated spectral response. For the very first time, advanced ab initio modeling is employed to interpret the high-resolution MAS NMR spectra recorded for the GMs embedded in MCM-41 (hereafter GMs@MCM-41; GMs = BA and 4-FBA). Our paper is constituted as follows: (i) First, we present an exhaustive ssNMR study of BA embedded in MCM-41 using the MeLo procedure, and we solidify our observations by extending the research toward a fluorinated BA analog (4FBA). (ii) In a further step, we present the interpretation of the spectral peculiarities based on advanced time-dependent ab initio modeling focused on BA@MCM-41. (iii) Finally, we examine the impact of the chemical and physical stimuli on the behavior of the encapsulated molecules. II. METHODS II.1. Sample Preparation. MCM-41 was purchased from Sigma-Aldrich and calcined at 300 °C for 1 h to remove water. BA and 4-FBA were used as received from Tokyo Chemical Industry (TCI). The physical mixtures of BA and calcined silica MCM-41 were mechanically ground (30 min at 450 rpm) in a planetary ball mill (Retsch PM 200) using the stainlesssteel jar (12 mL) and 10 stainless-steel balls with o.d. = 4 mm). The MeLo procedure was performed at various weight ratios (from 1:10 up to 1:1). A similar experimental procedure was applied to prepare the 4-FBA@MCM-41 samples. The physicochemical characterization of BA@MCM-41 mixtures (thermogravimetric analysis (TGA), differential scanning calorimetry (DSC), and gas-sorption analysis) is further presented as Supporting Information; see Sections S1.1. and S1.2. II.2. NMR Spectroscopy. The MAS NMR experiments were performed on a Bruker Avance III 600 spectrometer, operating at 600.13 MHz for 1H, and 150.90 MHz for 13C and employing a 1.3 mm MAS dual-channel 1H/BB (15N−31P) probe head working with 1.3 mm ZrO2rotors. The measurements were conducted under 60 kHz MAS conditions. The collected spectra were externally referenced to TSP (trimethylsilylpropanoic acid) (0 ppm) and the carbonyl atom of glycine (δ(C=O) = 176.50 ppm) for 1H and 13C, respectively. 19F MAS NMR spectra were recorded on a 600-MHz Bruker Avance III spectrometer, equipped with a 1H, 19F, and 13C 2.5 mm triple-resonance MAS probe head operating at 564.69 MHz. The 90°pulse duration was set to 2.25 μs. For the experiment with 19F decoupling, a π−pulse scheme was used with a 180°pulse equal to 6 μs. The collected spectra were externally referenced to CFCl3(0 ppm). For the 1H MAS spectra, 32k data points FIDs were collected with a spectral width of 36.5 kHz (61 ppm) and 90° pulse of 1.5 μs. Sixteen scans were accumulated per experiment with a relaxation delay of 3 s. 2D RFDR MAS experiments were performed with mixing time pulses synchronized with rotor cycles. Here, the spectra were accumulated with a spectral width of 13.2 kHz (22 ppm) in both dimensions and 90°pulse of 2.35 μs. Usually, 128 or 256 FIDs, 14k data points, and 16 scans were accumulated per experiment, with repetition delays of 3 and 60 s. Mixing times were set between 1.6 and 512 ms. For temperature stabilization, the temperature of the bearing gas was stabilized and controlled with the BCU II cooling unit, controlled with the BVT3200 variable temperature unit. The spectral data were acquired and processed using the TopSpin 3.1 program. High-Resolution MAS NMR experiments were performed on a BRUKER Avance III 400 Wide Bore spectrometer, operating at 400.13 for 1H. The 4 mm triple resonance 1H/13C/31P HR MAS probe head with deuterium (2H) lock channel and magic angle gradient coil were used. The 2D data were collected in 4k ×256 (F1 ×F2) data points matrix with eight scans per experiment, and a relaxation delay of 2 s. Spectral width was set to 5600 Hz (14 ppm) in both dimensions with 1H90°. The pulse of 7.2 μs before Fourier transformation data was apodized with a square sine (QSINE) bell shifted function in both dimensions, and zero-filled up to 4 ×1k data points matrix. Mixing times were set between 50 and 350 ms. For temperature regulation, the bearing gas was cooled by a liquid nitrogen heat exchanger. The temperature was controlled and stabilized with the VTU 3000 temperature unit. Before starting the measurements, the sample was kept at the desired temperature for at least 10 min. The spectral data were acquired and processed using the TopSpin 2.6 program. II.3. Theoretical Calculations. Theoretical calculations were performed in Periodic Boundary Conditions (PBC) using a“hybrid computational scheme”, where ab initio structural modeling (geometry optimization and molecular dynamics (MD) simulations) was performed using CP2K, 60,61 and on top of these calculations, CASTEP was used to model the NMR response. 62 Three cases were considered: (i) static calculations with MCM-41 filled with various content of BA (relying on an extended structural model of MCM-41 by Ugliengo et al. 58 ), used to account for sorption properties of the studied host matrix (presented as the Supporting Information; see Section S2.1); (ii) the reference NMR calculations, where the BA crystal (CCDC identifier: BENZAC02) 63 along with discrete molecular models (centered in a large, cubic box with a=35 Å) were used to test the numerical accuracy of the “hybrid computational scheme”(presented as Supporting Information); and (iii) the ab initio MD (AIMD) simulations, relying on a simplified two-dimensional slab-model of amorphous silica according to Corno et al., 59 used for a detailed interpretation of the ssNMR data presented in the main text. For the ab initio modeling the computational details are given as follows. In CP2K, exchange and correlation were approximated with a generalized-gradient-approximation-type (GGA) functional from Perdew−Burke−Ernzerhof (PBE), 64 The Journal of Physical Chemistry C pubs.acs.org/JPCC Article https://doi.org/10.1021/acs.jpcc.1c01675 J. Phys. Chem. C 2021, 125, 10096−10109 10098 augmented with Grimme’s D3 dispersion corrections. 65,66 A hybrid basis set combining plane-waves (PWs, with 400 Ry cutoff) and Gaussian atom-centered functions (DZVPMOLOPT) along with standard Goedecker−Teter−Hutter (GTH) pseudopotentials (PPs). The self-consistent field (SCF) was converged to an accuracy better than 1 ×10−12 Ha/atom, and the geometry optimization was performed to converge the residual Hellmann−Feynman forces down to 1 × 10−7Ha/bohr. The CP2K modeling also involved extensive AIMD simulations, keeping the same numerical conditions except for the SCF convergence, which was reduced down to 1 ×10−7Ha/atom. An integration step was set to 0.5 fs. A 10 ps equilibration was performed before collecting the 50 ps long production runs. The CASTEP calculations on top of the CP2K outputs imposed the PWs/PPs formulation of DFT, using the PBE-TS functional. 64,67 The output CP2K trajectory was sampled every 0.5 ps, and the NMR parameters were calculated from 100 snapshots. The calculations of the NMR shielding tensors were performed with the gauge including projector-augmented wave (GIPAW) method. 68 The allelectron information was reconstructed with scalar relativistic UltraSoft PPs (USPPs) generated on the fly. 69 The electronic wave functions were defined using a PW basis set with a 600 eV kinetic energy cutoff. The k-point sampling was limited to the Γ-point, unless otherwise stated. The SCF was converged to 5 ×10−10 eV/atom. The GIPAW outputs were processed with the help of the MagresView code by Sturnollio et al. 70 The calculated isotropic components of the magnetic shielding tensors for 13C, and 1H nuclei were referenced to standard values of 170 and 30 ppm, respectively, to obtain the chemical shifts. For the 19F NMR calculations, the results were referenced to 128 ppm. III. RESULTS AND DISCUSSION III.1. Mechanoloading (MeLo) of the GMs into MCM41. Application of the MeLo route allows monitoring evolution of the 1H NMR response in function of the guest−host weight-to-weight proportions. Figure 3 depicts the set of 1H MAS NMR (60 kHz) spectra recorded for the physical mixtures of BA with MCM-41 ground at different weight ratios (from 1:10 to 1:1). By inspection of the figure, one can note striking dual contributions to the 1H NMR spectra in the range 8−6 ppm, which is due to phenyl ring (Ph) protons marked as Xand Yfor the two different fractions detected at this point, respectively (with an obvious assignment to ortho-,para-, and metaspecies). At this point, it has to be stressed out that such a puzzling effect was also observed by Azais et al. for the BA encapsulated in MCM-100. 49 The authors assumed a bimodal distribution of the large pores trying to explain this finding, which, however, can be ruled out in our case (see Section S1 of the Supporting Information). Intriguingly, an inversion of the intensity relations for the fractions Xand Yis observed here at the ratio of 1:3, which evolves into a unimodal distribution for sample 1:2. Prior to this transition, one observes a continuous evolution of the spectral features. The 1H NMR spectra for the mixtures 1:10−1:4 are of remarkable similarity but manifest a progressive change in the spectral range below 6.5 ppm upon the increase of the BA content. As a possible explanation, we assume that the progressive downfield shift is related to the interaction of the GMs with the Si−OH groups via hydrogen bonding. The phenyl ring proton signals evolve toward lower ppm, whereas the opposite trend is observed for the broad band centered at ca. 5 ppm. This synergistic behavior can be easily understood by assuming that the Xmolecules are directly linked to the pore surface. Anticipating that sample 1:10 barely affects the native hydrogen-bond structure of MCM-41, these data can be treated as the reference. However, once the filling increases, more Si−OH groups contribute to the hydrogen bonding with the GMs, eventually reaching the equilibrium at a weight ratio of 1:5. We note an abrupt vanishing of the silanol contributions for sample 1:3, where the signals are barely detectable. Once the Si−OH groups are hardly observed, the contribution from the COOH protons above 8 ppm becomes visible. This can be linked to a quantitative effect, where the prevalent amount of the BA protons hampers observation of the Si−OH signals from the matrix. To gain further quantitative insights, we have performed nonlinear curve fitting to the NMR data. The representative results are displayed in Figure 4. In each case, the spectrum was described as the sum of Lorentzian functions. Since the resulting spectra are very complex, this process became very challenging, calling for the use of some constraints to reduce the number of the fitted parameters. Therefore, for each fraction considered, the full width at half-maximum (FWHM) was kept the same for phenyl protons (i.e., FWHMo= FWHMm= FWHMp). Furthermore, these FWHM’s were fixed to the value derived from the spectrum of the mixture for which a given fraction dominates (i.e., 1:10 for fraction X; and 1:2 for fraction Y). Also, in each case, the relation between the Lorentzian integrals (A) was defined as Ao=A m=2Ap= 2ACOOH. A detailed inspection of the data highlights the presence of Xand Yfractions. However, a careful decomposition of the spectra has proven the presence of a third fraction (hereafter, Z), which manifests with visibly broader signals (with nearly four-times larger FHWM that for Xand Y), which cannot be ignored. To understand the origin of each fraction, we refer to theoretical analysis of the BA-uptake capacity of MCM-41 presented in Section S2 of the Supporting Information. Figure 3. High-resolution 1H MAS (60 kHz) NMR spectra of physical mixtures of BA and MCM-41 (30 Å), varying in the considered weight ratios (from 1:10 to 1:1). The Journal of Physical Chemistry C pubs.acs.org/JPCC Article https://doi.org/10.1021/acs.jpcc.1c01675 J. Phys. Chem. C 2021, 125, 10096−10109 10099 According to this analysis, low content of GMs in the mixture 1:10 corresponds to BA molecules defined as separated units at the pore surface, and so the assignment of the Xfraction (nearly solely observed in the corresponding 1H NMR spectrum) to the low-density surface fraction becomes straightforward. Consequently, we can assume that the BA content in the mixtures 1:10−1:4 is up to the level of monolayer formation (see Figure S8 in Section S2.1. and related discussion), while for higher BA concentration, a denser fraction is formed (Y). To understand a peculiar picture of sample 1:1 (see Figure 4a), we note that the amount of BA used here largely exceeds the sorption capacity of the internal volume of MCM-41 (see the discussion in Section S2). Consequently, there must be a fraction of molecules adsorbed at the outer surface. Indeed, a detailed inspection reveals a visible broadening of the 1H NMR spectrum for sample 1:1, which suggests the presence of an unloaded, bulk fraction (Z). Further spectral analysis indicates the presence of both crystalline and amorphous contributions (as deduced from the position of the carboxyl signals in Figure 4a). The presence of the glass fraction is not surprising, since ball-milling is considered as well-established method for sample amorphization. At this point, we recall the recent findings from Limbach and co-workers presented in ref 46. As mentioned above, the external surface gives rise to the interstitial space between the randomly ordered crystallites of MCM-41. This becomes intriguing when confronted with the well-known instability of amorphous BA, clearly suggesting that the interstitial surface has the potential ability to stabilize amorphous BA. As found by Smirnova et al. in the studies of crystallization of BA inside the pores of aerogels, 71 strong interactions with the surface favor the amorphous form, whereas weak interactions favor crystalline particles. We can anticipate that a rough outer surface of the MCM-41 particles, subjected to mechanical grinding, preserves its polar character. Interestingly, the crystalline fraction was only detected for sample 1:1, as confirmed by both NMR and DSC analysis (see the Supporting Information), which might suggest that the amorphous fraction in the interstitial space eventually crystallizes when the BA is present in a considerable excess, following the mechanism proposed by Smirnova et al. 71 The Zfraction is observed for all other BA@MCM-41 samples analyzed here, without any trails of crystalline phase, and here onward the Zfraction will be ascribed to the amorphous phase. We also note a similarity in the peak positions of the amorphous fraction Zand the fraction Y. Since both phases are similar in terms of the signal distribution but largely differing in mobility, it is logical to assume that the Figure 4. Results of the nonlinear curve fitting to the representative 1H MAS (60 kHz) NMR spectra for the physical mixtures of BA and MCM-41 (30 Å), varying in the considered weight ratios, 1:1 (a), 1:2 (b), 1:3 (c), and 1:4 (d). Mind that the sample 1:1 includes a considerable fraction of unloaded species; and note an abrupt inversion of the signal intensities for the samples 1:3 and 1:4. The fractions Xand Yare presented as black and red areas, being assigned to the surface (interface) and fluid species, respectively. The blue shaded curves correspond to the amorphous phase (Z). See main text for further details. The “excess”curve denotes an extra function added to improve the fit and account for the high asymmetry of the COOH band. The Journal of Physical Chemistry C pubs.acs.org/JPCC Article https://doi.org/10.1021/acs.jpcc.1c01675 J. Phys. Chem. C 2021, 125, 10096−10109 10100 fraction Yis characterized by a similar local environment, yet, representing a more mobile ensemble. A downward shift on the ppm scale suggests that the fraction Yis slightly denser than the amorphous one, and so, most likely, more spatially confined. We, therefore, call the fraction Yafluid phase confined inside the MCM-41 pores. According to Figure 4 b, sample 1:2 emerges a nearly sole fraction Y. Hence, the mixture is considered as the one of possibly maximal loading (around 40−50 molecules per unit volume according to an extended Ugliengo model), with the BA condensed far beyond the monolayer level (see Figure S9). As seen in parts c and d of Figures 4, the experimental spectra of the samples with a lower BA content can be successfully described by accounting for the presence of all three phases. Nevertheless, we note that the quality of the fit becomes slightly imperfect for samples 1:4−1:10, which, according to Ugliengo’s model, would correspond to the half-filling of a pore (up to the formation of the BA monolayer). This manifests by a slight broadening of the Ysignals, which was not accounted for in the fitting procedure as leading to numerical instabilities. To solidify our findings, we have repeated the experiments for an analogous GM, 4-FBA@MCM-41. Selected weight-toweight proportions were used as follows: 1:1, 1:2, 1:3, 1:5, and 1:10. The results are presented in Figure 5. According to Ugliengo’s model, the number of embedded molecules per unit volume would be estimated as 68, 34, 23, 14, and 7, respectively, i.e., slightly lower than for the BA@MCM-41 (see the Supporting Information). Figure 5 a presents the 1H MAS (60 kHz) NMR spectra of the 4-FBA@MCM-41 mixtures, whereas Figure 5b augments the analysis with 19F NMR spectroscopy. Both results confirm that our previous findings on BA@MCM-41 do not constitute an isolated case, providing the same trends. The upward shift of the 1H resonances on the ppm scale may indicate that the observed fractions are less dense and the molecules are more spatially separated. This can be supported by the results of the studies of solid solutions of BA and 4-FBA, 72 showing the visible expansion of the cell volume with growing content of 4-FBA (up to 6%). For 4-FBA@MCM-41, an inversion of the signal distribution is observed for the mixture 1:5, and the nearly full loading is observed at the critical weight ratio of 1:2. As shown in parts c and d of Figure 5, the spectra were successfully fitted with the contribution from all three fractions. We also note that one could not detect the COOH signal for the fraction X, Figure 5. (Upper panels) High-resolution MAS (60 kHz) (a) 1H NMR and (b) 19F NMR spectra of physical mixtures of 4-FBA and MCM-41, varying in the considered weight ratios (1:1 to 1:10). The mixtures were mechanically ground in the planetary ball mill Retsch PM 200. (Bottom panels) Results of the nonlinear curve fitting for the representative 1H MAS NMR spectra of the mixtures 1:3 (c) and 1:5 (d). The fractions Xand Yare presented as black and red areas, being assigned to the surface (interface) and fluid species, respectively. The blue shaded curves correspond to the amorphous phase (Z). See main text for further details. The Journal of Physical Chemistry C pubs.acs.org/JPCC Article https://doi.org/10.1021/acs.jpcc.1c01675 J. Phys. Chem. C 2021, 125, 10096−10109 10101 which may stem from the water-promoted fast proton exchange with the hydroxylated silica surface. 50 Finally, we provide a quantitative analysis of percentage contributions from each fraction for both GMs@MCM-41 studied here (see Figure 6). While analyzing these contributions quantitatively, it should be noted with caution that each fraction manifests its presence with a different half-width of the respective signal, which depends on the molecular mobility due to different degrees of averaging the dipolar interactions. Hence, one needs to correct for different FWHM’s by reducing the values to a common denominator (see Figures 4 and 5, for the FWHM ratios). From inspection of Figure 6, it is found that the outer-space fraction Zcontributes with the content below 20% for all the samples considered but the 1:1 mixture, where the amount of the GMs largely exceeds the sorption capacity of MCM-41 (30 Å). We, however, note that in the special case of sample 1:1, relative contributions might be, somehow, questionable due to overwhelming broadening of the signal. For a typical MCM-41 material, the mesopores surface area (Smeso) is estimated as 1049 m2/g, whereas the external surface area (Sexternal) reaches 89 m2/g. 73 We can, hence, estimate the maximal in-pore adsorption as Smeso/Stotal = 92%. According to Figure 6, the ball-milling with the BA@ MCM-41 weight ratio of 1:2 allows one to reach this limit, with the Ycontribution found at the level of ca. 90%. In this case, the contributions from Xwere not detected. We can, hence, treat the 1:2 sample as a uniform fluid phase Y(in-pore confinement far beyond the formation of monolayer). An abrupt formation of the pure Yfraction indicates that the mass ratio 1:2 corresponds to some critical concentration of GMs, for which the affinity to condense far beyond the monolayer level is energetically favorable. For the other samples, the content of the inner-surface fraction Xdecreases linearly prior to the full pore filling. Finally, there is a remaining question if the fractions Xand Y coexist within the same pore in the mixtures with lower content of GMs (e.g., 1:3 and 1:5). Gubbins and co-workers earlier presented grand canonical Monte Carlo (GC-MC) simulations of benzene adsorption in a cylindrical silica nanopore (36 Å in diameter). 74 The results indicate that for a partially hydrogenated system the aggregates formed beyond the monolayer level appear before fully filling the core space. It is, however, questionable in the case of BA and 4-FBA, according to a larger size of the GMs. In such alignment, the distance between Xand Ymolecules in MCM-41 pore with diameter ca. 30 Å should be short, i.e., below 5 Å. Such distance can be detectable by measurement of 1H−1H dipolar interactions using two-dimensional ssNMR techniques. The 2D 1H−1H radio frequency driven recoupling (RFDR) MAS NMR spectra recorded for BA@MCM-41 and 4-FBA@MCM41 (1:3 weight ratio) are projected in parts a and b of Figure 7, respectively. We have chosen this technique because during the mixing time the magnetization transfer is much more efficient in the RFDR sequence as compared to nuclear Overhauser effect spectroscopy (NOESY). We have carried a series of measurements with a broad range of mixing times, varying from 1.6 to 512 ms. In no single case the contacts between molecules Xand Ywere detected, which suggests that both ensembles are clearly spatially separated. A similar conclusion is also valid for NOESY MAS experiments, performed with the spinning rate of 60 kHz. Figure 6. Quantitative analysis of the percentage contributions from each GMs fraction present in the mixtures of BA (a) and 4-FBA (b) with MCM-41(30 Å), varying in the considered weight ratios (1:1 to 1:10). The percentage contributions from the fractions Xand Yare presented as black and red bars, being assigned to the surface (interface) and fluid species, respectively. The blue bars denote percent contributions from the amorphous phase (mainly crystalline for the special case of sample 1:1, see main text for details). Figure 7. 1H−1H RFDR MAS spectra of the 1:3 BA@MCM-41 (panel a) and 1:3 4-FBA@MCM-41 (panel b) samples, recorded with a mixing time of 64 ms. The MAS frequency was 60 kHz for the BA@ MCM-41 sample and 24 kHz for the 4F-BA@MCM-4 one. 19F−19F RFDR MAS spectrum of the 1:3 4-FBA@MCM-41 sample recorded at 24 kHz MAS, with a mixing time of 50 ms (panel c). The Journal of Physical Chemistry C pubs.acs.org/JPCC Article https://doi.org/10.1021/acs.jpcc.1c01675 J. Phys. Chem. C 2021, 125, 10096−10109 10102 On the other hand, one has to be reminded that the efficiency of dipolar interactions strongly depends on the molecular dynamics. Fast in the NMR time-scale, largeamplitude motions significantly reduce the dipolar interactions. Usually, such motions are thermally dependent, where higher temperatures accelerate the reorientation of molecules. The very fast sample spinning generates an increase in temperature in the measured sample due to the friction effect. This phenomenon can be responsible for the lack of correlation peaks confirming Xand Yinteractions. Nevertheless, one could expect that for the 4-FBA@MCM-41 sample, the 19F−19F contacts might be less sensitive to some types of motion due to local molecular symmetry (e.g., rotation about the long molecular axis), while the mass effect due to the presence of fluorine atoms might partially slow down other types of motions (e.g., rotation about the short molecular axis). To this end, we have recorded a two-dimensional 19F−19F NMR spectrum (see Figure 7c), where, however, also no off-diagonal correlation peaks can be found. It is, hence, reasonable to assume that the abrupt inversion of the intensity distribution for the mixtures 1:2 can be linked to the transition from monolayer into full loading, and so the fractions Xand Y observed here simultaneously are present in different pores. An alternative explanation would assume that the lack of correlation in the 2D spectra may be due to different localization of the GMs in the MCM-41 matrix as shown in Figure 1. For instance, with the Xfraction being adsorbed solely at the external surface, where the Ysignals corresponding solely to the inside pore space. However, such an assumption would not clarify the presence of the amorphous fraction Z, and we do not find further arguments to support such a hypothesis. III.2. Insights from Ab Initio Modeling. Further shreds of evidence on the origin of the bimodal Xand Ycontributions have been derived from theoretical calculations. The largest model from Ugliengo used to analyze the sorption ability consists of 1344 atoms (see the Supporting Information), putting the NMR calculations to the technical limit. In this case, calculations of NMR shielding tensors required 2 days with ca. 500 processors at the high-performance computing facility PROMETHEUS. 74 This allowed us to further use this model for more advanced calculations. Instead of using the cylindrical pore model from Ugliengo et al., we, hence, turn to a simplified slab structure proposed by Corno et al. 59 The adopted model is defined by the wall thickness of 7.2 Å, the unit cell composition of H22O63Si26,andtheunitcell parameters: a= 12.60 Å, b= 12.83 Å, α=β= 90.0°, and γ = 83.1°. The model has been fully relaxed at the atmospheric pressure and the separation of silica layers was fixed to ca. 30 Å (c-axis). The silanol density is defined as 4.5 ×OH per nm2, which corresponds to eight SiOH surface groups per unit cell. To account for motional averaging, we have employed ab initio MD simulations (AIMD) at room-temperature conditions. The simplified slab model was filled by low (8 molecules) and high (22 molecules) content of BA. These structures are depicted in Figure 8. In the case of 4-FBA, the simplified AIMD simulations were performed to solely focus on the 19F signals. The starting model was build based on the last snapshot from the BA@MCM-41 trajectory and by replacing the proton in the para-position with the fluorine atoms. The short-time scale of the AIMD simulations was able to cover the vibrational dynamics along with fast molecular librations. The low-content model of BA@MCM-41 (Figure 8 a) represents solely the Xfraction, whereas the high content model describes both Xand Yensembles (note that the interface molecules were excluded from subsequent analysis of the Yfraction). Although the adopted models do not describe the cylindrical shape of the pore, however, they do account for the mutual and surface interactions of the BA molecules. From inspection of the right panels in Figure 8, it becomes evident that the main difference between Xand Yfractions stems from the local environment. The open space is interfaced with the X ensemble, whereas the Yfraction is fully embedded in the field of the surrounding molecules. This is of key importance for understanding the difference in the NMR response of both phases observed experimentally. The final evidence is provided by the analysis of the distribution of chemical shifts calculated through AIMD (see Figure 9). The chemical shifts were calculated every 0.5 ps, giving 100 spectra for each model considered. Eight (see Figure 8a) and 15 (Figure 8b) BA/4-FBA molecules were used to provide the statistic means for the Xand Yensembles, respectively (at least 800 counts for each nucleus were considered). The analysis of the chemical shift distributions presented in Figure 9 directly supports our assignment of the X and Yfractions. The relative differences in the predicted 1H NMR signals for the protons in the ortho-position (BA; at the level of ca. 0.3 ppm) and the 19F chemical shifts (4-FBA; at the level of 3 ppm) are corroborated by NMR experiments. The calculations are, however, unable to discern the hydrogens in the meta-andparapositions (BA), providing virtually identical chemical shifts (within 0.1 ppm). This might be ascribed to somehow limited accuracy of the adopted computational scheme (e.g., double-ζquality basis set used Figure 8. Models of 8 (a) and 22 (b) BA molecules embedded in MCM-41, as used in the AIMD simulations. These were based on a simplified model of amorphous silica according to Corno et al. 59 In each case, the unit cell (see the dotted squares in left panels) has been doubled for an ease of visualization. In the left panels, the first snapshot from the 50 ps production run is displayed in black and superimposed with the snapshots of the whole run, sampled every 2.5 ps. In the right panels this has been augmented by superposition of the intermolecular reduced density gradient (RDG) isosurfaces, 76 calculated for each of the time-evolved structure of BA. The RDG maps were calculated utilizing pro-molecular electron density and scaled by 0.4 for an ease of visualization. The significantly attracting forces (hydrogen bonds) are colored in blue, whereas the repulsive interactions are shown in red. Weak van der Waals forces are marked in green and yellow, following their repulsive nature. The surface/ interface and the fluid fractions are denoted as Xand Y, respectively, and marked with horizontal arrows. The Journal of Physical Chemistry C pubs.acs.org/JPCC Article https://doi.org/10.1021/acs.jpcc.1c01675 J. Phys. Chem. C 2021, 125, 10096−10109 10103 in the CP2K simulations). We are, hence, unable to tackle such subtle chemical shift differences, which would, perhaps require much larger statistical averages and a significantly larger basisset used in electronic structure calculations. III.3. Influence of Chemical and Physical Stimuli on the Behavior of GMs Trapped Inside the MCM-41 Pores. At the final stage, the physical properties of the BA@MCM-41 mixture are further studied by analyzing the influence of internal and external stimuli (variable temperature and diffusing solvents) on the behavior of BA embedded in the pores. First, we examine the above-discussed influence of temperature on the sample 1:3 BA@MCM-41. We have carried out the NMR measurements in a broad range of temperatures (253−353 K). Figure 10 (left column) shows the 1H NMR MAS spectra recorded at temperatures 283, 273, 263, and 253 K, respectively. From inspection of the spectra, it is apparent that with decreasing temperature the width of the resonance lines increases. At 253 K, the signals for Xand Ycomponents are very broad and partially overlapped, which supports vitrification at low-temperatures. 48,50 This effect is due to restriction in molecular mobility of both components, their “frozen”localization and, in consequence, dispersion of chemical shifts. The freezing of molecular motion prompted us to test the intermolecular contacts via analysis of the dipolar interactions. Figure 10 (right column) displays the NOESY spectrum of the 1:3 complex (BA@MCM-41) recorded at 263 K. As in the previous case (Figure 7), we did not observe any contacts between the Xand Ymolecules, which supports the reasoning that both fractions are spatially separated. Going further with inspection of the external stimuli on the behavior of BA embedded in the MCM-41 pores, we examine the influence of diffusing solvents. To this end, we select the BA@MCM-41 sample in a 1:5 weight ratio, with similar contribution from both Xand Ycomponents. Figure 11 Figure 9. Histograms presenting the distributions of the 1H (a) and 19F (b) chemical shifts for the BA and 4-FBA molecules embedded in MCM41, respectively. The chemical shifts were calculated for snapshots extracted every 500 fs from the 50 ps long AIMD trajectories. The data are binned with 50 points (shaded bars), and the normal distributions are presented as shaded areas (a) and solid lines (b), respectively, with the labels indicating the mean values obtained from the statistical average. Figure 10. (Left) 1H MAS NMR spectra of BA@MCM-41 sample with 1:5 weight ratio recorded at low temperatures as indicated in the figures. (Right) 1H−1H NOESY MAS spectrum recorded at temperature 263 K with mixing time 250 ms. Figure 11. 1H (60 kHz) MAS NMR spectra of the ground physical mixture of BA and MCM-41 (1:5 weight ratio) processed differently. According to the labels, the pristine sample (BA@MCM-41 1:5) was treated either with the vapors of (a) deuterated (EtOH-d6), (b) hydrogenous ethanol (EtOH), or (c) deuterated water (D2O). The resulting samples were thermally treated in the oven at 80 °C (a) or dried over P2O5(c). Each spectrum has been labeled according to the dominant component (X,Y, and S(EtOH) or S′(D2O) for the solvated species). The Journal of Physical Chemistry C pubs.acs.org/JPCC Article https://doi.org/10.1021/acs.jpcc.1c01675 J. Phys. Chem. C 2021, 125, 10096−10109 10104