An atomistic view on the uptake of aromatic compounds by cyclodextrin immobilized on mesoporous silica

The effect of immobilized β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upbeta$$\end{document}-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.


Introduction
Cyclodextrins (CDs) are cyclic oligosaccharides with conical 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 chiral selectors [3] or in (asymmetric) supramolecular catalysis [4][5][6][7][8][9][10] among many other fields of applications including the use of cyclodextrins as drug delivery agents or as building 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 scaffold allows to introduce additional catalytic groups or binding 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 β-cyclodextrinsilane in the presence of a structure directing agent [17][18][19][20], resulting in silica-based materials that possess a uniform framework mesoporosity with defined nanoscaled cavities. The ability of removing aromatic compounds from an aqueous 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 fundamental understanding at the molecular level. Gas adsorption in microporous materials such as zeolites [21,22] or other nanostructured solids [23] has been studied for more than 45 years in particular using classical molecular simulations. The combination of molecular simulation techniques with experimental measurements allows to examine in detail fundamental diffusion 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 (multicomponent) mixtures at a solid surface functionalized for example with alkyl chains [25][26][27][28][29]. Other studies used coarsegrained 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 aqueous 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 isotherms is discussed in view of the underlying molecular level picture.

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 counting (DC), the occurrences of bound, N b , and unbound, N u , 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 V 0 = 1.661 nm 3 and average volume of the simulation box V box . In addition, average bound ⟨t b ⟩ and unbound ⟨t u ⟩ residence times can be calculated, yielding association k on and dissociation k of f rates [34] (1) with guest molecule concentration C g . The binding free enthalpy can then be determined using these rate constants (RC) with standard state concentration C 0 = 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 circle 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 different duration. Averaging over these instances yields time averages ⟨t b ⟩ and ⟨t u ⟩ . 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 G M→M � u is calculated, resulting from decoupling the unbound state u, i.e., turning off the 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 G M 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 G M→M � tor and lastly turning off the initial restraints G M � 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 G M→M � b→tor = G M b→tor + G M→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 summation over the whole cycle is equal to zero thus resulting into (4) 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 r aA,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 K i . 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]

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 structures of the linkers are illustrated in Fig. 3.

Simulated systems
For calculating the binding free enthalpy by direct counting and via rate constants, long unbiased simulations were conducted in the NpT ensemble with one guest molecule, benzene or p-nitrophenol, respectively, and one host molecule, β-cyclodextrin, solvated in 1000 water molecules. Simulations 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 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 reservoir 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 experimentally 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 oxygen atoms and silanol groups are taken from Gulmen and Thompson [44] and are shown in the Supplementary Material 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, 298 K and 350 K and two cyclodextrin surface densities. For p-nitrophenol an adsorption isotherm was calculated at 350 K using the pore with the higher cyclodextrin surface density.

Fig. 3
Atom names and partial atomic charges (purple) of cyclodextrin (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 Trofymchuk et al. [20], referred to as L2 in the present work, with a total charge of − 0.32e (Color figure online) 6

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 specific 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 topologies 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 Lennard-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 intermolecular interactions of the guest molecule with its environment 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 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.

Analysis
Distances and angles between the reference systems during 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 during 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 r i of slice i. The number density i is then determined by counting the number of molecules N i within each slice during the simulation, resulting into Adsorption isotherms describe the amount of solute molecules adsorbed on the surface N pore as a function of the amount of solute in the bulk phase N bulk and are therefore, similarly to the density, determined from counting the number of molecules within the pore and within the bulk reservoirs normalized by the number of frames with the number of molecules N j within the pore or bulk phase at frame j = 1, … , N F . These values are then converted 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 displacement i (t) over an observation time of 4-20 ps within each bin with a tolerance of ±1 bins 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.

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 simulation 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 q max denotes the cyclodextrin density inside the pore and (10) (16) q ads = q max K ⋅ C 1 + K ⋅ C with C 0 as the concentration of the standard state, i.e. 1 mol 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 simulation 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.

Bulk phase simulations
In order to verify the quality of the force field for the solute molecules p-nitrophenol and benzene, validating simulations 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 G hyd 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 constants 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 gyration around 0.6 nm [71], the spatial cutoff distance was chosen 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 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 depending on the temporal resolution of the measurement device. For benzene the effect on G is minor, while for p-nitrophenol an effect is only visible between temporal cut-offs of 0 ps and 100 ps. For both molecules the k on -and k of f -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 298 K 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 Material in Fig. S1. Therefore only the double decoupling results are reported at this temperature. At 350 K unbound occurrences are more likely due to the higher temperature, allowing 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 indicates that longer unbound instances are still infrequent. The reason for the larger binding free enthalpy at a higher temperature for benzene in the unbiased approaches is due to a Table 2 Liquid densities ( kg m −3 ) for p-nitrophenol (p-NP) and benzene (BEN) compared with experimental data [69] Determined at a 387 K , b 298 K . Hydration free enthalpy G hyd ( kJ mol −1 ) at 298 K 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]   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 distributions are in good agreement with those reported by Gebhardt et al. [53] using the same force field. Additional dihedralangle 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). 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 200 mmol 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 300 K the experimental self-diffusion coefficient of water in the bulk phase, 2.3 × 10 −9 m 2 s −1 , is reduced to 2.0 × 10 −9 m 2 s −1 in SBA-15 with a pore mean pore diameter of 6.6 nm and to 1.7 × 10 −9 m 2 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 m 2 s −1 at a temperature of 298 K and 5.63 × 10 −9 m 2 s −1 at a higher temperature of 350 K . This yields a factor 1.70 at 298 K and 1.74 at 350 K by comparing the pore diffusion to the bulk diffusion of the same system.

Cyclodextrin-functionalized silica pores
By functionalizing the surface, benzene molecules have a significantly higher concentration at the configurational 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 298 K and a high cyclodextrin concentration at 298 K and 350 K . For p-nitrophenol only the larger surface concentration at 350 K 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 ( kJ mol −1 ) with experimental (Exp) data of β-cyclodextrin-ligand complexes, with association k on ( 10 8 dm 3 mol −1 s −1 ) and dissociation k of f ( 10 6 s −1 ) rates at a temporal cut-off of 1 ns calculated at temperature T ( K) Mean value from a DC, DD, and b DC, 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) 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 cyclodextrin 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 7 mmol l −1 , i.e., roughly one third of the solubility limit of 23.8 mmol 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 60 mmol l −1 were used to improve the statistical sampling. However, the amount adsorbed did not exceed the 1:1 binding capacity 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. 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 association 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 capacities significantly larger than those corresponding to a 1:1 binding scenario and explained this finding with hydrogen 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 concentration 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.

Conclusion and outlook
An efficient model building is a prerequisite for computational studies of functionalized mesoporous silica materials. The Python package PoreMS introduced previously [42,43] was complemented by two additional program packages for preparing MD simulations of porous materials with GROMACS, PoreSim [45], and for analyzing the simulation trajectories, PoreAna [66], the latter providing results such as density and diffusion profiles, thereby reducing the overhead for system preparation to analysis. The selected case study of adsorption of aromatic molecules 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 adsorption isotherms show a more complex behavior than predicted 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 interactions on the surface as well as more unspecific host-guest interactions have an influence on the shape of the isotherm. 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 investigate 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 298 K (dashed lines) and 350 K (solid lines). The grey line represents a hypothetical maximum of 1:1 binding with 11 cyclodextrin molecules (Color figure online)