Vol.:(0123456789)1 3 Adsorption (2022) 28:125–136 https://doi.org/10.1007/s10450-022-00356-w An atomistic view on the uptake of aromatic compounds by cyclodextrin immobilized on mesoporous silica Hamzeh Kraus1 · Niels Hansen1 Received: 1 October 2021 / Revised: 18 January 2022 / Accepted: 24 January 2022 / Published online: 7 March 2022 © The Author(s) 2022 Abstract The effect of immobilized β-cyclodextrin (bCD) molecules inside a mesoporous silica support on the uptake of benzene and p-nitrophenol from aqueous solution was investigated using all-atom molecular dynamics (MD) simulations. The calculated adsorption isotherms are discussed with respect to the free energies of binding for a 1:1 complex of bCD and the aromatic guest molecule. The adsorption capacity of the bCD-containing material significantly exceeds the amount corresponding to a 1:1 binding scenario, in agreement with experimental observations. Beside the formation of 1:2 and, to a lesser extent, 1:3 host:guest complexes, also host–host interactions on the surface as well as more unspecific host–guest interactions govern the adsorption process. The demonstrated feasibility of classical all-atom MD simulations to calculate liquid phase adsorption isotherms paves the way to a molecular interpretation of experimental data that are too complex to be described by empirical models. Keywords Mesoporous silica · Cyclodextrin · Adsorption · Molecular dynamics simulation 1 Introduction Cyclodextrins (CDs) are cyclic oligosaccharides with coni- cal shape built of (1-4)-linked α-d-glucopyranoside units. Their hydrophilic outer surface renders them soluble in water, while the hydrophobic interior can accommodate hydrophobic compounds. Therefore, they are employed in the removal of pollutants from aqueous media [1, 2], as chi- ral selectors [3] or in (asymmetric) supramolecular catalysis [4–10] among many other fields of applications including the use of cyclodextrins as drug delivery agents or as build- ing blocks in the design of dynamic and adaptive materials [11, 12]. By tuning the size and shape of the hydrophobic cavity through derivatization of the native cyclodextrins the selectivity towards the target compounds can be increased [13]. In addition, the ease of functionalization of their scaf- fold allows to introduce additional catalytic groups or bind- ing sites in specific positions [14]. Together with the solvent used and a possible immobilization of the cyclodextrins on a solid support, a complex multidimensional design problem emerges. Several approaches can be envisioned to produce CD-functionalized materials. The first involves cross-linking of CDs into polymers using C-OH linkers [15]. A second approach utilizes the coating or grafting of CD moieties onto a stationary phase such as organic polymers or silica gel [16]. A third type of material is based on a mesoporous silica support. Mesoporous silicas with chemically attached macrocyclic moieties were successfully prepared by sol-gel condensation of tetraethyl orthosilicate and β-cyclodextrin- silane in the presence of a structure directing agent [17–20], resulting in silica-based materials that possess a uniform framework mesoporosity with defined nanoscaled cavities. The ability of removing aromatic compounds from an aque- ous phase was investigated, and it was concluded that the synthesized materials is promising for this purpose [17, 18, 20]. Molecular modelling approaches are particularly feasible to study the delicate balance between the various intermolecular forces determining macroscopic behavior and allow a funda- mental understanding at the molecular level. Gas adsorption in microporous materials such as zeolites [21, 22] or other nano- structured solids [23] has been studied for more than 45 years in particular using classical molecular simulations. The com- bination of molecular simulation techniques with experimental * Niels Hansen hansen@itt.uni-stuttgart.de 1 Institute of Thermodynamics and Thermal Process Engineering, University of Stuttgart, 70569 Stuttgart, Germany http://orcid.org/0000-0003-4366-6120 http://crossmark.crossref.org/dialog/?doi=10.1007/s10450-022-00356-w&domain=pdf 126 Adsorption (2022) 28:125–136 1 3 measurements allows to examine in detail fundamental diffu- sion processes within nanoporous solids and allows to better understand nano-confinement effects [24]. The investigation of liquid phase systems was largely driven by studies related to liquid chromatography, i.e., partitioning behavior of (mul- ticomponent) mixtures at a solid surface functionalized for example with alkyl chains [25–29]. Other studies used coarse- grained molecular dynamics (MD) simulations to investigate protein adsorption on solid supports [30, 31] including the calculation of adsorption isotherms. An atomistic simulation study of enantiomeric separation of (R)- and (S)-ibuprofen in methanol solvent by means of immobilized cyclodextrins in a slit shaped model pore was reported recently [32]. Therefore, molecular simulation has evolved toward a promising tool to study liquid phase separation processes that are too complex to be described by phenomenological models [33]. The aim of the present work is to investigate the feasibility of all-atom classical MD simulations to reconcile liquid phase adsorption experiments with theoretical predictions. For this purpose adsorption of benzene and p-nitrophenol from aque- ous solutions onto cyclodextrin-functionalized mesoporous silica support is modelled and analyzed. Binding free enthalpy calculations in bulk solvent are related to the Henry regime of the adsorption isotherm on the functionalized material. The interpretation of experimentally observed adsorption iso- therms is discussed in view of the underlying molecular level picture. 2 Methodology 2.1 Calculation of binding free enthalpies and rate constants in bulk solution Two approaches were employed to calculate binding free enthalpies, an unbiased (counting) one and a biased (double decoupling) one. In the unbiased approach, referred to as direct count- ing (DC), the occurrences of bound, Nb , and unbound, Nu , instances during long ( t ≥ 10 μs ) standard MD simulations of one host–guest pair solvated in a box of water are counted. The binding free enthalpy is then obtained from [34] with standard state volume V0 = 1.661 nm3 and average vol- ume of the simulation box Vbox . In addition, average bound ⟨tb⟩ and unbound ⟨tu⟩ residence times can be calculated, yielding association kon and dissociation koff rates [34] (1)�GDC = −RT ln Nb Nu − RT ln Vbox V0 (2)kon = 1 ⟨tu⟩Cg , koff = 1 ⟨tb⟩ with guest molecule concentration Cg . The binding free enthalpy can then be determined using these rate constants (RC) with standard state concentration C0 = mol l−1. In order to identify bound and unbound instances, the host and guest molecule structures need to be geometrically reduced to comparable reference points. Using the different glucopyranose units, the conical shape of cyclodextrin was first reduced to three main circles, one running through the oxygen atoms of the primary hydroxyl groups, one central circle passing through the central carbon atoms, and one cir- cle connecting the oxygen atoms of the secondary hydroxyl groups. This system was then further reduced to a three point system based on the centers of mass of the different circles. Similarly, the two guest molecules were also reduced to a three point system as illustrated in Fig. 1. To determine whether a configuration is in a bound state, a spatial cut-off was defined for the central distance within which the guest molecule is assumed to be bound to the host molecule. Monitoring the angle of the two three point systems provides the orientation in which the guest molecule is bound to the host. Observing the bound and unbound state over time results into bound and unbound instances of dif- ferent duration. Averaging over these instances yields time averages ⟨tb⟩ and ⟨tu⟩ . Optionally, a temporal cut-off can be introduced that removes bound and unbound instances with a smaller residence time from the averaging. The impact of the temporal cut-off is shown in the Supplementary Material in Table S1 and discussed below. The second method used for calculating the binding free enthalpy was double decoupling (DD) [35]. As illustrated in Fig. 2 the process is divided into two parts. First, the free enthalpy difference �GM→M� u is calculated, resulting from decoupling the unbound state u, i.e., turning off the (3) �GRC = −RT ln konC 0 koff = −RT ln ⟨tb⟩ ⟨tu⟩ − RT ln Vbox V0 Fig. 1 Geometrical representation of the cyclodextrin and guest mol- ecule structures by a three point system 127Adsorption (2022) 28:125–136 1 3 intermolecular interactions with the environment ( M → M� ) of the guest molecule in a box of water while preserving intramolecular interactions. �GM→M� u is equal to the negative hydration free enthalpy −�GM hyd . Second, the free enthalpy difference �GM→M� b is calculated by decoupling the bound state (b), i.e., turning off intermolecular interactions of the guest molecule with the environment in a simulation with a host–guest complex in solvent. The latter is divided into three contributions, the difference �GM b→tor by turning on translational and orientational restraints (tor) between host and guest in order to guarantee that the guest molecule stays within the host when turning off intermolecular interactions �GM→M� tor and lastly turning off the initial restraints �GM� tor . The first two free enthalpy differences can be determined through molecular dynamics simulations by gradually turning on restraints and then turning off the interactions. This was done in one continuous simulation resulting into a combined value �GM→M� b→tor = �GM b→tor + �GM→M� tor . The free enthalpy from turning off the restraints from the non-interacting guest molecule can be calculated analytically. According to the thermodynamic cycle [36, 37] shown in Fig. 2, the summa- tion over the whole cycle is equal to zero thus resulting into (4) �GDD = �GM u→b = −�GM hyd − �GM b→tor − �GM→M� tor − �GM� tor = −�GM hyd − �GM→M� b→tor − �GM� tor . The procedure of decoupling the bound ligand was adopted from Boresch and Karplus [38]. Six atoms were chosen, three from the host a, b, c and three from the guest molecule A, B, C, see Fig. S2 in the Supplementary Material. Using these atoms a total of six harmonic restraints have been applied, a distance raA,0 , two angles �A,0 , �B,0 , and three dihedral angles �A,0 , �B,0 , �C,0 . The values for the reference distance and angles have been determined by averaging the distances and angles of host–guest complexes throughout the unbiased simulations. The analytical part for turning off these restraints is then given by with force constants Ki . Depending on the guest molecule, distinct binding configurations may become apparent. In this case the decoupling has to be performed for each of those states k. The total binding free enthalpy is then calculated by a logarithmic mean [39] 2.2 Immobilization of cyclodextrins For immobilizing cyclodextrin in a silica-pore, two linker concepts used by Huq et al. [17] and Trofymchuk et al. [20] were utilized for generating the molecule structure for simulation. Since the native cyclodextrin molecules were described by the Amber-compatible q4md-CD force field [40], the linker molecules were parametrized via Amber- Tools20 [41]. The parametrized linker structures were then appended to the cyclodextrin topology while accounting for additional connectivity parameters which are listed in the Supplementary Material in Table S3. The molecular struc- tures of the linkers are illustrated in Fig. 3. 2.3 Simulated systems For calculating the binding free enthalpy by direct counting and via rate constants, long unbiased simulations were con- ducted in the NpT ensemble with one guest molecule, ben- zene or p-nitrophenol, respectively, and one host molecule, β-cyclodextrin, solvated in 1000 water molecules. Simula- tions for determining the binding free enthalpy by double decoupling were initialized from a configuration of the long unbiased simulations containing a host–guest complex in the bound state. The hydration free enthalpy was calculated from a similar system without a host molecule. The binding (5) �GM� tor = −RT ln ⎡ ⎢⎢⎣ 8�V0(KrK�A K�B K�A K�B K�C ) 1 2 r2 aA,0 sin �A,0 sin �B,0(2�RT) 3 ⎤ ⎥⎥⎦ (6)𝛥ḠM u→b = −RT ln [∑ k exp ( 𝛥GM u→b,k RT )] . Fig. 2 Thermodynamic cycle for calculating the binding free enthalpy of a host–guest system. Hereby M indicates that the guest molecule is possessing full intermolecular (environmental) interac- tion (full triangle), while M′ denotes the guest molecule possessing only intramolecular interactions (framed triangle). Index u stands for an unbound state, b for the bound state, and tor for translational and orientational restraints 128 Adsorption (2022) 28:125–136 1 3 free enthalpy was calculated for two temperatures, 298 K and 350 K at a constant pressure of 1 bar. In order to generate and functionalize silica-pores as shown in Fig. 4, the PoreMS Python package [42] version 0.2.2 [43] was utilized. The systems are composed of a cylindrical pore of 4.05 nm diameter carved out of a (6.07, 6.14, 10.08) nm (x, y, z) β-cristobalite block. A bulk reser- voir with the length of 5 nm was attached on each side of the pore structure. The internal surface was functionalized with 0.07 and 0.14 μmol m−2 β-cyclodextrin, respectively. For the simulations representing the system by Trofymchuk et al. [20], additional 0.11 and 0.22 μmol m−2 of the propylamine group was added to the internal surface, which was experi- mentally used to immobilize cyclodextrin. This resulted into a total hydroxylation density (OH-groups on the surface) of 8.56 μmol m−2 on the internal surface and 8.82 μmol m−2 on the external surface. Further properties of the pore are listed Table 1. The topology parameters for the lattice silicon and oxy- gen atoms and silanol groups are taken from Gulmen and Thompson [44] and are shown in the Supplementary Mate- rial in Table S4. These systems were simulated in the NVT ensemble by achieving the desired density and pressure in a system with a specific concentration of guest molecules, through iteratively filling the simulation box with solute molecules until the difference between the density within the bulk-reservoir of the pore system and a preliminary NpT simulation at a pressure of 1 bar with the same solvent concentration, is less than one percent. Benzene adsorption isotherms were determined for two temperatures, 298K and 350K and two cyclodextrin surface densities. For p-nitro- phenol an adsorption isotherm was calculated at 350K using the pore with the higher cyclodextrin surface density. Na+ Fig. 3 Atom names and partial atomic charges (purple) of cyclodex- trin (red) with a pore surface linker (blue) used by a Huq et al. [17] and referred to as L1 in the present work, with a total charge beyond the Si1-atom (not counting the Na+ ) of − 1.32e and by b Trofym- chuk et  al. [20], referred to as L2 in the present work, with a total charge of − 0.32e (Color figure online) 6.07 nm 4.05 nm z x 5 nm 5 nm10.08 nm y x Fig. 4 Pore system functionalized with β-cyclodextrin using the sur- face linker used by Trofymchuk et al. [20]. a Side view of the sim- ulation box indicating the length of the central silica block and the solvent reservoirs. b Front view of the pierced silica block contain- ing the 4 nm diameter pore. The chemistry of the exterior surface is based on the (111) face of β-cristobalite silica. Bonded-phase groups are randomly distributed on the silica surface. Ligand densities, resid- ual surface hydroxylation, and further details are specified in Table 1. Colour code: Si, yellow lines; O, red lines; β-cyclodextrin, blue; pro- pylamine groups, magenta; surface silanol groups, yellow (Color fig- ure online) 129Adsorption (2022) 28:125–136 1 3 2.4 Simulation parameters Simulations were prepared using the open source package PoreSim [45] which generates folder structures and other practical scripts for pore simulations. The simulation suite GROMACS 2016.5 [46, 47] was used for all simulations, while PLUMED 2.5 [48, 49] was utilized to extract spe- cific distances and angles. Based on earlier work [50, 51], the general Amber force field (GAFF) [52] was chosen to describe intra- and intermolecular interactions. This is further backed up by the quality of the GAFF-compatible force field q4md-CD for cyclodextrin simulations [40, 53]. In order to verify the quality of the force-fields for the solute molecules p-nitrophenol and benzene, which were studied in the work of Huq et al. [17] and Trofymchuk et al. [20], validating simulations were carried out based on topolo- gies provided by Mobley et al. [54]. Water was described by the TIP3P model [55] as the relatively large pore diamter of 4.05 nm is not expected to lead to a strong confinement effect at ambient pressure [56]. All MD simulations were performed under periodic boundary conditions. Temperature was controlled using the Nosé–Hoover thermostat [57, 58] with a coupling constant of 1 ps, while pressure for simulations in the NpT ensemble was controlled by the Parrinello-Rahman barostat [59, 60] with a coupling constant of 5.0 ps and compressibility of 4.5 × 10−5 bar−1 . Bond lenghts between heavy atoms and hydrogens were constrained with the LINCS algorithm [61, 62] with an order of 4. Short-range electrostatic and Len- nard–Jones parameters were evaluated up to a cutoff distance of 1.4 nm. Analytical dispersion corrections for energy and pressure were included. Long-range electrostatic interactions were treated with the particle-mesh Ewald algorithm [63, 64]. The long unbiased simulations were run for 10 μs with a time-step of 2 fs after a total equilibration time of 3 ns. Decoupling simulations used a total of 25 intermediate states ( �-points) for a slow equidistant deactivation of the intermo- lecular interactions of the guest molecule with its environ- ment while preserving intramolecular interactions. Each � -point was run for 1 ns with a time-step of 2 fs. During the simulation of the pore systems, silicon and oxygen grid atoms were frozen in their position to preserve the original pore shape, this includes the silicon atom of surface groups. Table 1 Properties of the cylindrical mesopore model before and after functionalization (func.), generated as described by Trofymchuk et al. [20], with surface densities ( μmol m−2 ) and pore dimensions (nm) a Calculated as the standard deviation of the shortest distances between the central pore axis and the surface Si atoms b Pore variants with 5 and 11 β-cyclodextrin groups were prepared in this work Interior Exterior Silica block xyz-dim. 6.07; 6.14; 10.08 Pore drilling direction z Pore diameter 4.05 Surface roughnessa 0.08 0.00 Solvent reservoir z-dim. 5.00 Simulation box xyz-dim. 6.07; 6.14; 20.08 Pore volume 129.74 Solvent reservoir volume 2 × 186.36 Surface area 128.16 2 × 24.39 Surface chemistry before func.  Num. of single silanol groups 553 205  Num. of geminal silanol groups 68 27  Num. of siloxane bridges 0 0  Total number of OH groups 689 259  Overall hydroxylation 8.93 8.82 Surface chemistry after func.  Num. of β-cyclodextrin groupsb 11 0  β-cyclodextrin density 0.14 0  Num. of propylamine groups 17 0  Propylamine density 0.22 0  Bonded-phase density 0.36 0.00  Num. of residual OH groups 661 259  Residual hydroxylation 8.56 8.82 130 Adsorption (2022) 28:125–136 1 3 For these systems a trajectory length of 1 μs was generated with a time-step of 1 fs and a total equilibration time of 100 ns. 2.5 Analysis Distances and angles between the reference systems dur- ing the long unbiased simulations were extracted using PLUMED. The determination of bound and unbound states, and the calculation of the binding free enthalpy using the direct counting and rate constants method was conducted by in-house Python scripts. For determining the free enthalpy differences from decoupling simulations thermodynamic integration [65] was utilized. The density profiles and adsorption isotherms were calculated using the PoreAna package version 0.2.0 [66], developed dur- ing this project in object oriented Python 3 to complement the PoreMS package. Radial density and diffusion profiles within the pore were calculated by dividing the cylindrical shape into radial volume slices with pore length z and radius ri of slice i. The number density �i is then determined by counting the number of molecules Ni within each slice during the simulation, resulting into Adsorption isotherms describe the amount of solute mol- ecules adsorbed on the surface Npore as a function of the amount of solute in the bulk phase Nbulk and are therefore, similarly to the density, determined from counting the num- ber of molecules within the pore and within the bulk reser- voirs normalized by the number of frames with the number of molecules Nj within the pore or bulk phase at frame j = 1,… ,NF . These values are then con- verted to surface and volume concentrations respectively based on the volume of the pore system and inner pore surface. The diffusion coefficient D‖,i parallel to the pore surface was determined from the slope of the mean square dis- placement �i(t) over an observation time of 4-20 ps within each bin with a tolerance of ±1 bins (7)Vi = �z(r2 i − r2 i−1 ) (8)�i = Ni Vi = Ni �z 1 r2 i − r2 i−1 . (9) ⟨N⟩pore = 1 NF NF� j=1 N pore j , ⟨N⟩bulk = 1 NF NF� j=1 Nbulk j . By weighting the axial diffusion profile D‖,i with the density profile �i along the radius r, a mean diffusion coefficient ⟨D‖⟩ can be calculated with the cross-sectional area of bin i. 2.6 Langmuir model Considering a simulation set-up as shown in Fig. 4 but with only a single cyclodextrin molecule bound to the inner pore surface and a single guest molecule present in the simula- tion box, the direct counting method (Eq. 1) would result in the same standard binding free enthalpy as in the bulk phase simulation given that the pore walls do not interact substantially with the guest molecule. Therefore, the ratio of bound to unbound samples in the pore system would be where the sub- or superscript ’pore’ refers to the entire accessible volume of the simulation box, containing the pore and the solvent reservoirs. Rearranging to and replacing the two rightmost terms by means of Eq. 1 results in The first term on the right-hand side can be identified with the bulk concentration of the guest molecule such that the equation has the form of the Henry isotherm. If we assume that each cyclodextrin can only accommodate one guest molecule and relate the amount adsorbed to the inner pore surface we obtain the Langmuir form [67] where qmax denotes the cyclodextrin density inside the pore and (10)D‖,i = 1 2 d�i(t) dt . (11)⟨D‖⟩ = ∑ �iD‖,iAi∑ �iAi = ∑ �iD‖,i(r2i − r2 i−1 ) ∑ �i(r 2 i − r2 i−1 ) (12)Ai = �(r2 i − r2 i−1 ) (13) ( Nb Nu ) pore = ( Nb Nu ) bulk ⋅ Vbulk box V pore box (14)N pore b = ( Nu Vbox ) pore ⋅ ( Nb Nu ) bulk ⋅ Vbulk box (15)N pore b = ( Nu Vbox ) pore ⋅ V0 ⋅ exp ( − �GDC RT ) . (16)qads = qmax K ⋅ C 1 + K ⋅ C 131Adsorption (2022) 28:125–136 1 3 with C0 as the concentration of the standard state, i.e. 1mol l−1 and �G as the standard binding free enthalpy obtained via double decoupling, direct counting or any other suitable computational or experimental approach. In that way, a Langmuir isotherm can be computed and compared to isotherms obtained from experiments or molecular simu- lation to assess whether other processes such as binding of multiple guest molecules to one cyclodextrin or cooperative effects of cyclodextrin molecules in close vicinity on the surface are likely to occur. 3 Results and discussion 3.1 Bulk phase simulations In order to verify the quality of the force field for the sol- ute molecules p-nitrophenol and benzene, validating simu- lations were carried out based on topologies provided by Mobley et al. [54] resulting in liquid density values in good agreement with experiment, see Table 2. Hydration free enthalpies �Ghyd reported by the Mobley group [68] were reproduced by decoupling simulations of these molecules in water, which in turn are in fair (p-nitrophenol) or good (benzene) agreement with the experimental values. Table 3 shows the binding free enthalpies and rate con- stants of CD + p-nitrophenol and CD + benzene complexes. For the long unbiased simulations first, a spatial cutoff had to be defined in order to differentiate between bound and unbound states. Since β-Cyclodextrin has a radius of gyra- tion around 0.6 nm [71], the spatial cutoff distance was cho- sen at 0.7 nm to account for binding on the inner edge of the host molecule. This assumption is further strengthened by Fig. 5 where the histogram maximum of the center of mass distances between host and guest molecules diminishes at the chosen cutoff and a smaller cluster emerges that indicates the unbound states. The effect of the temporal cut-off between 0 and 1 ns on the binding free enthalpy is summarized in Table S1 in (17)K = 1 C0 exp ( − �G RT ) the Supplementary Material. The consideration behind the temporal cut-off is that the average residence time of the guest molecule inside the host is usually shorter than an average experimentally determined survival time because in an experiment several short-term events are likely to be missed. This effect has to be kept in mind when comparing to experimentally determined rates which may differ depend- ing on the temporal resolution of the measurement device. For benzene the effect on �G is minor, while for p-nitrophe- nol an effect is only visible between temporal cut-offs of 0 ps and 100 ps. For both molecules the kon - and koff-rates calculated with 0 ps cut-off are an order of magnitude larger compared to the values calculated with a larger cut-off. In the present work only bound/unbound periods that lasted longer than 1 ns were counted as one bound/unbound event, in agreement with previous works [34, 72]. For CD + benzene all three approaches used to calculate the binding free enthalpy yield a consistent value at both temperatures. For CD + p-nitrophenol the counting approaches did not provide reliable values at 298K due to the rather strong binding and thus low occurrences of unbound instances. The time evolution of bound and unbound instances are provided in the Supplementary Mate- rial in Fig. S1. Therefore only the double decoupling results are reported at this temperature. At 350K unbound occur- rences are more likely due to the higher temperature, allow- ing an improved sampling which leads to results for the direct counting approach that is in good agreement with the double decoupling method. Nevertheless, the drop of the binding free enthalpy in the rate constants approach indi- cates that longer unbound instances are still infrequent. The reason for the larger binding free enthalpy at a higher tem- perature for benzene in the unbiased approaches is due to a Table 2 Liquid densities � ( kgm−3 ) for p-nitrophenol (p-NP) and benzene (BEN) compared with experimental data [69] Determined at a387K , b298K . Hydration free enthalpy �Ghyd ( kJmol−1 ) at 298K was calculated via simulation (sim) and compared to simulation values reported by the Mobley group [68] (mo) and experimental values (exp) from Cabani et al. [70] �sim �exp �Gsim hyd �Gmo hyd �G exp hyd p-NP 1274.75a 1283.15a − 35.35 − 35.44 − 44.08 BEN 861.24b 873.49b − 3.23 − 3.39 − 3.62 Fig. 5 Histogram of the central distances C-C with a marked area of the bound states (grey) of the reference systems for host and guest molecules as illustrated in Fig. 1 for benzene (blue) and p-nitrophenol (orange) at T = 350K with the corresponding image showing ben- zene in the bound state (Color figure online) 132 Adsorption (2022) 28:125–136 1 3 larger volumetric correction in Eq. ( 1) at RT 298 ln Vbox V0 = 7.28 kJ mol −1 and, RT350 ln Vbox V0 = 8.69 kJ mol−1 , which is not entirely compensated by the decrease of the ratio of bound and unbound instances RT298 ln Nb Nu = 5.98 kJ mol−1 and, RT350 ln Nb Nu = 5.41 kJ mol−1 , as they stay almost similar. For assessing the properties of the parametrized CD molecule with an attached linker for connecting to the pore surface, NpT simulations in TIP3P water were carried out and the resulting dihedral-angle distributions were analyzed. Figure 6a shows that for the native cyclodextrin the distribu- tions are in good agreement with those reported by Gebhardt et al. [53] using the same force field. Additional dihedral- angle distributions for α - and γ-cyclodextrin are shown in the Supplementary Material in Fig. S3. Attaching the linker affects the corresponding dihedral angle describing rotation around the C6–O6 bond, see Fig. 3. With the L1-variant essentially one rotamer around 60◦  is populated (Fig. 6b) while the L2-variant has two rotamers populated (Fig. 6c). 3.2 Unfunctionalized silica pores Figure 7 shows density profiles of water, benzene, and p-nitrophenol as well as axial diffusion coefficient profiles in a system containing 60 solute molecules which corresponds to a concentration of 200mmol l−1 . Hardly any adsorption for benzene or p-nitrophenol is visible, in agreement with experiments [17, 20]. The number densities converges rapidly towards the density value at the pore center. The slowdown of water diffusion in confinement relative to the bulk phase can be compared with recent experimental data reported by Jani et al. [76]. At a temperature of 300K the experimental self-diffusion coefficient of water in the bulk phase, 2.3 × 10−9 m2 s−1 , is reduced to 2.0 × 10−9 m2 s−1 in SBA-15 with a pore mean pore diameter of 6.6 nm and to 1.7 × 10−9 m2 s−1 in MCM-41 with a pore diameter of 3.8 nm, close to the simulated pore diameter of 4 nm. This accounts to change of factor 1.15 for the SBA-15 system and 1.35 in the MCM-41 pore. The mean diffusion ⟨D‖⟩ of water in the simulated unfunctionalized pore is 3.07 × 10−9 m2 s−1 at a temperature of 298K and 5.63 × 10−9 m2 s−1 at a higher temperature of 350K . This yields a factor 1.70 at 298K and 1.74 at 350 K by comparing the pore diffusion to the bulk diffusion of the same system. 3.3 Cyclodextrin‑functionalized silica pores By functionalizing the surface, benzene molecules have a significantly higher concentration at the configura- tional space of the cyclodextrin center of mass, indicating host–guest interaction, see Fig. 8. Benzene adsorption isotherms were calculated for a low cyclodextrin surface concentration at 298K and a high cyclodextrin concentration at 298K and 350K . For p-nitro- phenol only the larger surface concentration at 350K was considered due to the strong binding affinity that leads to unreliable statistics at the lower temperature. The volumetric and the surface concentration of the guest molecules were assessed. The volumetric concentration in the reservoir was converted from the average number of molecules, which was determined by counting the solute molecule occurrences inside the reservoir normalized by the number of frames. Similarly, a surface specific concentration was determined by counting the appearances of the solute molecules inside the pore. In order to obtain excess adsorption, the adsorption value of the solute molecules within a non-functionalized pore has been subtracted. Repeating the pore simulation for Table 3 Comparison of direct counting (DC), rate constants (RC) and double decoupling (DD) methods for determining the binding free enthalpy �G ( kJmol−1 ) with experimental (Exp) data of β-cyclodextrin-ligand complexes, with association kon ( 108 dm 3 mol−1 s−1 ) and dissociation koff ( 106 s−1 ) rates at a temporal cut-off of 1 ns calculated at temperature T ( K) Mean value from aDC, DD, and bDC, RC, DD. The chosen spatial cut-off for differentiating between bound and unbound states is 0.7 nm for both guest molecules p-nitrophenol (p-NP) and benzene (BEN) T 298 350 −�G kon koff −�G kon koff p-NP  DC 24.22  RC 17.01 11.6 3.34  DD 25.79 24.07  Meana 24.15  Exp [73, 74] 10–17 BEN  DC 13.25 14.11  RC 13.27 5.41 2.55 13.46 8.58 8.41  DD 13.12 12.15  Meanb 13.21 13.24  Exp [75] 12–13 133Adsorption (2022) 28:125–136 1 3 different solute concentrations in the system, results into adsorption isotherms shown in Fig. 9. Further isotherms are provided in Fig. S4 of the Supplementary Material. Similar to experimental observation a simple relation between the amount adsorbed and the number of cyclo- dextrin molecules attached to the surface could only be observed at small concentration, following Langmuir behavior. For benzene, Trofymchuk et al. [20] reported adsorption isotherms for materials with different amounts of cyclodextrin molecules up to benzene concentrations of 7mmol l−1 , i.e., roughly one third of the solubility limit of 23.8mmol l−1 at 25 °C [77]. The corresponding amount adsorbed exceeded the capacity of a 1:1 binding by more than a factor of 10 and the isotherms showed dual-site Langmuir type behavior. In the present work higher benzene concentrations of up to 60mmol l−1 were used to improve the statistical sampling. However, the amount adsorbed did not exceed the 1:1 binding capac- ity by more than a factor of two. The visual inspection of the trajectories shows that for the lower cyclodextrin density (i.e. five molecules attached to the pore surface) up to three benzene molecules may be associated to one cyclodextrin molecule. For the higher cyclodextrin density some cyclodextrin molecules may also be trapped in the space between the pore surface and the outer surface of the cyclodextrin, thus enhancing the adsorption capacity. Fig. 6 a Distributions of the seven individual dihedral angles H6O– O6–C6–H61 for the native β-cyclodextrin resulting from an NpT simulation in TIP3P water. The different colours denote the different glucopyranose units. b Distribution of six dihedral angles H6O–O6– C6–H61 for β-cyclodextrin connected with the L1 linker [17]. The blue line represents the distribution of the C5(tail)–O6–C6–H61 dihe- dral angle. c Analogous to b but for the L2 linker [20] with the blue line representing the distribution of the C4(tail)–O6–C6–H61 dihe- dral angle (Color figure online) Fig. 7 a Density and b axial diffusion profiles of water (blue), ben- zene (BEN) (orange), and p-nitrophenol (p-NP) (green) in a non- functionalized pore at temperatures 298K (dashed lines) and 350K (solid lines). The shaded area denotes the configurational space of the silanol group oxygen atoms (yellow) (Color figure online) Fig. 8 Radial density profile of the benzene centre of mass in a non- functionalized pore (dashed lines) and in a functionalized pore using the L2-variant cyclodextrin (L2-CD) [20] (solid lines) simulated for 1 μs . Shaded areas denote the configurational space of the cyclodex- trin centre of mass (blue) and the silanol groups oxygen atoms (yel- low) (Color figure online) 134 Adsorption (2022) 28:125–136 1 3 Moreover, cyclodextrin molecules in close vicinity may associate and form rather long-living complexes that encapsulate the solute molecules. In addition, a benzene molecule dissociating from one cyclodextrin is very likely associating with the next one close by instead of leaving the pore. For p-nitrophenol experimental adsorption isotherms were reported for very low bulk concentrations below 0.15 mM, compared to the aqueous solubility of 115 mM at 25 °C [78]. In this regime Langmuir behavior was observed with the Langmuir constant resembling the 1:1 associa- tion free enthalpy, as expected [18]. Shen et al. performed adsorption measurements at higher initial concentration of up to 28 mM p-nitrophenol and found adsorption capaci- ties significantly larger than those corresponding to a 1:1 binding scenario and explained this finding with hydro- gen bonds that are formed between the polar groups of p-nitrophenol and the hydroxyl groups of cyclodextrin and the amine groups of the functionalized silica, respectively [79]. In the simulations, a large range of bulk-phase con- centration has been studied. The binding free enthalpy for the 1:1 complex is overestimated by the force field, leading to rather strong increase in the amount adsorbed at low bulk phase concentration that first follows the Langmuir isotherm but exceeds the 1:1 binding capacity by more than a factor of three due to the association of up to three nitrophenol molecules with one cyclodextrin molecule, the trapping between cyclodextrin and the pore surface and the formation of hydrogen bonds between the hydroxyl groups at the rims of cyclodextrin and p-nitrophenol in the solvent phase. 4 Conclusion and outlook An efficient model building is a prerequisite for computa- tional studies of functionalized mesoporous silica materi- als. The Python package PoreMS introduced previously [42, 43] was complemented by two additional program packages for preparing MD simulations of porous materi- als with GROMACS, PoreSim [45], and for analyzing the simulation trajectories, PoreAna [66], the latter provid- ing results such as density and diffusion profiles, thereby reducing the overhead for system preparation to analysis. The selected case study of adsorption of aromatic mol- ecules in cyclodextrin-functionalized silica mesopores shows that current moderate computational resources allow an atomistically resolved model (~ 65,000 atoms) to be propagated to the μs time scale. The calculated adsorp- tion isotherms show a more complex behavior than pre- dicted by a simple Langmuir model corresponding to a 1:1 host:guest binding complex. Beside the formation of 1:2 and even 1:3 host:guest complexes also host–host interac- tions on the surface as well as more unspecific host–guest interactions have an influence on the shape of the iso- therm. The information is relevant, for example, for the prediction of band broadening in liquid chromatography [80]. While the aim of the present work was to investi- gate the feasibility of calculating liquid phase adsorption isotherms by all-atom MD, future work may address the study of multicomponent mixtures, the influence of the solvent or the investigation of chiral stationary phases, e.g by attaching cyclodextrin derivatives. Fig. 9 a Excess adsorption isotherms of pore simulations utilizing an L2-variant [20] cyclodextrin functionalized surface (blue) compared to the analytical solution shown in equation (16) (orange) at varying amounts of benzene. The binding free enthalpy used to evaluate Eq. (16) was calculated as the mean value of all utilized methods (DC, RC, DD) at the respective temperature. b Like a but with an L1-var- iant cyclodextrin [17] with p-nitrophenol and the mean free enthalpy value only from the DD and DC approach. Simulation were run for 1 μs at 298K (dashed lines) and 350K (solid lines). The grey line rep- resents a hypothetical maximum of 1:1 binding with 11 cyclodextrin molecules (Color figure online) 135Adsorption (2022) 28:125–136 1 3 Supplementary Information The online version contains supplemen- tary material available at https:// doi. org/ 10. 1007/ s10450- 022- 00356-w. Acknowledgements This work was financially supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foun- dation) - project-ID 358283783 - SFB 1333 and by the Ministry of Science, Research and the Arts of the State of Baden-Württemberg, Germany. Simulation were performed on the supercomputer ForHLR 2 and its successor HoreKa funded by the Ministry of Science, Research and the Arts Baden-Württemberg and by the Federal Ministry of Edu- cation and Research. Funding Open Access funding enabled and organized by Projekt DEAL. Declarations Conflict of interest The authors declare that they have no conflict of interest. Open Access This article is licensed under a Creative Commons Attri- bution 4.0 International License, which permits use, sharing, adapta- tion, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. References 1. Sawicki, R., Mercier, L.: Environ. Sci. Technol. 40(6), 1978 (2006). https:// doi. org/ 10. 1021/ es051 441r 2. Tian, B., Hua, S., Tian, Y., Liu, J.: Environ. Sci. Pollut. Res. 28(2), 1317 (2021). https:// doi. org/ 10. 1007/ s11356- 020- 11168-2 3. Lorenz, H., Seidel-Morgenstern, A.: Angew. Chem. Int. Ed. 53(5), 1218 (2014). https:// doi. org/ 10. 1002/ anie. 20130 2823 4. Breslow, R., Dong, S.D.: Chem. Rev. 98(5), 1997 (1998). https:// doi. org/ 10. 1021/ cr970 011j 5. Hu, S., Li, J., Xiang, J., Pan, J., Luo, S., Cheng, J.P.: J. Am. Chem. Soc. 132(20), 7216 (2010). https:// doi. org/ 10. 1021/ ja102 819g 6. van der Helm, M.P., Klemm, B., Eelkema, R.: Nat. Rev. Chem. 3(8), 491 (2019). https:// doi. org/ 10. 1038/ s41570- 019- 0116-0 7. De Rosa, M., La Manna, P., Talotta, C., Soriente, A., Gaeta, C., Neri, P.: Front. Chem. 6, 84 (2018). https:// doi. org/ 10. 3389/ fchem. 2018. 00084 8. Rai, P., Srivastava, M., Yadav, S., Singh, J., Singh, J.: Catal. Lett. 145(12), 2020 (2015). https:// doi. org/ 10. 1007/ s10562- 015- 1588-2 9. Molnár, A.: ChemCatChem 13(6), 1424 (2021). https:// doi. org/ 10. 1002/ cctc. 20200 1610 10. Gaeta, C., La Manna, P., De Rosa, M., Soriente, A., Talotta, C., Neri, P.: ChemCatChem 13(7), 1638 (2021). https:// doi. org/ 10. 1002/ cctc. 20200 1570 11. Del Valle, E.: Process Biochem. 39(9), 1033 (2004). https:// doi. org/ 10. 1016/ S0032- 9592(03) 00258-9 12. Schmidt, B.V.K.J., Barner-Kowollik, C.: Angew. Chem. Int. Ed. 56(29), 8350 (2017). https:// doi. org/ 10. 1002/ anie. 20161 2150 13. Gattuso, G., Nepogodiev, S.A., Stoddart, J.F.: Chem. Rev. 98(5), 1919 (1998). https:// doi. org/ 10. 1021/ cr960 133w 14. Khan, A.R., Forgo, P., Stine, K.J., D’Souza, V.T.: Chem. Rev. 98(5), 1977 (1998). https:// doi. org/ 10. 1021/ cr970 012b 15. Harada, A., Li, J., Kamachi, M.: Nature 364(6437), 516 (1993). https:// doi. org/ 10. 1038/ 36451 6a0 16. Ciucanu, I., König, W.A.: J. Chromatogr. A 685(1), 166 (1994). https:// doi. org/ 10. 1016/ 0021- 9673(94) 00564-8 17. Huq, R., Mercier, L., Kooyman, P.J.: Chem. Mater. 13(12), 4512 (2001). https:// doi. org/ 10. 1021/ cm010 171i 18. Bibby, A., Mercier, L.: Green Chem. 5, 15 (2003). https:// doi. org/ 10. 1039/ B2092 51B 19. Roik, N.V., Belyakova, L.A., Trofymchuk, I.M., Dziazko, M.O., Oranska, O.I.: J. Nanopart. Res. 19(9), 317 (2017). https:// doi. org/ 10. 1007/ s11051- 017- 3999-z 20. Trofymchuk, I., Roik, N., Belyakova, L.: Nanoscale Res. Lett. 12(1), 307 (2017). https:// doi. org/ 10. 1186/ s11671- 017- 2072-2 21. Derrah, R.I., Ruthven, D.M.: Can. J. Chem. 53(7), 996 (1975). https:// doi. org/ 10. 1139/ v75- 142 22. Smit, B., Maesen, T.L.M.: Chem. Rev. 108(10), 4125 (2008). https:// doi. org/ 10. 1021/ cr800 2642 23. Da̧browski, A.: Adv. Colloid. Interface Sci. 93(1–3), 135 (2001). https:// doi. org/ 10. 1016/ S0001- 8686(00) 00082-8 24. Bukowski, B.C., Keil, F.J., Ravikovitch, P.I., Sastre, G., Snurr, R.Q., Coppens, M.O.: Adsorption 27, 683 (2021). https:// doi. org/ 10. 1007/ s10450- 021- 00314-y 25. Lindsey, R.K., Rafferty, J.L., Eggimann, B.L., Siepmann, J.I., Schure, M.R.: J. Chromatogr. A 1287, 60 (2013). https:// doi. org/ 10. 1016/j. chroma. 2013. 02. 040 26. Melnikov, S.M., Höltzel, A., Seidel-Morgenstern, A., Tallarek, U.: Angew. Chem. Int. Ed. 51(25), 6251 (2012). https:// doi. org/ 10. 1002/ anie. 20120 1096 27. Rybka, J., Höltzel, A., Tallarek, U.: J. Phys. Chem. C 121(33), 17907 (2017). https:// doi. org/ 10. 1021/ acs. jpcc. 7b047 46 28. El Hage, K., Bemish, R.J., Meuwly, M.: Phys. Chem. Chem. Phys. 20(27), 18610 (2018). https:// doi. org/ 10. 1039/ C8CP0 2899K 29. Trebel, N., Höltzel, A., Lutz, J.K., Tallarek, U.: J. Phys. Chem. B 125(40), 11320 (2021). https:// doi. org/ 10. 1021/ acs. jpcb. 1c067 32 30. Basconi, J.E., Carta, G., Shirts, M.R.: AIChE J. 63, 4564 (2017). https:// doi. org/ 10. 1002/ aic. 15798 31. Jakobtorweihen, S., Heuer, J., Waluga, T.: J. Chromatogr. A 1620, 460940 (2020). https:// doi. org/ 10. 1016/j. chroma. 2020. 460940 32. Nezbeda, I., Škvára, J.: Mol. Simul. 47, 846 (2020). https:// doi. org/ 10. 1080/ 08927 022. 2020. 18285 84 33. Gritti, F.: Anal. Chem. 93, 5653 (2021). https:// doi. org/ 10. 1021/ acs. analc hem. 0c050 78 34. Baz, J., Gebhardt, J., Kraus, H., Markthaler, D., Hansen, N.: Chem. Ing. Tech. 90(11), 1864 (2018). https:// doi. org/ 10. 1002/ cite. 20180 0050 35. Gilson, M., Given, J., Bush, B., McCammon, J.: Biophys. J. 72(3), 1047 (1997). https:// doi. org/ 10. 1016/ S0006- 3495(97) 78756-3 36. van Gunsteren, W.F., Berendsen, H.J.C.: J. Comput. Aid. Mol. Des. 1(2), 171 (1987). https:// doi. org/ 10. 1007/ BF016 76960 37. Mugnai, M.L., Elber, R.: J. Chem. Theory Comput. 8(9), 3022 (2012). https:// doi. org/ 10. 1021/ ct300 3817 38. Boresch, S., Tettinger, F., Leitgeb, M., Karplus, M.: J. Phys. Chem. B 107(35), 9535 (2003). https:// doi. org/ 10. 1021/ jp021 7839 39. Mobley, D.L., Chodera, J.D., Dill, K.A.: J. Chem. Phys. 125(8), 084902 (2006). https:// doi. org/ 10. 1063/1. 22216 83 40. Cézard, C., Trivelli, X., Aubry, F., Djedaïni-Pilard, F., Dupradeau, F.Y.: Phys. Chem. Chem. Phys. 13(33), 15103 (2011). https:// doi. org/ 10. 1039/ c1cp2 0854c 41. Case, D.A., Belfon, K., Ben-Shalom, I.Y., Brozell, S.R., Cerutti, D.S., Cheatham, T.E., III., Cruzeiro, V.W.D., Darden, T.A., Duke, R.E., Giambasu, G., Gilson, M.K., Gohlke, H., Goetz, A.W., Har- ris, R., Izadi, S., Izmailov, S.A., Kasavajhala, K., Kovalenko, A., https://doi.org/10.1007/s10450-022-00356-w http://creativecommons.org/licenses/by/4.0/ https://doi.org/10.1021/es051441r https://doi.org/10.1007/s11356-020-11168-2 https://doi.org/10.1002/anie.201302823 https://doi.org/10.1021/cr970011j https://doi.org/10.1021/cr970011j https://doi.org/10.1021/ja102819g https://doi.org/10.1038/s41570-019-0116-0 https://doi.org/10.3389/fchem.2018.00084 https://doi.org/10.3389/fchem.2018.00084 https://doi.org/10.1007/s10562-015-1588-2 https://doi.org/10.1002/cctc.202001610 https://doi.org/10.1002/cctc.202001610 https://doi.org/10.1002/cctc.202001570 https://doi.org/10.1002/cctc.202001570 https://doi.org/10.1016/S0032-9592(03)00258-9 https://doi.org/10.1016/S0032-9592(03)00258-9 https://doi.org/10.1002/anie.201612150 https://doi.org/10.1021/cr960133w https://doi.org/10.1021/cr970012b https://doi.org/10.1038/364516a0 https://doi.org/10.1016/0021-9673(94)00564-8 https://doi.org/10.1021/cm010171i https://doi.org/10.1039/B209251B https://doi.org/10.1039/B209251B https://doi.org/10.1007/s11051-017-3999-z https://doi.org/10.1007/s11051-017-3999-z https://doi.org/10.1186/s11671-017-2072-2 https://doi.org/10.1139/v75-142 https://doi.org/10.1021/cr8002642 https://doi.org/10.1016/S0001-8686(00)00082-8 https://doi.org/10.1007/s10450-021-00314-y https://doi.org/10.1007/s10450-021-00314-y https://doi.org/10.1016/j.chroma.2013.02.040 https://doi.org/10.1016/j.chroma.2013.02.040 https://doi.org/10.1002/anie.201201096 https://doi.org/10.1002/anie.201201096 https://doi.org/10.1021/acs.jpcc.7b04746 https://doi.org/10.1039/C8CP02899K https://doi.org/10.1021/acs.jpcb.1c06732 https://doi.org/10.1002/aic.15798 https://doi.org/10.1016/j.chroma.2020.460940 https://doi.org/10.1080/08927022.2020.1828584 https://doi.org/10.1080/08927022.2020.1828584 https://doi.org/10.1021/acs.analchem.0c05078 https://doi.org/10.1021/acs.analchem.0c05078 https://doi.org/10.1002/cite.201800050 https://doi.org/10.1002/cite.201800050 https://doi.org/10.1016/S0006-3495(97)78756-3 https://doi.org/10.1007/BF01676960 https://doi.org/10.1021/ct3003817 https://doi.org/10.1021/jp0217839 https://doi.org/10.1063/1.2221683 https://doi.org/10.1039/c1cp20854c https://doi.org/10.1039/c1cp20854c 136 Adsorption (2022) 28:125–136 1 3 Krasny, R., Kurtzman, T., Lee, T.S., LeGrand, S., Li, P., Lin, C., Liu, J., Luchko, T., Luo, R., Man, V., Merz, K.M., Miao, Y., Mikhailovskii, O., Monard, G., Nguyen, H., Onufriev, A., Pan, F., Pantano, S., Qi, R., Roe, D.R., Roitberg, A., Sagui, C., Schott- Verdugo, S., Shen, J., Simmerling, C.L., Skrynnikov, N.R., Smith, J., Swails, J., Walker, R.C., Wang, J., Wilson, L., Wolf, R.M., Wu, X., Xiong, Y., Xue, Y., York, D.M., Kollman, P.A.: AMBER 2020 (2020) 42. Kraus, H., Rybka, J., Höltzel, A., Trebel, N., Tallarek, U., Hansen, N.: Mol. Simul. 47(4), 306 (2021). https:// doi. org/ 10. 1080/ 08927 022. 2020. 18714 78 43. Kraus, H., Hansen, N.: Porems: 0.2.2 (2021). https:// doi. org/ 10. 5281/ zenodo. 47380 36 44. Gulmen, T.S., Thompson, W.H.: Langmuir 22(26), 10919 (2006). https:// doi. org/ 10. 1021/ la062 285k 45. Kraus, H.: Poresim: 0.1.0 (2021). https:// doi. org/ 10. 5281/ zenodo. 45888 61 46. Abraham, M.J., Murtola, T., Schulz, R., Páll, S., Smith, J.C., Hess, B., Lindahl, E.: SoftwareX 1–2, 19 (2015). https:// doi. org/ 10. 1016/j. softx. 2015. 06. 001 47. Hess, B., Kutzner, C., van der Spoel, D., Lindahl, E.: J. Chem. Theory. Comput. 4(3), 435 (2008). https:// doi. org/ 10. 1021/ ct700 301q 48. Bonomi, M., Branduardi, D., Bussi, G., Camilloni, C., Provasi, D., Raiteri, P., Donadio, D., Marinelli, F., Pietrucci, F., Broglia, R.A., Parrinello, M.: Comput. Phys. Commun. 180(10), 1961 (2009). https:// doi. org/ 10. 1016/j. cpc. 2009. 05. 011 49. Tribello, G.A., Bonomi, M., Branduardi, D., Camilloni, C., Bussi, G.: Comput. Phys. Commun. 185(2), 604 (2014). https:// doi. org/ 10. 1016/j. cpc. 2013. 09. 018 50. Ziegler, F., Teske, J., Elser, I., Dyballa, M., Frey, W., Kraus, H., Hansen, N., Rybka, J., Tallarek, U., Buchmeiser, M.R.: J. Am. Chem. Soc. 141(48), 19014 (2019). https:// doi. org/ 10. 1021/ jacs. 9b087 76 51. Ziegler, F., Kraus, H., Benedikter, M.J., Wang, D., Bruckner, J.R., Nowakowski, M., Weißer, K., Solodenko, H., Schmitz, G., Bauer, M., Hansen, N., Buchmeiser, M.R.: ACS Catal. 11, 11570 (2021). https:// doi. org/ 10. 1021/ acsca tal. 1c030 57 52. Wang, J., Wolf, R.M., Caldwell, J.W., Kollman, P.A., Case, D.A.: J. Comput. Chem. 25(9), 1157 (2004). https:// doi. org/ 10. 1002/ jcc. 20035 53. Gebhardt, J., Kleist, C., Jakobtorweihen, S., Hansen, N.: J. Phys. Chem. B 122(5), 1608 (2018). https:// doi. org/ 10. 1021/ acs. jpcb. 7b118 08 54. Mobley, D.L., Shirts, M., Lim, N., Chodera, J., Beauchamp, K.: Lee-Ping. MobleyLab/FreeSolv: Version , 52 (2018). https:// doi. org/ 10. 5281/ zenodo. 11612 45 55. Jorgensen, W.L., Chandrasekhar, J., Madura, J.D., Impey, R.W., Klein, M.L.: J. Chem. Phys. 79(2), 926 (1983). https:// doi. org/ 10. 1063/1. 445869 56. Dix, J., Lue, L., Carbone, P.: J. Comput. Chem. 39(25), 2051 (2018). https:// doi. org/ 10. 1002/ jcc. 25369 57. Nosé, S.: Mol. Phys. 52(2), 255 (1984). https:// doi. org/ 10. 1080/ 00268 97840 01012 01 58. Hoover, W.G.: Phys. Rev. A 31(3), 1695 (1985). https:// doi. org/ 10. 1103/ physr eva. 31. 1695 59. Parrinello, M., Rahman, A.: Phys. Rev. Lett. 45(14), 1196 (1980). https:// doi. org/ 10. 1103/ physr evlett. 45. 1196 60. Parrinello, M., Rahman, A.: J. Appl. Phys. 52(12), 7182 (1981). https:// doi. org/ 10. 1063/1. 328693 61. Hess, B., Bekker, H., Berendsen, H.J.C., Fraaije, J.G.E.M.: J. Comput. Chem. 18(12), 1463 (1997). https:// doi. org/ 10. 1002/ (sici) 1096- 987X(199709) 18: 12% 3C3c1 463:: aid- jcc4% 3C3e3.0. co;2-h 62. Hess, B.: J. Chem. Theory Comput. 4(1), 116 (2008). https:// doi. org/ 10. 1021/ ct700 200b 63. Darden, T., York, D., Pedersen, L.: J. Chem. Phys. 98(12), 10089 (1993). https:// doi. org/ 10. 1063/1. 464397 64. Essmann, U., Perera, L., Berkowitz, M.L., Darden, T., Lee, H., Pedersen, L.G.: J. Chem. Phys. 103(19), 8577 (1995). https:// doi. org/ 10. 1063/1. 470117 65. Kirkwood, J.G.: J. Chem. Phys. 3(5), 300 (1935). https:// doi. org/ 10. 1063/1. 17496 57 66. Kraus, H., Hansen, N.: Poreana: 0.2.0 (2021). https:// doi. org/ 10. 5281/ zenodo. 47381 82 67. Langmuir, I.: J. Am. Chem. Soc. 40(9), 1361 (1918). https:// doi. org/ 10. 1021/ ja022 42a004 68. Duarte Ramos Matos, G., Kyu, D.Y., Loeffler, H.H., Chodera, J.D., Shirts, M.R., Mobley, D.L.: J. Chem. Eng. Data 62(5), 1559 (2017). https:// doi. org/ 10. 1021/ acs. jced. 7b001 04 69. Dortmund Data Bank. www. ddbst. com (2014) 70. Cabani, S., Gianni, P., Mollica, V., Lepori, L.: J. Solut. Chem. 10(8), 563 (1981). https:// doi. org/ 10. 1007/ BF006 46936 71. Al-Burtomani, S.K.S., Suliman, F.O.: RSC Adv. 7(16), 9888 (2017). https:// doi. org/ 10. 1039/ C6RA2 8638K 72. Tang, Z., Chang, C.A.: J. Chem. Theory Comput. 14(1), 303 (2018). https:// doi. org/ 10. 1021/ acs. jctc. 7b008 99 73. Rekharsky, M.V., Inoue, Y.: Chem. Rev. 98(5), 1875 (1998). https:// doi. org/ 10. 1021/ cr970 015o 74. Cramer, F., Saenger, W., Spatz, H.C.: J. Am. Chem. Soc. 89(1), 14 (1967). https:// doi. org/ 10. 1021/ ja009 77a003 75. Chen, W., Chang, C.E., Gilson, M.K.: Biophys. J. 87(5), 3035 (2004). https:// doi. org/ 10. 1529/ bioph ysj. 104. 049494 76. Jani, A., Busch, M., Mietner, J.B., Ollivier, J., Appel, M., Frick, B., Zanotti, J.M., Ghoufi, A., Huber, P., Fröba, M., Morineau, D.: J. Chem. Phys. 154(9), 094505 (2021). https:// doi. org/ 10. 1063/5. 00407 05 77. Schwarz, F.P.: Anal. Chem. 52(1), 10 (1980). https:// doi. org/ 10. 1021/ ac500 51a004 78. Achard, C., Jaoui, M., Schwing, M., Rogalski, M.: J. Chem. Eng. Data 41(3), 504 (1996). https:// doi. org/ 10. 1021/ je950 202o 79. Shen, H.M., Zhu, G.Y., Yu, W.B., Wu, H.K., Ji, H.B., Shi, H.X., Zheng, Y.F., She, Y.B.: RSC Adv. 5, 84410 (2015). https:// doi. org/ 10. 1039/ C5RA1 5592D 80. Karger, B.L., Snyder, L.R., Horvath, C.: An Introduction to Sepa- ration Science. Wiley, New York (1973) Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. https://doi.org/10.1080/08927022.2020.1871478 https://doi.org/10.1080/08927022.2020.1871478 https://doi.org/10.5281/zenodo.4738036 https://doi.org/10.5281/zenodo.4738036 https://doi.org/10.1021/la062285k https://doi.org/10.5281/zenodo.4588861 https://doi.org/10.5281/zenodo.4588861 https://doi.org/10.1016/j.softx.2015.06.001 https://doi.org/10.1016/j.softx.2015.06.001 https://doi.org/10.1021/ct700301q https://doi.org/10.1021/ct700301q https://doi.org/10.1016/j.cpc.2009.05.011 https://doi.org/10.1016/j.cpc.2013.09.018 https://doi.org/10.1016/j.cpc.2013.09.018 https://doi.org/10.1021/jacs.9b08776 https://doi.org/10.1021/jacs.9b08776 https://doi.org/10.1021/acscatal.1c03057 https://doi.org/10.1002/jcc.20035 https://doi.org/10.1002/jcc.20035 https://doi.org/10.1021/acs.jpcb.7b11808 https://doi.org/10.1021/acs.jpcb.7b11808 https://doi.org/10.5281/zenodo.1161245 https://doi.org/10.5281/zenodo.1161245 https://doi.org/10.1063/1.445869 https://doi.org/10.1063/1.445869 https://doi.org/10.1002/jcc.25369 https://doi.org/10.1080/00268978400101201 https://doi.org/10.1080/00268978400101201 https://doi.org/10.1103/physreva.31.1695 https://doi.org/10.1103/physreva.31.1695 https://doi.org/10.1103/physrevlett.45.1196 https://doi.org/10.1063/1.328693 https://doi.org/10.1002/(sici)1096-987X(199709)18:12%3C3c1463::aid-jcc4%3C3e3.0.co;2-h https://doi.org/10.1002/(sici)1096-987X(199709)18:12%3C3c1463::aid-jcc4%3C3e3.0.co;2-h https://doi.org/10.1002/(sici)1096-987X(199709)18:12%3C3c1463::aid-jcc4%3C3e3.0.co;2-h https://doi.org/10.1021/ct700200b https://doi.org/10.1021/ct700200b https://doi.org/10.1063/1.464397 https://doi.org/10.1063/1.470117 https://doi.org/10.1063/1.470117 https://doi.org/10.1063/1.1749657 https://doi.org/10.1063/1.1749657 https://doi.org/10.5281/zenodo.4738182 https://doi.org/10.5281/zenodo.4738182 https://doi.org/10.1021/ja02242a004 https://doi.org/10.1021/ja02242a004 https://doi.org/10.1021/acs.jced.7b00104 http://www.ddbst.com https://doi.org/10.1007/BF00646936 https://doi.org/10.1039/C6RA28638K https://doi.org/10.1021/acs.jctc.7b00899 https://doi.org/10.1021/cr970015o https://doi.org/10.1021/ja00977a003 https://doi.org/10.1529/biophysj.104.049494 https://doi.org/10.1063/5.0040705 https://doi.org/10.1063/5.0040705 https://doi.org/10.1021/ac50051a004 https://doi.org/10.1021/ac50051a004 https://doi.org/10.1021/je950202o https://doi.org/10.1039/C5RA15592D https://doi.org/10.1039/C5RA15592D An atomistic view on the uptake of aromatic compounds by cyclodextrin immobilized on mesoporous silica Abstract 1 Introduction 2 Methodology 2.1 Calculation of binding free enthalpies and rate constants in bulk solution 2.2 Immobilization of cyclodextrins 2.3 Simulated systems 2.4 Simulation parameters 2.5 Analysis 2.6 Langmuir model 3 Results and discussion 3.1 Bulk phase simulations 3.2 Unfunctionalized silica pores 3.3 Cyclodextrin-functionalized silica pores 4 Conclusion and outlook Acknowledgements References