Article pubs.acs.org/JPCB

Thermodynamic Driving Forces for Dye Molecule Position and Orientation in Nanoconfined Solvents Jacob A. Harvey and Ward H. Thompson* Department of Chemistry, University of Kansas, Lawrence, Kansas 66045, United States S Supporting Information *

ABSTRACT: The results of replica exchange molecular dynamics simulations of a coumarin 153 (C153) dye molecule dissolved in ethanol confined within a 2.4 nm hydrophilic amorphous silica pore are presented. The C153 dye position and orientation distributions provide insight into time-dependent fluorescence measurements in nanoconfined solvents as well as general features of chemistry in mesoporous materials. In addition to the distributions themselves, the free energy, internal energy, and entropic contributions have been calculated to explore the factors determining the distributions. The most likely location of C153 is found to be near the pore surface, but two possible hydrogen-bonding structures lead to differing orientations. Internal energy and entropy are found to be competing forces within the pore, with entropy playing a significant role with unexpected consequences. These results represent a crucial step in determining how the nanoconfining framework can affect measurements of solvation dynamics.

■

INTRODUCTION Nanoconfined liquids exist in a variety of materials including membranes, mesoporous silica, reverse micelles, and proteins. Understanding how the various confining framework properties, including the tunable surface chemistries, affect the properties of the confined liquid is of paramount importance in applications such as catalysis, drug delivery, and chemical sensors. Effects of nanoconfinement on liquids include changes in structure,1−8 diffusion coefficients,9−12 and, of particular relevance to this work, solvation dynamics.5,13−29 Time-dependent fluorescence (TDF) measurements are the primary experimental tool for probing solvation dynamics in nanoconfined solvents. Virtually all TDF measurements show significantly slower dynamics, sometimes appearing as new time scales not present in the bulk solvent. A number of theoretical models have been proposed to explain these results,30,31 including our suggestion that fluorescent probe solutes can diffuse when electronically excited from their ground state and that the diffusion may be observable as a time scale in the TDF signal.24,32−34 Thus, understanding the distribution of positions of a dye molecule, and how it is affected by the confining environment, is important for testing this hypothesis. More generally, this argument applies to all cases involving two species with different charge distributions, e.g., ground- and excited-state molecule, reactant or product, and thus it has broader implications for chemistry in mesoporous materials. Here we investigate the position distribution of a common fluorescent probe, coumarin 153 (C153), inside a nanoscale silica pore. In doing so we also examine the temperaturedependent free energies to elucidate the thermodynamic driving forces, particularly the roles of entropy and hydrogen bonding. © XXXX American Chemical Society

Experiments have been used extensively to study solvation and reorientation dynamics of nanoconfined liquids, and covering the large body of work on the subject is beyond the scope of this work. We therefore refer the reader to recent reviews and the references therein.33,35,36 However, one study in particular has motivated our work on the system considered here. Baumann et al.21 used TDF measurements of C153 to study the solvation dynamics of ethanol in 2.5 and 5.0 nm sol− gel glasses. Three solvation time scales were found for both the bulk and confined liquids; however the longest time scale increased from 26.3 ps in the bulk to 50.7 and 102.7 ps in the 5.0 and 2.5 nm pores, respectively. They suggested that there is a preferential alignment of polar solvent molecules at the pore interface and that the change in the polarization field would likely be the cause of the slower solvation dynamics.21 Fully atomistic simulations of large fluorescent probes in neat liquids are now relatively routine, but few studies involving nanoconfined systems have been reported.37−40 Those studies have generally involved 20−50 ns equilibrium simulations and/ or 30−200 ps for nonequilibrium dynamics simulations. However, our preliminary simulations suggested that a dye molecule such as C153 experiences significant barriers in a nanoconfined liquid, making statistical sampling a challenging issue. Thus, accelerated approaches are required.40 Therefore, in this work we use replica exchange molecular dynamics to determine the position and orientation distributions of groundstate C153 dissolved in ethanol in hydrophilic silica pores. This Special Issue: Branka M. Ladanyi Festschrift Received: September 7, 2014 Revised: October 8, 2014

A

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Simulation Details. Replica exchange molecular dynamics (REMD)51−53 simulations were carried out as implemented in the large-scale atomic/molecular massively parallel simulator (LAMMPS) molecular dynamics package.54 In this approach replicas of the system are simulated at different temperatures simultaneously. Attempts to exchange coordinates from adjacent (i.e., j = i ± 1) temperatures are made at regular intervals and allowed to exchange with probability

approach provides the additional advantage that the contributions of entropy and internal energy to the overall free energy curve can be determined.

■

METHODS Pore Model. All simulations have been performed with silica pore models of diameter approximately 2.4 nm. A brief description of the system will be repeated here; however a more detailed discussion can be found in our previous work.41−43 The pore framework consists of 1042 Si atoms and 2144 O atoms which remained fixed during the simulation. The pore surface is functionalized with 36 silanol, SiOH, groups and 6 geminal, Si(OH)2, groups. Surface O−H bond lengths were held constant using the SHAKE algorithm,44 but the OH groups are allowed to rotate with Si−O−H angles described by an harmonic potential. All pore force field parameters can be found in the Supporting Information. The pore was filled with ethanol using a grand canonical Monte Carlo simulation. A C153 molecule was then inserted at random and any overlapping ethanol molecules were removed, which yielded a final system consisting of 121 ethanol molecules and 1 C153 molecule that was used for all of the simulations described later. A total of 3844 atoms were contained in a 44 × 44 × 30 Å3 simulation cell. A snapshot of the pore structure is shown in Figure 1.

A = min{1, exp[+ (βi − βj)(E(riN ) − E(r jN ))]}

(1)

where is the total energy at inverse temperature βi = 1/ kBTi. Configurations at higher temperatures are able to more freely sample the free energy landscape while at lower temperatures the C153 molecule will likely be constrained to a local free energy well. Thus, if each configuration completes enough round trips from the lowest temperature to the highest temperature all free energy wells should be sampled. The choice of temperatures is nontrivial and deserves some discussion here. The lowest temperature was chosen to be the temperature of interest, 298 K, while the highest temperature ideally allows for uninhibited sampling of the free energy surface. It is difficult to know, a priori, how large the barriers in the system might be. To explore this, we ran 10 ns simulations at 598, 648, 698, 748, and 798 K. While the C153 molecule can still reside at the surface for extended times, e.g., ∼2 ns, at the higher temperatures, these periods are significantly shorter and less frequent than at 298 K. Based on these results and the determination of intermediate temperatures (see below), a maximum temperature of 688 K was chosen as a reasonable trade-off between required computing time and the sampling ability of the maximum temperature. Intermediate temperatures were chosen according to the procedure of Rathore et al.,55 where it was shown that, assuming Gaussian distributions of energy, the acceptance probability between temperatures is dependent on the ratio E(rNi )

⟨E(Ti )⟩ − ⟨E(Tj)⟩ (σ(Ti ) + σ(Tj))/2

(2)

where ⟨E(Ti)⟩ and σ(Ti) are the average and standard deviation, respectively, of the total energy at temperature Ti. An analytical function for ⟨E(T)⟩ and σ(T) was determined by shorter simulations at a few temperatures with both quantities found to vary linearly with T. The Ti values were then determined by iteratively solving

Figure 1. Structure of the pore and a C153 molecule (orange) dissolved in ethanol (magenta). The four atomistic sites (O16, F18, C23, and C24) that are the focus of this work, two hydrogen-bond acceptors and two representative ring sites, are circled for emphasis. Carbon sites C23 (blue) and C24 (red) are distinguished with different colors. Note that the number following the atom type indicates the atom index within C153; results for the remaining atoms are provided in the Supporting Information.

ΔE σm

The C153 molecule is described atomistically and is fully flexible except that all C−H covalent bonds were constrained to their equilibrium values using SHAKE. Force field parameters for C153 were taken from the generalized Amber force field (GAFF)45 which has been used recently in simulations of C153 in ionic liquids.46 A ground-state set of charges was adopted from the work of Li et al.47 Ethanol was modeled using the optimized potentials for liquid simulations (OPLS) force field48−50 with the slight modification of using equivalent Lennard-Jones parameters for the ethanol hydrogen and the pore silanol hydrogens. This has been used in our previous work to avoid spurious solvent−pore interactions while maintaining the same liquid structure.34 All force field parameters used for C153 and ethanol can be found in the Supporting Information.

Ti

⎡ ΔE ⎤ =⎢ ⎣ σ ⎥⎦target

(3)

where ΔE = ⟨E(Ti)⟩ − ⟨E(Ti−1)⟩, σm = [σ(Ti) + σ(Ti−1)]/2, and [ΔE/σ]target can be estimated from Figure 6 of ref 55. A target acceptance probability of 20% was chosen, which has been proposed to give the best performance.55 Using this approach 23 temperatures from 298 to 688 K were found. Plots of ⟨E(T)⟩ and σ(T) as well as the temperatures used for our simulations are given in the Supporting Information. Each replica was run in the NVT ensemble using a Berendsen thermostat with a damping parameter of 0.1 ps. The same starting structure was used for each temperature. The REMD simulations consisted of 25 ns of equilibration time followed by 100 ns of data collection time for each temperature for a total of 2.3 μs of simulation time. Exchanges between temperatures are attempted every 1 ps. Cartesian coordinates were printed B

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

every 100 fs. Error bars on all calculations were calculated as 95% confidence intervals using block averaging with 20 blocks and the student t-distribution.56 Free Energy Curves and Volume Normalization. The distance, d, from the pore wall is defined as the distance to the nearest oxygen atom. This coordinate has been shown to offer a much simpler interpretation of the structure within the pore.8 However, surface roughness, which can cause a variation in the available volume at a given point in the pore, can make interpreting distributions in this coordinate difficult. The variation in volume can be accounted for by normalizing all distributions by the available free volume determined using a Monte Carlo approach.8 That is, the location of a given solute or solvent site, x, can be described by the number density ρx (d) =

Nx(d) Vx(d)

Figure 2. Density as a function of the C23 distance from the pore wall, ρC23(d), calculated after 25 (red), 50 (blue), 75 (violet), and 100 ns (black) of simulation time.

(4)

used in interpreting shorter simulations ( 7.0 Å. Using this definition we found that O16 switched from surface contact to dissolved 1015 distinct times. Given that there are favored locations on the pore surface for O16, it is interesting to consider the driving forces for this result, i.e., the structure of the pore at those locations. Using the image on the right in Figure 3, which shows the underlying pore atoms with circles highlighting the most likely positions, we see that at some points O16 is interacting with surface OH groups while at others it is, perhaps surprisingly, interacting with Si atoms; the latter actually represents the most likely location. This qualitatively indicates that hydrogen bonding plays an important role in the preferential positioning of C153 but it is not the only factor. We attempt to quantify the magnitude of its role below. Role of Hydrogen Bonding. We investigated the role of hydrogen bonding in determining the C153 position distribution using the following procedure. Densities, ρx(d), for all the C153 sites, x, were calculated as outlined above. These densities were further decomposed into the contributions when a hydrogen bond was formed with the pore surface at either the carbonyl oxygen, ρOH x (d), or any of the fluorine

where Nx(d) is the number of times site x is found at d and Vx(d) is the available volume at d. Free energies as a function of position for each atom in C153 are then determined from this density as ⎡ ρ (d ) ⎤ ⎥ ΔAx (d) = −kBT ln⎢ x ⎢⎣ ρx (do) ⎥⎦

(5)

where ρx(d0) is the density at a chosen reference distance d0. The free energy for each atom was calculated at all temperatures, and then the entropy was extracted from the relationship ⎛ ∂ΔAx (d) ⎞ ΔSx(d) = −⎜ ⎟ ⎝ ∂T ⎠V , N

(6)

A representative plot of ΔAx(d) depicting its linearity with respect to temperature is shown in the Supporting Information. Fits were performed using all but the highest three temperatures. Then the internal energy, ΔUx(d), was obtained as ΔUx(d) = ΔAx (d) + T ΔSx(d)

(7)

with uncertainties determined by error propagation.

■

RESULTS AND DISCUSSION In this section we report the results of the REMD simulations of C153 in a hydrophilic silica pore. In particular, the position and orientation distributions are presented with the corresponding free energy surfaces. REMD Convergence. Before examining the C153 distributions, it is useful to consider the convergence of the replica exchange MD simulations. To test the convergence, we calculated the volume normalized distribution of distances to the nearest pore oxygen after 25, 50, 75, and 100 ns for the C23 atomistic site. The results are shown in Figure 2. After 25 ns we see distinct peaks near 4 and 8 Å; however upon further averaging the magnitude of these peaks is greatly diminished. Similarly there is a peak near 6.5 Å that is not observed after only 25 ns. Sampling of large molecules dissolved in nanoconfined solvents represents an obvious challenge as our distributions continue to change up to 100 ns of simulation; however the curves remain within error bars after 50 ns. We emphasize that an enhanced sampling approach, such as REMD, combined with long simulation times is necessary to collect meaningful statistics,40 suggesting that caution should be C

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 3. Contact map of the C153 carbonyl oxygen on the pore surface. The surface indicates the residence time at that point during the simulation. This residence time is indicated, in decreasing order, from blue to white to red. The image on the right depicts the same part of the pore but shows the pore atoms beneath the surface. Circles are drawn to highlight the most likely positions of C153 on the pore.

Figure 4. Density as a function of the distance from the pore wall, ρx(d), for atoms (a) C23, (b) C24, (c) O16, and (d) F18 (shown in black). For each position the time a hydrogen bond is formed with the pore wall at the carbonyl oxygen (blue) or any of the fluorines (red) is also depicted. All curves are normalized by the available free volume.

atoms, ρFH x (d). All curves are normalized by the available free volume as detailed earlier and defined so that ρx(d) = ρOH x (d) + FH no HB OH−FH no HB (d) − ρx (d), where ρx (d) and ρx (d) + ρx (d) are the densities for C153 when no hydrogen ρOH−FH x bond with the pore is present and when hydrogen bonds are simultaneously formed with the pore at both the O16 and F sites, respectively. Hydrogen bonds are determined using a geometric definition of ROO,OF ≤ 3.5 Å, rHO,HF ≤ 2.4 Å, and θHOO,HFO ≤ 20°. Results of this calculation for C23, C24, F18, and O16 are shown in Figure 4. Similar figures for all of the other atoms can be found in the Supporting Information.

All sites depicted here have a dominant peak near 3 Å which more conclusively indicates that C153 preferentially sits near the pore surface. The C23 site, Figure 4a, has a second peak around 6.5 Å, while C24 and F18, Figure 4b,d, have shoulders that extend past their main peaks near the pore surface. The distribution for the O16 site has a peak that is closer to the pore wall than all of the other sites and its peak is significantly narrower, with no secondary peaks. Note that hydrogen bonding with the fluorine atoms is nearly negligible. The secondary features for C23, C24, and F18 can be attributed to hydrogen bonding at the O16 site. That is, the peak at 6.5 Å for C23 is almost entirely accounted for by D

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

hydrogen bonding as ρC23(d) ≈ ρOH C23(d) in that region and the shoulders for C24 and F18 are also located where ρOH x (d) peaks. This suggests that hydrogen bonding with the pore is important in determining the favored locations and orientations; however we observe that hydrogen bonding plays a more diminished role at other locations in the position distributions. For example, the carbonyl oxygen has just one peak in its density profile at 2.8 Å where hydrogen bonding constitutes only ∼50% of the position peak height. The O16 density has an extended tail that cannot be accounted for by hydrogen bonding to the pore surface given our hydrogen-bonding criteria (however, see later text). Similarly, when the other sites are near the pore wall, there is significantly less hydrogen bonding at the O16 site. In total, the C153 molecule is hydrogen-bonded to the pore wall at the O16 site during ∼34% of the simulation. Thus, pore hydrogen bonding is a key component determining the C153 position, but other interactions, such as solvent interactions (see later text), are also playing a critical role. We have shown that C153 is preferentially found near the pore surface and that hydrogen bonding is a key factor. However, we can also use the hydrogen-bonding analysis to qualitatively infer possible orientations of the solute with respect to the pore wall. Specifically, C23 and O16 are on opposite sides of the molecule, yet when C23 is far from the pore wall, we see that hydrogen bonding with the pore is likely at the O16 site; i.e., O16 is near the pore wall. A relatively perpendicular arrangement with the pore wall would be necessary for this effect to be possible. We also see that both O16 and C23 have peaks near the pore wall while only C23 has a peak further away from the wall. Therefore, it is likely that the molecule is more parallel to the pore when C23 is close to the surface or else O16 would have a corresponding peak at large d. This indicates that there are two molecular orientations of the C153 molecule next to the pore surface. Solute Orientational Distribution. A more quantitative approach was taken by considering the angle θ that the solute forms with the pore wall. An angle with the pore wall was calculated between the molecular axis and a radial vector that extends from the pore axis to the pore wall through the O16 site. The molecular axis 1 from Figure 5 was used for this

analysis. The probability distribution of cos θ and the O16 distance from the pore wall is shown in Figure 6 alongside vertical slices of the three main peaks from the full distributions. When O16 is close to the pore wall (d = 2.7 Å) a sharp peak is present at cos θ = 1 which suggests an almost perpendicular orientation. However, when O16 moves farther from the surface (d = 3.2 Å), the main peak shifts to cos θ = 0.5, θ = 60°, indicating a more parallel orientation relative to the pore wall. As O16 moves into the interior (d = 5.24 Å) the main peak returns to cos θ = 1; however beyond that point the distribution has a long flat tail that extends all the way to cos θ = −1. This change in orientational preference shows a strong correlation between position and orientation, the origins of which are discussed in more detail in later text. Free Energy Profiles. With the position and orientation densities of C153 in hand, we considered the thermodynamic driving forces that shape them. Free energy curves as a function of distance from the pore wall were calculated from the position densities, ρx(d), as defined in eq 4. This was done for all atomic sites at all temperatures. Figure 7 depicts the positiondependent density and corresponding free energy profile for C23 at a selected set of temperatures (the importance of free volume normalization is evident in this figure as the position distribution at 298 K would go to zero in the interior of the pore due to the diminishing volume at large d). These data show that a higher temperature increases the likelihood of finding C23 in the interior of the pore. Calculated free energies, internal energies, and entropies for C23, C24, O16, and F18 are shown in Figure 8. Entropy is plotted as −TΔSx to emphasize its contribution to the total free energy. The C23 site has two free energy minima, separated by a barrier of ∼1.3 kcal/mol, while the other three sites exhibit only one. The barrier to move into the pore is the largest for the O16 site at ∼2.5 kcal/mol, while the C24 and F18 sites have barriers of ∼2.0 and ∼1.5 kcal/mol. The calculated barriers are relatively small which is surprising given the known slow dynamics of the C153 molecule. This suggests that the long time scales for C153 rotation and diffusion is better attributed to the inhibited dynamics of the solvent rather than the free energy barriers. The effect of hydrogen bonding can be seen by comparing Figure 8 to Figure 4 as the minima in the internal energy for C23 at ∼6 Å and at ∼4 Å for C24 and F18 correlate with significant hydrogen bonding. We note that a second minimum in the internal energy for O16 appears at ∼3.2 Å which is likely caused by hydrogen bonding with the solvent. We revisit this point in later text, but it remains consistent with the suggestion made earlier that interactions with the solvent become important as O16 moves away from the pore wall. The free energy profiles reflect a competition between the internal energy and entropy. For example, for the O16 site, Figure 8c, the internal energy has a broad minimum that extends, within error bars, from 2.8 to 4.0 Å. However, the free energy has a well-defined minimum at 2.8 Å, an effect created by an increase in −TΔSO16 at larger d values. A similar effect is seen for C24, Figure 8b, for which there is a minimum in −TΔSC24 and ΔAC24 at ∼3.2 Å, whereas the minimum in the internal energy is at ∼4.2 Å. In general, ΔUx has a broad minimum or two minima that are close in energy. These are created by two possible orientations with respect to the pore wall and two hydrogen-bonding states, namely, to the pore or the solvent (see later text). The entropy favors one orientation such that ΔAx has a single, well-defined minimum.

Figure 5. Molecular axes used to calculate the angle θ with the pore wall. θ represents the angle between the chosen molecular axis and a vector that extends from the pore axis to the pore wall through a given atomistic site. E

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 6. Normalized probability distribution of cos θ for axis 1 at each distance from the pore wall for the O16 atomistic site.

■

DISCUSSION Orientational Entropy. It is interesting to consider the origins of the entropic contributions. We have previously found that entropy, and rotational entropy in particular, was a key factor in determining the location of model solutes in simpler confining frameworks.59 Thus, the important role of entropy here is not totally unexpected. It is interesting then to compare the nature of these entropic effects in this more complex system. The orientational entropy can be estimated using60 1

ΔSorient(d) = −kB

∫−1 P(cos θ ; d) ln P(cos θ ; d) d cos θ (8)

where P(cos θ;d) is the probability distribution in cos θ, displayed in Figure 6, for a molecule to have an orientation defined by cos θ when it is at a distance d from the nearest pore oxygen. Both molecular axes from Figure 5 were considered here for the calculation of θ. Rotational entropies using the two molecular axes with the total entropy are shown in Figure 9. The choice of molecular axis used to calculate θ appears to have little effect on the estimated orientational entropy. We observe that ΔSorient(d) is always at a minimum (i.e., −TΔSorient increases) near the pore wall, again caused by the pore wall restricting the rotation of the large C153 molecule. For the O16 site the orientational entropy is flat except for this change at the pore wall. However, for C23, ΔSorient(d) has another minimum at d ≈ 5.8 Å, a point at which hydrogen bonding is probable. This suggests that hydrogen bonding can restrict the molecules configuration, thereby decreasing the orientational entropy. Also apparent is that orientational entropy is a dominant factor in the total entropy for some portions of the molecule; especially near the pore wall. Most notably for C23, Figure 9a, where ΔSorient(d) closely matches the total entropy up to almost d = 7 Å. However, the orientational entropy does not correlate well with the total entropy for O16. Another surprising result is that orientational entropy does not increase like the total entropy in the interior of the pore. One might expect that as the molecule moves toward the interior there should be more room to rotate and thus we might see in increase in rotational entropy. Instead these results suggest different possibilities, e.g., increased translational entropy of C153 or reduced solvent ordering when C153 is in the pore interior. Solvent Interactions. The broad position distribution and the disagreement between rotational entropy and total entropy

Figure 7. (a) Position densities, ρC23(d), and (b) free energy, ΔAC23(d), for the C23 site at T = 298 (black), 357 (red), 436 (blue), and 529 K (violet).

All sites have a minimum in −TΔSx near the pore wall; however closest to the wall we see an increase in −TΔSx. This is an intuitive result as the large C153 molecule is constrained by the curvature of the pore wall. We stress that entropy plays a large role in determining the overall free energy curve, particularly near the pore wall, where there is considerable excluded volume. Crowding has been shown to increase the entropic contributions to protein stabilization,57,58 and it is possible that an analogous effect is observed here. However, what remains to be explored is the discrepancy between the location of the minimum in the free energy and the entropy of O16. Thus, below we examine the role of different types of entropies contributing to the total entropy. F

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 8. Free energy, ΔAx(d) (red), internal energy, ΔUx(d) (black), and entropy, plotted as −TΔSx(d) (blue), are shown for C153 atoms (a) C23, (b) C24, (c) O16, and (d) F18. The −TΔSx(d) and ΔUx(d) are shifted by 1 and 2 kcal/mol, respectively, for clarity.

■

SUMMARY AND CONCLUDING REMARKS We have used replica exchange molecular dynamics to determine the position and orientation of the widely used C153 fluorescent probe confined in a ∼2.4 nm hydrophilic silica pore. The C153 was dissolved in ethanol and modeled in its electronic ground state. The present analysis provides detailed insight into the properties of a C153 molecule before excitation in a time-dependent fluorescence experiment. More generally, the results illustrate general features of the behavior of a large solutes in nanoconfined liquids. The simulations found that, while the C153 carbonyl oxygen is almost always found at the pore surface, other atomistic sites of the molecule can be found near the surface as well. Two primary orientations were thus indicated: one in which the carbonyl oxygen points at the pore wall with C153 lying almost perpendicular and another where the carbonyl oxygen moves away from the pore wall and the whole molecule lies relatively flat against the surface. The free energy barriers between these orientations were found to be comparatively small, a result that is particularly surprising in light of the slow dynamics of confined C153. An interesting connection to TDF experiments can be made based on this result which suggests that displacement of solvent molecules can play an important role. We have shown that for simpler systems solute diffusion can occur after the excitation of fluorescent probes.24,32−34 However, the C153 molecules in these two orientations would likely diffuse at different rates. The flat orientation has a larger effective surface area to move into the interior than the perpendicular arrangement and would, therefore, perhaps diffuse more slowly. Thus, if solute diffusion is observable as an additional time scale, the orientation of the solute at the time of excitation must also

of O16 indicate that additional interactions are important. Earlier we speculated that interactions with the solvent might play a role, and we examine that hypothesis in greater detail here. Expanding upon the data in Figure 4, we calculated for a given position, d, the fraction of time a hydrogen bond is formed at the carbonyl oxygen with an ethanol molecule in addition to hydrogen bonds formed with the pore. Subsequently we are able to determine the percentage of time that a hydrogen bond is not formed with either the solvent or pore wall. The results of this calculation are shown in Figure 10. Interactions with the solvent are clearly important as the tail in the position distribution that extends from 3 to 4 Å is almost entirely represented by hydrogen bonding to ethanol, and this creates a minimum in the internal energy at ∼3.25 Å. This also matches the region at which the cos θ distribution changes from having a perpendicular arrangement with the pore wall to a parallel arrangement (Figure 6). Thus, it is reasonable to suggest that hydrogen bonds with the ethanol energetically stabilize the O16 atom at large distances from the pore wall as the entire C153 molecule lays more parallel to the pore wall. However, this minimum at ∼3.25 Å is not observed in the free energy due to an increase in −TΔSO16 which also corresponds to the point at which the position density is entirely accounted for by hydrogen bonding with the solvent. An additional minimum in the internal energy exists at ∼2.9 Å, which is near the free energy minimum. This occurs at a point where hydrogen bonding with the pore becomes important. These results suggest that hydrogen bonding can have a substantial effect on the internal energy, entropy, and free energy and that perhaps a range of hydrogen-bonding motifs at the pore wall are important. G

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

The two orientations differ in the identity of the hydrogenbonding donor to the carbonyl oxygen; the perpendicular orientation is hydrogen-bonded to the pore while the parallel orientation is hydrogen-bonded to the ethanol. This has broader implications for design criteria of porous materials. While hydrogen bonding at the pore surface clearly plays an important role in the position of reactant complexes, the solvent can also be used to control position. Moreover we have shown that entropy near the pore surface has a large influence on the free energy landscape. This effect is created by excluded volume, similar to what is seen in protein stabilization by crowding. This effect could be tuned by varying the roughness of the pore to create more pockets at the pore surface. Simulations have revealed that pore heterogeneity can strongly affect the dynamics of confined water,61,62 suggesting another possible design route. In order for solute diffusion to play a role in time-dependent fluorescence experiments, the ground- and excited-state equilibrium position (or orientation) distributions must differ. A condition which, if met, would cause the solute to diffuse or reorient itself from the ground-state to the excited-state distribution after excitation. Thus, analogous simulations are currently underway for the C153 excited electronic state. In addition, the effect of pore chemistry on the dye properties is being examined through simulations with hydrophobic pores and ones with no hydrogen-bond donors, to connect with previous experiments.21

■

Figure 9. Orientational entropy, ΔSorient(d), for atoms (a) C23 and (b) O16. Two different molecular axes, shown in Figure 5, were considered in calculating θ. The total entropy (black) is included for comparison.

ASSOCIATED CONTENT

S Supporting Information *

Text describing additional information for the study and accompanying references, tables listing charges and LennardJones parameters, bond and angle parameters, dihedral parameters, GAFF atom types, charges, and Lennard-Jones parameters for C153 atoms, and bonding, angular, and dihedral parameters used for C153, and figures showing the 2D structure of coumarin 153, ⟨E(T)⟩ and σ(T) calculated at various tempteratures, ΔAx(d) as a function of temperature, a contact map for the half of the pore not depicted in the paper, and position distributions, free energy, internal energy, and entropy curves for all other atoms not depicted here. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

Figure 10. Bottom: The total density as a function of O16 position, ρO16(d) (black) along with its decomposition into contributions where a hydrogen bond is formed with the pore wall (red) or with the solvent (blue) and where no hydrogen bond is formed (violet). Top: Same as Figure 8, showing the alignment between the internal energy minima (indicated by vertical dashed lines) as well as −TΔS and the peak in the density for O16 with no hydrogen bond. Energies are in kilocalories per mole.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Support for this work was provided by the Chemical Sciences, Geosciences, and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U.S. Department of Energy (Grant No. DE-FG02-05ER15708). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1053575. All calculations were performed on the supercomputer “Kraken” located at the National Institute for Computational Science at the University of Tennesse Knoxville.

be considered. This possibly represents a new difficulty in interpreting results and assigning the origins of time scales from TDF experiments of confined liquids. Preliminary results indicate that both orientations exist in the excited state which suggests a conversion between the two states might occur after excitation. In this context, it will be interesting to consider whether the dynamics depend on the starting orientation of the molecule. H

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

■

Article

(23) Shirota, H.; Segawa, H. Solvation Dynamics of Formamide and N,N-Dimethylformamide in Aerosol OT Reverse Micelles. Langmuir 2004, 20, 329−335. (24) Thompson, W. H. Simulations of Time-Dependent Fluorescence in Nano-Confined Solvents. J. Chem. Phys. 2004, 120, 8125− 8133. (25) Yamaguchi, A.; Amino, Y.; Shima, K.; Suzuki, S.; Yamashita, T.; Teramae, N. Local Environments of Coumarin Dyes within Mesostructured Silica-Surfactant Nanocomposites. J. Phys. Chem. B 2006, 110, 3910−3916. (26) Kumbhakar, M.; Goel, T.; Nath, S.; Mukherjee, T.; Pal, H. Microenvironment in the Corona Region of Triblock Copolymer Micelles: Temperature Dependent Solvation and Rotational Relaxatioin Dynamics of Coumarin Dyes. J. Phys. Chem. B 2006, 110, 25646−25655. (27) Kumbhakar, M.; Ganguly, R. Influence of Electrolytes on the Microenvironment of F127 Triblock Copolymer Micelles: A Solvation and Rotation Dynamics Study of Coumarin Dyes. J. Phys. Chem. B 2007, 111, 3935−3942. (28) Kamijo, T.; Yamaguchi, A.; Suzuki, S.; Teramae, N.; Itoh, T.; Ikeda, T. Solvation Dynamics of Coumarin 153 in Alcohols Confined in Silica Nanochannels. J. Phys. Chem. A 2008, 112, 11535−11542. (29) Das, A.; Chakrabarti, J. Microscopic Mechanisms of Confinement-Induced Slow Solvation. J. Phys. Chem. A 2013, 117, 10571− 10575. (30) Wong, M.; Thomas, J. K.; Grätzel, M. Fluorescence Probing of Inverterd Micelles. The State of Solubilized Water Clusters in Alkane/ Diisoocytl Sulfosuccinate (Aerosol OT) Solution. J. Am. Chem. Soc. 1976, 98, 2391−2397. (31) Bhattacharyya, K.; Bagchi, B. Slow Dynamics of Constrained Water in Complex Geometries. J. Phys. Chem. A 2000, 104, 10603− 10613. (32) Thompson, W. H. A Monte Carlo Study of Spectroscopy in Nanoconfined Solvents. J. Chem. Phys. 2002, 117, 6618−6628. (33) Thompson, W. H. Solvation Dynamics and Proton Transfer in Nanoconfined Liquids. Annu. Rev. Phys. Chem. 2011, 62, 599−619. (34) Vartia, A. A.; Thompson, W. H. Solvation and Spectra of a Charge Transfer Solute in Ethanol Confined Within Nanoscale Silica Pores. J. Phys. Chem. B 2012, 116, 5414−5424. (35) Levinger, N. E.; Swafford, L. A. Ultrafast Dynamics in Reverse Micelles. Annu. Rev. Phys. Chem. 2009, 60, 385−406. (36) Fayer, M. D. Dynamics of Water Interacting with Interfaces, Molecules, and Ions. Acc. Chem. Res. 2012, 45, 3−14. (37) Elola, M. D.; Rodriguez, J.; Laria, D. Liquid Methanol Confined within Functionalized Silica Nanopores. 2. Solvation Dynamics of Coumarin 153. J. Phys. Chem. B 2011, 115, 12859−12867. (38) Elola, M. D.; Rodriguez, J. Solvation of Coumarin 480 within Nano-Confining Environments: Structure and Dynamics. J. Chem. Phys. 2014, 140, 034702. (39) Guchhait, B.; Biswas, R. Solute and Solvent Dynamics in Confined Equal-Sized Aqueous Environments of Charged and Neutral Reverse Micelles: A Combined Dynamic Fluorescence and All-Atom Molecular Dynamics Simulation Study. J. Phys. Chem. B 2013, 117, 3345−3361. (40) Palonćyová, M.; Berka, K.; Otyepka, M. Convergence of Free Energy Profile of Coumarin in Lipid Bilayer. J. Chem. Theory Comput. 2012, 8, 1200−1211. (41) Brodka, A.; Zerda, T. W. Properties of Liquid Acetone in Silica Pores: Molecular Dynamics Simulations. J. Chem. Phys. 1996, 104, 6319−6326. (42) Gulmen, T. S.; Thompson, W. H. Dynamics in Small Confining Systems VIII. In Materials Research Society Symposium Preceedings 899E; Fourkas, J. T., Levitz, P., Overney, R., Urbakh, M., Eds.; Materials Research Society: Warrendale, PA, USA, 2005. (43) Gulmen, T. S.; Thompson, W. H. Testing a Two-State Model of Nanoconfined Liquids: Conformational Equilibrium of Ethylene Glycol in Amorphous Silica Pores. Langmuir 2006, 22, 10919−10923. (44) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System With

REFERENCES

(1) Gordillo, M. C.; Martí, J. Hydrogen Bond Structure of Liquid Water Confined in Nanotubes. Chem. Phys. Lett. 2000, 329, 341−345. (2) Shao, Q.; Huang, L.; Zhou, J.; Lu, L.; Zhang, L.; Lu, X.; Jiang, S.; Gubbins, K. E.; Zhu, Y.; Shen, W. Molecular Dynamics Study on Diameter Effect in Structure of Ethanol Molecules Confined in SingleWalled Carbon Nanotubes. J. Phys. Chem. C 2007, 111, 15677−15685. (3) Morales, C. M.; Thompson, W. H. Simulations of Infrared Spectra of Nanoconfined Liquids: Acetonitrile Confined in Nanoscale, Hydrophilic Silica Pores. J. Phys. Chem. A 2009, 113, 1922−1933. (4) Chowdhary, J.; Ladanyi, B. M. Molecular Dynamics Simulation of Aerosol-OT Reverse Micelles. J. Phys. Chem. B 2009, 113, 15029− 15039. (5) Milischuk, A. A.; Ladanyi, B. M. Polarizability Anisotropy Relaxation in Nanoconfinement: Molecular Simulation Study of Acetonitrile in Silica Pores. J. Phys. Chem. B 2013, 117, 15729−15740. (6) Rivera, C. A.; Bender, J. S.; Manfred, K.; Fourkas, J. T. Persistence of Acetonitrile Bilayers at the Interface of Acetonitrile/ Water Mixtures with Silica. J. Phys. Chem. A 2013, 117, 12060−12066. (7) Moutain, R. D. Molecular Dynamics Simulation of WaterAcetonitrile Mixtures in a Silica Slit. J. Phys. Chem. C 2013, 117, 3923− 3929. (8) Thompson, W. H. Structure, Dynamics and Hydrogen Bonding of Acetonitrile in Nanoscale Silica Pores. Mol. Simul. 2014, http://dx. doi.org/10.1080/08927022.2014.926550. (9) D’Orazio, F.; Bhattacharja, S.; Halperin, W. P. Enhanced SelfDiffusion of Water in Restricted Geometry. Phys. Rev. Lett. 1989, 63, 43−46. (10) Koone, N.; Shao, Y.; Zerda, T. W. Diffusion of Simple Liquids of Porous Sol−Gel Glass. J. Phys. Chem. 1995, 99, 16976−16981. (11) Kittaka, S.; Iwashita, T.; Serizawa, A.; Kranishi, M.; Takahara, S. Low Temperature Properties of Acetonitrile Confined in MCM-41. J. Phys. Chem. B 2005, 109, 23162−23169. (12) Norton, C. D.; Thompson, W. H. On the Diffusion of Acetonitrile in Nanoscale Amorphous Silica Pores. Understanding Anisotropy and the Effects of Hydrogen Bonding. J. Phys. Chem. C 2013, 117, 19107−19114. (13) Liu, G.; Li, Y.; Jonas, J. Confined Geometry Effects on Reorientational Dynamics of Molecular Liquids in Porous Silica Glasses. J. Chem. Phys. 1991, 95, 6892−6901. (14) Sarkar, N. Solvation Dynamics of Coumarin 480 in Reverse Micelles. Slow Relaxation of Water Molecules. J. Phys. Chem. 1996, 100, 10523−10527. (15) Pant, D.; Riter, R. E.; Levinger, N. E. Influence of Restricted Environment and Ionic Interactions on Water Solvation Dynamics. J. Chem. Phys. 1998, 109, 9995−10003. (16) Riter, R. E.; Willard, D. M.; Levinger, N. E. Water Immobilization at Surfactant Interfaces in Reverse Micelles. J. Phys. Chem. B 1998, 102, 2705−2714. (17) Senapati, S.; Chandra, A. Molecular Dynamics Simulations of Simple Dipolar Liquids in Spherical Cavity: Effects of Confinement on Structural, Dielectric, and Dynamical Properties. J. Chem. Phys. 1999, 111, 1223−1230. (18) Shirota, H. Solvation Dynamics in Nonaqueous Reverse Micelles. J. Phys. Chem. B 1999, 103, 1437−1443. (19) Pal, S. K.; Sukul, D.; Mandal, D.; Sen, S.; Bhattacharyya, K. Solvation Dynamics of Coumarin 480 in Sol−Gel Matrix. J. Phys. Chem. B 2000, 104, 2613−2616. (20) Baumann, R.; Ferrante, C.; Deeg, F. W.; Bräuchle, C. Solvation Dynamics of Nile Blue in Ethanol Confined in Porous Sol-Gel Glasses. J. Chem. Phys. 2001, 114, 5781−5791. (21) Baumann, R.; Ferrante, C.; Kneuper, E.; Deeg, F.-W.; Bräuchle, C. Influence of Confinement on the Solvation and Rotational Dynamics of Coumarin 153 in Ethanol. J. Phys. Chem. A 2003, 107, 2422−2430. (22) Chowdhury, P. K.; Halder, M.; Sanders, L.; Arnold, R. A.; Liu, Y.; Armstrong, D. W.; Kundu, S.; Hargrove, M. S.; Song, X.; Petrich, J. W. The Complex of Apomyoglobin with the Fluorescent Dye Coumarin 153. Photochem. Photobiol. 2004, 79, 440−446. I

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (45) Wang, J.; Wolf, R.; Caldwell, J.; Kollman, P.; Case, D. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (46) Terranova, Z. L.; Corcelli, S. A. On the Mechanism of Solvation Dynamics in Imidazolium-Based Ionic Liquids. J. Phys. Chem. B 2013, 117, 15659−15666. (47) Li, H.; Arzhantsev, S.; Maroncelli, M. Solvation and Solvatochromism in CO2-Expanded Liquids. 2. Experiment−Simulation Comparisons of Preferential Solvation in Three Prototypical Mixtures. J. Phys. Chem. B 2007, 111, 3208−3221. (48) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc. 1984, 106, 6638−6646. (49) Jorgensen, W. L. Optimized Intermolecular Potential Functions for Liquid Alcohols. J. Phys. Chem. 1986, 90, 1276−1284. (50) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. An All Atom Force Field for Simulations of Proteins and Nucleic Acids. J. Comput. Chem. 1986, 7, 230−252. (51) Swendsen, R. H.; Wang, J. Replica Monte Carlo Simulation of Spin-Glasses. Phys. Rev. Lett. 1986, 57, 2607−2609. (52) Hukushima, K.; Nemoto, K. Exchange Monte Carlo Method and Application to Spin Glass Simulations. J. Phys. Soc. Jpn. 1996, 65, 1604−1608. (53) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 199, 314, 141−151. (54) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (55) Rathore, N.; Chopra, M.; de Pablo, J. J. Optimal Allocation of Replicas in Parallel Tempering Simulations. J. Chem. Phys. 2005, 122, No. 023111. (56) Garland, C. W.; Nibler, J. W.; Shoemaker, D. P. Experiments in Physical Chemistry; McGraw-Hill: New York, NY, USA, 2009. (57) Minton, A. P. Implications of Macromolecular Crowding for Protein Assembly. Curr. Opin. Struct. Biol. 2000, 10, 34−39. (58) Cho, S. S.; Reddy, G.; Straub, J. E.; Thirumalai, D. Entropic Stabilization of Proteins by TMAO. J. Phys. Chem. B 2011, 115, 13401−13407. (59) Mitchell-Koch, K. R.; Thompson, W. H. How Important is Entropy in Determing the Position-Dependent Free Energy of a Solute in a Nanoconfined Solvent. J. Phys. Chem. C 2007, 111, 11991− 12001. (60) Ashbaugh, H. S.; Paulaitis, M. E. Entropy of Hydrophobic Hydration: Extension to Hydrophobic Chains. J. Phys. Chem. 1996, 100, 1900−1913. (61) Laage, D.; Thompson, W. H. Reorientation Dynamics of Nanoconfined Water: Power-Law Decay, Hydrogen-Bond Jumps and Test of a Two-State Model. J. Chem. Phys. 2012, 136, No. 044513. (62) Fogarty, A. C.; Duboué-Dijon, E.; Laage, D.; Thompson, W. H. Origins of the Non-Exponential Reorientation Dynamics of Nanoconfined Water. J. Chem. Phys. 2014, 141, No. 18C523.

J

dx.doi.org/10.1021/jp509051n | J. Phys. Chem. B XXXX, XXX, XXX−XXX