1 Introduction

Cyclodextrins (CDs) are cyclic oligosaccharides with conical shape built of (1-4)-linked \(\upalpha\)-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 \(\upbeta\)-cyclodextrin-silane 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 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 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.

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 counting (DC), the occurrences of bound, \(N_\mathrm {b}\), and unbound, \(N_\mathrm {u}\), instances during long (\(t\ge 10\,\upmu \hbox {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]

$$\begin{aligned} \varDelta G^\text {DC}=-RT\ln \frac{N_\mathrm {b}}{N_\mathrm {u}}-RT\ln \frac{V_\mathrm {box}}{V^0} \end{aligned}$$
(1)

with standard state volume \(V^0=1.661 \hbox { nm}^{3}\) and average volume of the simulation box \(V_\mathrm {box}\). In addition, average bound \(\langle t_b\rangle\) and unbound \(\langle t_u\rangle\) residence times can be calculated, yielding association \(k_\mathrm {on}\) and dissociation \(k_\mathrm {off}\) rates [34]

$$\begin{aligned} k_\mathrm {on}=\frac{1}{\langle t_\mathrm {u}\rangle C_\mathrm {g}},\ \ k_\mathrm {off}=\frac{1}{\langle t_\mathrm {b}\rangle } \end{aligned}$$
(2)

with guest molecule concentration \(C_\mathrm {g}\). The binding free enthalpy can then be determined using these rate constants (RC)

$$\begin{aligned} \begin{aligned} \varDelta G^\text {RC}&=-RT\ln \frac{k_\mathrm {on}C^0}{k_\mathrm {off}}\\&=-RT\ln \frac{\langle t_\mathrm {b}\rangle }{\langle t_\mathrm {u}\rangle }-RT\ln \frac{V_\text {box}}{V^0} \end{aligned} \end{aligned}$$
(3)

with standard state concentration \(C^0=\hbox {mol }\hbox{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.

Fig. 1
figure 1

Geometrical representation of the cyclodextrin and guest molecule structures by a three point system

Full size image

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 \(\langle t_\mathrm {b}\rangle\) and \(\langle t_\mathrm {u}\rangle\). 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 \(\varDelta G_\mathrm {u}^{\mathrm {M}\rightarrow {\mathrm {M}^{\prime }}}\) is calculated, resulting from decoupling the unbound state u, i.e., turning off the intermolecular interactions with the environment (\(\mathrm {M}\rightarrow {\mathrm {M}^{\prime }}\)) of the guest molecule in a box of water while preserving intramolecular interactions. \(\varDelta G_\mathrm {u}^{\mathrm {M}\rightarrow {\mathrm {M}^{\prime }}}\) is equal to the negative hydration free enthalpy \(-\varDelta G_\mathrm {hyd}^\mathrm {M}\). Second, the free enthalpy difference \(\varDelta G_\mathrm {b}^{\mathrm {M}\rightarrow {\mathrm {M}^{\prime }}}\) 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 \(\varDelta G_{\mathrm {b}\rightarrow \mathrm {tor}}^\mathrm {M}\) 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 \(\varDelta G_{\mathrm {tor}}^{\mathrm {M}\rightarrow {\mathrm {M}^{\prime }}}\) and lastly turning off the initial restraints \(\varDelta G_{\mathrm {tor}}^{\mathrm {M}^{\prime }}\). 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 \(\varDelta G_{\mathrm {b}\rightarrow \mathrm {tor}}^{\mathrm {M}\rightarrow {\mathrm {M}^{\prime }}}=\varDelta G_{\mathrm {b}\rightarrow \mathrm {tor}}^\mathrm {M}+\varDelta G_{\mathrm {tor}}^{\mathrm {M}\rightarrow {\mathrm {M}^{\prime }}}\). 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

$$\begin{aligned} \begin{aligned}\varDelta G^\text {DD}&=\varDelta G_{\mathrm {u\rightarrow b}}^\mathrm {M}\\&= -\varDelta G_\mathrm {hyd}^\mathrm {M}-\varDelta G_{\mathrm {b\rightarrow tor}}^\mathrm {M}-\varDelta G_{\mathrm {tor}}^{\mathrm {M}\rightarrow {\mathrm {M}^{\prime }}}-\varDelta G_{\mathrm {tor}}^{\mathrm {M}^{\prime }}\\& =-\varDelta G_\mathrm {hyd}^\mathrm {M}-\varDelta G_{\mathrm {b\rightarrow tor}}^{\mathrm {M}\rightarrow {\mathrm {M}^{\prime }}}-\varDelta G_{\mathrm {tor}}^{\mathrm {M}^{\prime }}. \end{aligned} \end{aligned}$$
(4)
Fig. 2
figure 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) interaction (full triangle), while \(\text{M}^{\prime }\) 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

Full size image

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 \(\theta _{A,0}\), \(\theta _{B,0}\), and three dihedral angles \(\phi _{A,0}\), \(\phi _{B,0}\), \(\phi _{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

$$\begin{aligned} \begin{aligned}&\varDelta G_{\mathrm {tor}}^{\mathrm {M}^{\prime }}\\&\quad = -RT\ln \left[ \frac{8\pi V^0(K_rK_{\theta _A}K_{\theta _B}K_{\phi _A}K_{\phi _B}K_{\phi _C})^\frac{1}{2}}{r_{aA,0}^2\sin \theta _{A,0}\sin \theta _{B,0}(2\pi RT)^3}\right] \end{aligned} \end{aligned}$$
(5)

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]

$$\begin{aligned} \varDelta {\bar{G}}_{\mathrm {u\rightarrow b}}^\mathrm {M}=-RT\ln \left[ \sum _k\exp \left( \frac{\varDelta G_{\mathrm {u\rightarrow b},k}^\mathrm {M}}{RT}\right) \right] . \end{aligned}$$
(6)

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 AmberTools20 [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.

Fig. 3
figure 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 \(\hbox {Na}^+\)) of \(-\,1.32\hbox{e}\) and by b Trofymchuk et al. [20], referred to as L2 in the present work, with a total charge of \(-\,0.32\hbox{e}\) (Color figure online)

Full size image

2.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, \(\upbeta\)-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\hbox { K}\) and \(350\hbox { K}\) at a constant pressure of \(1\hbox { 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) \(\upbeta\)-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\,\upmu \hbox {mol }\hbox {m}^{{-2}}\) \(\upbeta\)-cyclodextrin, respectively. For the simulations representing the system by Trofymchuk et al. [20], additional 0.11 and \(0.22\,\upmu \hbox {mol }\hbox {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\,\upmu \hbox {mol }\hbox {m}^{{-2}}\) on the internal surface and \(8.82\,\upmu \hbox {mol }\hbox {m}^{{-2}}\) on the external surface. Further properties of the pore are listed Table 1.

Fig. 4
figure 4

Pore system functionalized with \(\upbeta\)-cyclodextrin using the surface linker used by Trofymchuk et al. [20]. a Side view of the simulation box indicating the length of the central silica block and the solvent reservoirs. b Front view of the pierced silica block containing the 4 nm diameter pore. The chemistry of the exterior surface is based on the (111) face of \(\upbeta\)-cristobalite silica. Bonded-phase groups are randomly distributed on the silica surface. Ligand densities, residual surface hydroxylation, and further details are specified in Table 1. Colour code: Si, yellow lines; O, red lines; \(\upbeta\)-cyclodextrin, blue; propylamine groups, magenta; surface silanol groups, yellow (Color figure online)

Full size image
Table 1 Properties of the cylindrical mesopore model before and after functionalization (func.), generated as described by Trofymchuk et al. [20], with surface densities (\(\upmu \hbox {mol }\hbox {m}^{{-2}}\)) and pore dimensions (nm)
Full size table

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\,\hbox {K}\) and \(350\,\hbox {K}\) and two cyclodextrin surface densities. For p-nitrophenol an adsorption isotherm was calculated at \(350\,\hbox {K}\) using the pore with the higher cyclodextrin surface density.

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 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 \times 10^{-5}\,\hbox {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\,\upmu \hbox {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 (\(\lambda\)-points) for a slow equidistant deactivation of the intermolecular interactions of the guest molecule with its environment while preserving intramolecular interactions. Each \(\lambda\)-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. For these systems a trajectory length of \(1\,\upmu \hbox {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 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

$$\begin{aligned} V_i=\pi z(r_i^2-r_{i-1}^2) \end{aligned}$$
(7)

with pore length z and radius \(r_i\) of slice i. The number density \(\rho _i\) is then determined by counting the number of molecules \(N_i\) within each slice during the simulation, resulting into

$$\begin{aligned} \rho _i=\frac{N_i}{V_i}=\frac{N_i}{\pi z}\frac{1}{r_i^2-r_{i-1}^2}. \end{aligned}$$
(8)

Adsorption isotherms describe the amount of solute molecules adsorbed on the surface \(N^\mathrm {pore}\) as a function of the amount of solute in the bulk phase \(N^\mathrm {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

$$\begin{aligned} \begin{aligned}&\langle N\rangle ^\mathrm {pore}=\frac{1}{N_\mathrm {F}}\sum _{j=1}^{N_\mathrm {F}}N_j^\mathrm {pore},\\&\langle N\rangle ^\mathrm {bulk}=\frac{1}{N_\mathrm {F}}\sum _{j=1}^{N_\mathrm {F}}N_j^\mathrm {bulk}. \end{aligned} \end{aligned}$$
(9)

with the number of molecules \(N_j\) within the pore or bulk phase at frame \(j=1,\dots ,N_\mathrm {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_{\Vert ,i}\) parallel to the pore surface was determined from the slope of the mean square displacement \(\varDelta _i(t)\) over an observation time of 4-20 ps within each bin with a tolerance of \(\pm 1\) bins

$$\begin{aligned} D_{\Vert ,i}=\frac{1}{2}\frac{d\varDelta _i(t)}{dt}. \end{aligned}$$
(10)

By weighting the axial diffusion profile \(D_{\Vert ,i}\) with the density profile \(\rho _i\) along the radius r, a mean diffusion coefficient \(\langle D_{\Vert }\rangle\) can be calculated

$$\begin{aligned} \langle D_{\Vert }\rangle =\frac{\sum \rho _iD_{\Vert ,i}A_i}{\sum \rho _iA_i} =\frac{\sum \rho _iD_{\Vert ,i}(r_i^2-r_{i-1}^2)}{\sum \rho _i(r_i^2-r_{i-1}^2)} \end{aligned}$$
(11)

with the cross-sectional area

$$\begin{aligned} A_i=\pi (r_i^2-r_{i-1}^2) \end{aligned}$$
(12)

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 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

$$\begin{aligned} \left( \frac{N_\mathrm {b}}{N_\mathrm {u}}\right) _\mathrm {pore} = \left( \frac{N_\mathrm {b}}{N_\mathrm {u}}\right) _\mathrm {bulk}\cdot \frac{V_\mathrm {box}^\mathrm {bulk}}{V_\mathrm {box}^\mathrm {pore}} \end{aligned}$$
(13)

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

$$\begin{aligned} N_\mathrm {b}^\mathrm {pore}=\left( \frac{N_\mathrm {u}}{V_\mathrm {box}}\right) _\mathrm {pore}\cdot \left( \frac{N_\mathrm {b}}{N_\mathrm {u}}\right) _\mathrm {bulk}\cdot V_\mathrm {box}^\mathrm {bulk} \end{aligned}$$
(14)

and replacing the two rightmost terms by means of Eq. 1 results in

$$\begin{aligned} N_\mathrm {b}^\mathrm {pore}=\left( \frac{N_\mathrm {u}}{V_\mathrm {box}}\right) _\mathrm {pore}\cdot V^0\cdot \exp \left( -\frac{\varDelta G^\mathrm {DC}}{RT}\right) . \end{aligned}$$
(15)

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]

$$\begin{aligned} q_\mathrm {ads}=q_\mathrm {max}\frac{K\cdot C}{1 + K\cdot C} \end{aligned}$$
(16)

where \(q_\mathrm {max}\) denotes the cyclodextrin density inside the pore and

$$\begin{aligned} K= \frac{1}{C^0}\exp \left( -\frac{\varDelta G}{RT}\right) \end{aligned}$$
(17)

with \(C^0\) as the concentration of the standard state, i.e. \(1\,\mathrm {mol\,l^{-1}}\) and \(\varDelta 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.

3 Results and discussion

3.1 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 \(\varDelta G_\mathrm {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 2 Liquid densities \(\rho\) (\(\hbox {kg}\,\hbox{m}^{-3}\)) for p-nitrophenol (p-NP) and benzene (BEN) compared with experimental data [69]
Full size table

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 \(\upbeta\)-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.

Fig. 5
figure 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=350\,\hbox {K}\) with the corresponding image showing benzene in the bound state (Color figure online)

Full size image

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 \(\varDelta 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_{\mathrm {on}}\)- and \(k_{\mathrm {off}}\)-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].

Table 3 Comparison of direct counting (DC), rate constants (RC) and double decoupling (DD) methods for determining the binding free enthalpy \(\varDelta G\) (\(\hbox {kJ}\,\mathrm {mol}^{-1}\)) with experimental (Exp) data of \(\upbeta\)-cyclodextrin-ligand complexes, with association \(k_\mathrm {on}\) (\(10^{8}\hbox { dm}^{3}\hbox { mol}^{-1}\hbox { s}^{-1}\)) and dissociation \(k_\mathrm {off}\) (\(10^{6}\hbox { s}^{-1}\)) rates at a temporal cut-off of 1 ns calculated at temperature T (\(\hbox {K}\))
Full size table

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\,\hbox {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\,\hbox {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 larger volumetric correction in Eq. ( 1) at \(RT_{298}\ln \frac{V_\text {box}}{V^0}=7.28\hbox { kJ mol}^{-1}\) and, \(RT_{350}\ln \frac{V_\text {box}}{V^0}=8.69\hbox { kJ mol}^{-1}\), which is not entirely compensated by the decrease of the ratio of bound and unbound instances \(RT_{298}\ln \frac{N_\text {b}}{N_\text {u}}=5.98\hbox { kJ mol}^{-1}\) and, \(RT_{350}\ln \frac{N_\text {b}}{N_\text {u}}=5.41\hbox { kJ mol}^{-1}\), as they stay almost similar.

Fig. 6
figure 6

a Distributions of the seven individual dihedral angles H6O–O6–C6–H61 for the native \(\upbeta\)-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 \(\upbeta\)-cyclodextrin connected with the L1 linker [17]. The blue line represents the distribution of the C5(tail)–O6–C6–H61 dihedral 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 dihedral angle (Color figure online)

Full size image

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 dihedral-angle distributions for \(\upalpha\)- and \(\upgamma\)-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^{\circ }\)  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 \(200\,\hbox {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\,\hbox {K}\) the experimental self-diffusion coefficient of water in the bulk phase, \(2.3 \times 10^{-9}\hbox { m}^{2}\hbox { s}^{-1}\), is reduced to \(2.0 \times 10^{-9}\hbox { m}^{2}\hbox { s}^{-1}\) in SBA-15 with a pore mean pore diameter of 6.6 nm and to \(1.7\times 10^{-9}\hbox { m}^{2}\hbox { 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 \(\langle D_\Vert \rangle\) of water in the simulated unfunctionalized pore is \(3.07\times 10^{-9}\hbox { m}^{2}\hbox { s}^{-1}\) at a temperature of \(298\,\hbox {K}\) and \(5.63\times 10^{-9}\hbox { m}^{2}\hbox { s}^{-1}\) at a higher temperature of \(350\,\hbox {K}\). This yields a factor 1.70 at \(298\,\hbox {K}\) and 1.74 at 350 K by comparing the pore diffusion to the bulk diffusion of the same system.

Fig. 7
figure 7

a Density and b axial diffusion profiles of water (blue), benzene (BEN) (orange), and p-nitrophenol (p-NP) (green) in a non-functionalized pore at temperatures \(298\,\hbox {K}\) (dashed lines) and \(350\,\hbox {K}\) (solid lines). The shaded area denotes the configurational space of the silanol group oxygen atoms (yellow) (Color figure online)

Full size image

3.3 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.

Fig. 8
figure 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\,\upmu \hbox {s}\). Shaded areas denote the configurational space of the cyclodextrin centre of mass (blue) and the silanol groups oxygen atoms (yellow) (Color figure online)

Full size image
Fig. 9
figure 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-variant cyclodextrin [17] with p-nitrophenol and the mean free enthalpy value only from the DD and DC approach. Simulation were run for \(1\,\upmu \hbox {s}\) at \(298\,\hbox {K}\) (dashed lines) and \(350\,\hbox {K}\) (solid lines). The grey line represents a hypothetical maximum of 1:1 binding with 11 cyclodextrin molecules (Color figure online)

Full size image

Benzene adsorption isotherms were calculated for a low cyclodextrin surface concentration at \(298\,\hbox {K}\) and a high cyclodextrin concentration at \(298\,\hbox {K}\) and \(350\,\hbox {K}\). For p-nitrophenol only the larger surface concentration at \(350\,\hbox {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 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\,\hbox {mmol l}^{-1}\), i.e., roughly one third of the solubility limit of \(23.8\,\hbox {mmol}\,\mathrm {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\,\hbox {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.

4 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 \(\upmu \hbox {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.