Solute location in a nanoconfined liquid depends on charge distribution

Mesoporous materials have attracted intense interest for applications such as catalysis, drug delivery, and separations.^1–8 Many of these systems involve processes that include charge transfer or a change in the electronic distribution, for example, mesoporous materials used for photo-activated release of drugs,² solar energy conversion and storage,^5,6 or catalytic frameworks.^6,7 An improved picture of the effects of changes in molecular charge distributions within nanoconfined liquids is thus key to advancing the rational design of many applications.

At the same time, these systems represent an exciting platform for exploring fundamental issues of complex liquid dynamics that emerge at interfaces and in confinement.^9–11 An impressive diversity of frameworks has been synthesized that can confine liquids on nanometer length scales, including silicate glasses, reverse micelles, metal-organic frameworks, biological systems, and supramolecular assemblies. Each of these systems presents varied and tunable characteristics—e.g., size, shape, surface chemistry, flexibility—that modify the structure and dynamics of the confined liquid. Indeed, this variety is matched by the broad range of properties affected by nanoconfinement, such as molecular diffusion and reorientation, solvation dynamics, conformational equilibria, and reaction mechanisms.

The effects of nanoscale confinement on chemistry are thus substantial and numerous. A key class of reactions is that involving charge transfer, including proton and electron transfer reactions. These are important processes in a wide range of synthetic chemistry as well as in biological systems. At the same time, they are also strongly affected by nanoconfinement as the reaction coordinate is a collective one involving the reorientation and translation of the molecules surrounding the reaction complex, motions that are strongly modified in a confined liquid. Indeed, proton transfer rate constants, for example, are generally (though not always) found to be reduced in nanoconfined solvents.¹¹ While many factors that influence these changes have been discussed—e.g., reduced solvent polarity, greater effective viscosity, slowed reorientation of interfacial solvent, and interactions with surface groups—a clear, predictive picture of charge transfer reactivity in nanoconfined solvents has yet to emerge.

In this paper, we examine one such factor that is still not well understood: the location of the reaction complex within the nanoconfining framework.^11–21 We have previously found in simple models that the position distributions of a proton transfer reactant complex and product complex differ and that this difference can play a role in the reaction mechanism.^14,15 Interestingly, Sedgwick et al. have observed different locations of a photoacid in reverse micelles constructed with different surfactant molecules;¹⁸ in reverse micelles, a solute can move not only within the water pool but also intercalate into the surfactant layer that represents the confining framework. However, we are unaware of any measurements that probe the relative positions before and after acid ionization of a proton transfer complex (or photoacid).

One way to model the fate of proton (or electron) transfer reactants and products is through a dye molecule with a charge transfer electronic transition; such a system can be easier to both probe experimentally and simulate computationally. Typically, the ground electronic state has a low dipole moment similar to an acid-base, A–H⋯B, reactant complex while the excited electronic state can exhibit a large dipole moment (e.g., >10 D) analogous to the product complex, A⁻⋯H–B⁺. The question then becomes are the ground- and excited-state position distributions of such a dye molecule in a nanoconfined solvent different? As just outlined, the answer has implications beyond the behavior of an artificial dye molecule and is relevant to a wide scope of chemical reactions that involve charge transfer, i.e., reactants and products that differ significantly in molecular charge distribution.

It is now well-established that the structure and dynamics of a nanoconfined liquid are often dramatically modified relative to the bulk liquid and vary strongly with location in the confining framework, even over length scales of only a few angstroms.^{9–11,22–24} Thus, even small changes in reaction complex location can have important consequences for the reaction mechanism, barriers, and rate constants that are relevant to applications of mesoporous materials, for example, for catalysis or sensing. This provides a strong impetus to understand the relationship between solute charge distribution and location in nanoconfined solvents, including the key factors that influence this connection.

In the remainder of this paper, we describe the results of atomistic simulations of a coumarin 153 (C153) dye molecule dissolved in ethanol confined in an approximately 2.4 nm hydrophilic (–OH terminated) silica pore. The focus is on examining how a change in the charge distribution of a solute molecule affects its position. This system was the subject of a previous experimental study,²⁵ but information about the solute position was not obtained. The C153 molecule possesses the properties required to model a charge transfer process as discussed above. In particular, it possesses a charge transfer electronic transition that increases its dipole moment from 6 D in the ground state to 13 D in the excited state.^26–28 In previous simulations, we have shown that ground-state C153 is preferentially localized at the pore interface.²⁹ Here, we compare those results to simulations of the excited-state solute. We find that the location of the molecule is not the same for the two charge distributions and we examine the generality and origin of these differences in detail.

The system considered here, illustrated in Fig. 1, is an amorphous ∼2.4 nm diameter silica pore model³⁰ filled with 121 ethanol molecules and 1 C153 molecule. The pore framework consists of 1042 Si atoms and 2144 O atoms which are held fixed during the simulation. The pore surface is functionalized with 36 silanol, SiOH, groups, and 6 geminal, Si(OH)₂, groups. Surface O–H bonds are held fixed using the SHAKE algorithm³¹ but the OH groups are allowed to rotate. The Si–O–H angle is described by a harmonic potential. The C153 molecule is described atomistically and is fully flexible except that SHAKE is used to hold the C–H bond lengths constant. Force field parameters for the pore, provided in the supplementary material,⁴⁷ are described in detail elsewhere,^30,32,33 while parameters for the solvent and solute are adapted from the Optimized Potentials for Liquid Simulations (OPLS) force field^34–36 and the Generalized AMBER Force Field (GAFF),³⁷ respectively. Ethanol was added to the pore via a grand canonical Monte Carlo simulation and the C153 molecule was inserted in a random position and orientation followed by removal of any overlapping ethanol molecules. The final system used in all simulations consists of 121 ethanol molecules and 1 C153 molecule all contained within a 44 × 30 × 40 ų simulation box. Atomistic models for C153 were taken from the work of Maroncelli et al.³⁸ The ground- and excited-state parameters differ only in the atomic charges, giving μ_g = 6 D and μ_e = 13 D.^26–28 

Replica exchange molecular dynamics³⁹ (REMD) simulations for the excited-state C153 were carried out using LAMMPS⁴⁰ at 20 different temperatures for a total of 150 ns after a 5 ns equilibration. The results are compared to those for ground-state C153 from a 100 ns simulation, reported previously,²⁹ to explore the electronic state dependence. Replicas were chosen using the procedure of Rathore et al.,⁴¹ in which 23 temperatures were found to be needed for the ground-state simulations²⁹ and 20 temperatures for the excited-state simulations. The temperatures used are given in the supplementary material.⁴⁷ Replicas in the ground-state simulations were found to make an average of 8.3 round trips from the lowest temperature to the highest temperature and back, with an average acceptance ratio of 33.8%. Excited-state replicas made an average of 21.1 round trips with an average acceptance probability of 23.3%.

Each replica is run in the NVT ensemble using a Berendsen thermostat with a temperature damping parameter of 0.1 ps. The details of the ground-state simulations are given in Ref. 29. Briefly, ground-state simulations were initialized with an identical starting structure at each temperature and followed by 25 ns of REMD equilibration and a 100 ns REMD data collection stage. Each replica in the excited-state simulations was initialized from the final configuration of a 5 ns equilibration run at the desired temperature. Data were collected from 150 ns of REMD. A 1 fs time step was used in the ground-state simulations while a 2 fs time step was used for the excited-state. In both simulations swaps between configurations were attempted every 1 ps, and configurations were printed every 100 fs.

Insight into the location of the C153 molecule in the silica pore is provided in Fig. 2, which shows the distribution of the carbonyl oxygen in the xy plane (averaged over the position along the pore, z, axis); data for other sites on the molecule yield the same qualitative results (see the supplementary material⁴⁷) but here we focus on the carbonyl oxygen as it is particularly instructive in understanding the C153 arrangement.²⁹ Distributions are shown for both the ground and excited electronic states of C153. It is clear from Fig. 2 that the position distribution of C153 changes significantly with its electronic state and hence charge distribution. As we have noted elsewhere²⁹ and is illustrated in Fig. 2, ground-state C153 sits preferentially at the pore surface with significant affinity for a small number of sites. In contrast, excited-state C153 has a weaker attraction to the pore interface. This gives rise to both a broader distribution of positions at the pore surface and a greater likelihood of C153 being found in the pore interior. Note that while it is common (and tempting) to talk about “the” location of a solute in a nanoconfined liquid, there is generally a distribution of positions. This is true even for the more localized ground-state C153, but particularly so for the excited-state molecule.

To better understand what shapes the ground- and excited-state C153 distributions, we next focus on the underlying driving forces. In this context, it is useful to examine the (Helmholtz) free energy as a function of C153 position, measured as the distance, d, to the nearest pore oxygen atom. The free energy for the carbonyl oxygen, ΔA_O(d), is obtained from the probability distribution for the site calculated from the REMD simulation, P_O(d), and normalized by the available volume at each distance d, V_O(d), to give the number density,²⁹ 

The available volume is obtained from a Monte Carlo calculation in which the site, here the oxygen atom, is inserted into the unfilled silica pore. An attempted insertion is rejected if the atom overlaps (based on the Lennard-Jones parameters) with a pore atom and accepted otherwise. The accepted insertions are binned according to their distance, d, to yield V_O(d).

The free energy is then calculated from the density as

where d_∘ is a reference distance that sets the zero of energy. Calculating ΔA_O(d) for the different temperatures used in the REMD simulations yields the position-dependent entropy from the slope of the free energy with temperature, i.e., from ΔA_O(d) = ΔU_O(d) − TΔS_O(d), for each value of d. Finally, the same equation is used to obtain the internal energy profile, ΔU_O(d), from the free energy and entropy. This procedure²⁹ has been carried out for each site on the C153 molecule (see the supplementary material⁴⁷), but we focus on the carbonyl oxygen to illustrate the results.

The free energy, internal energy, and entropic contribution, −TΔS_O(d), are shown in Fig. 3 for both ground- and excited-state C153. The free energy surface for both electronic states has a global minimum around 2.8 Å from the interface. However, this favorable interaction of the carbonyl oxygen with the pore surface is stronger for ground-state C153 than for the excited state; the free energy required to move from the minimum to the pore interior is decreased by ∼1 kcal/mol in the excited state. Thus, excited-state C153 is three times more likely than the ground-state molecule, 25% versus 8%, to be found in the pore interior (defined as d ≥ 4 Å) in addition to its more diffuse distribution of positions at the pore surface illustrated in Fig. 2. Examination of the free energy contributions indicate that this difference can be attributed to the internal energy. The dependence of the entropy on C153 position is not substantially different in the two electronic states; the excited-state entropic contribution actually slightly disfavors C153 in the interior relative to the ground state, in opposition to the internal energy. The net effect is that electronic excitation of C153 lowers the internal (and hence free) energy in the pore interior leading to a broader distribution of positions.

It is important to consider the larger implications for the result that a change in charge distribution affects the solute location. In the system considered here, the ground- and excited-state C153 molecules differ only in the atomic partial charges. Yet, this difference in charge distribution gives rise to significant changes in the solute disposition in the pore. Such changes in the electrostatic properties of a molecule are ubiquitous in chemical processes, occurring, e.g., upon molecular binding, conformational rearrangement, and chemical reaction. Moreover, the fact that the solute position depends on its charge distribution is almost surely general, though the specific nature of the effect, i.e., how the solute position changes upon electronic excitation will surely vary from system to system.

Here, we consider the generality of a dependence of solute position distribution on its charge distribution. The main idea is that a change in the arrangement of molecular charges modifies the interactions with the confining framework so that the position distribution relative to the framework is determined by a different Hamiltonian. Consistent with the C153 system considered here, in the following, the two charge arrangements are labeled by the electronic state as “g” (ground) and “e” (excited), but the result applies for any two charge states, e.g., reactants and products of a chemical reaction.

The two solute position distributions can be expressed at a particular value r₀ of a solute molecule coordinate, r, as

where α = g or e, V_α is the potential energy function for state α, β = 1/k_bT, and q = (r, q′) is the complete set of coordinates of the system that includes the solute coordinate, r, and all others, q′. Defining the difference between the potential energy in the two states as ΔE(r, q′) = V_e(r, q′) − V_g(r, q′), we can write the excited-state distribution as

Then, denoting the partition function for state α as

we have

using the properties of the δ-function. Finally, using Eq. (3) for P_g(r₀), we have that

Note that the right-hand-side can only be zero if e^{−βΔE(r₀,q′)}(Q_g/Q_e) = 1 for the solute at position r₀. Since Q_e and Q_g are independent of r₀, this is only true for all solute molecule positions if ΔE(r₀, q′) = V_e(r₀) − V_g(r₀) is itself independent of r₀, i.e., the potentials for the two charge states depend on the solute position in exactly the same way. However, the two potential energies will not be equal as a function of molecular positions relative to the confining interface except under quite special conditions, particularly when the two states differ significantly in the molecular charge distribution. Thus, while the difference in the two position distributions may be slight for some systems with weaker solute-pore interactions, for which ΔE is naturally small because V_g and V_e are each individually small, it will in general not vanish and may be substantial for any chemical process involving a change in charge distribution.

We next investigate the origin of the change in the internal energy upon a shift in the C153 charge distribution. We have previously found that the position of ground-state C153 is strongly affected by hydrogen bonding (H-bonding) of the C153 carbonyl oxygen with the pore silanol (SiOH) groups, which drives its localization at the interface.²⁹ Thus, it is instructive to compare the H-bonding status for the two electronic states; this is done in Fig. 4(a) using a geometric H-bond definition of R_O⋯O ≤ 3.5 Å, r_H⋯O ≤ 2.45 Å, and θ_HOO ≤ 30°. The percentage of simulation time a H-bond is formed between the C153 carbonyl oxygen and a pore silanol group, an ethanol solvent molecule, or both donors simultaneously are compared along with the percentage of time that no H-bond is formed.

The results presented in Fig. 4(a) show that C153 accepts a hydrogen bond ∼80% of the time in both its ground and excited states, predominantly (∼50%) from solvent ethanol molecules. However, a dramatic difference is observed in the H-bonding with the pore surface, with which the ground-state molecule H-bonds ∼33% of the time; this decreases to only 7% in the excited state. While part of this difference is compensated by an increase in the fraction of time that excited-state C153 spends accepting H-bonds from both a silanol and ethanol, the net result is that a clear decrease in H-bonding with the pore for the excited state is observed.

The reduced affinity for H-bonding with the silica surface is likely a result of a competition with interactions with the solvent rather than a cause of the changes in the free energy profile. The C153 carbonyl oxygen charge is essentially unchanged upon excitation,³⁸ Δq = 0.0011, as illustrated in Fig. 4(c), so these effects cannot be attributed to any reduction in the H-bond accepting ability. At the same time, a compensating increase in H-bonding with ethanol would be expected if H-bonding is the driving force for the change in solute location upon excitation, but this is not observed. Instead, we can note that the change in the C153 charge distribution upon excitation primarily involves the aromatic rings as shown in Fig. 4(c) (see also Table S5⁴⁷). The electrostatic interactions of these ring sites with the pore and the solvent are best described as non-specific solvation, and the results in Fig. 4 thus implicate these as the origin of the change in position distribution upon excitation.

The role of non-specific solvation in increasing the probability of finding the C153 solute in the pore interior upon excitation was investigated further by calculating the potential energy of interaction between the entire C153 molecule and the pore atoms or the solvent molecules. These results are shown in Fig. 4(b) for ground- and excited-state C153. For the ground-state solute, the average C153-pore interaction energy is nearly equal to the C153-solvent energy. However, upon excitation, the average excited-state C153-solvent interactions are significantly more favorable than for the ground-state while the C153-pore interactions are weaker. Combined with the H-bonding analysis, these results indicate that better solvation of excited-state C153 by ethanol is the key factor in reducing the internal energy for the molecule in the pore interior relative to its ground-state form.

Finally, it is interesting to consider the thermodynamic results discussed above in the context of the known slow dynamics of large solutes within confining frameworks. While the work presented here provides free energies and densities as a function of solute position, the results do not speak to the dynamics or time scales of solute diffusion, with one exception. The calculated free energies do show that a C153 solute must diffuse after excitation from the ground electronic state to the excited state because the free energies (and densities) for the two states differ. The time scale for that diffusion must be obtained by unbiased, nonequilibrium trajectories that are currently in progress in our laboratory.

Preliminary, unbiased MD simulations of both ground- and excited-state C153 in ethanol-filled silica pores show that the solute exhibits slow dynamics, often “sticking” at the surface for up to tens of nanoseconds. It is this slow dynamics that makes replica exchange MD necessary for mapping the solute position distributions and free energies.⁴² At the same time, in light of these long (though still ill-determined) time scales for solute diffusion, it is perhaps surprising that the solute free energy landscape has relatively modest fluctuations, e.g., variations of only ∼2 kcal/mol for the ground-state C153 and ∼1 kcal/mol for the excited-state molecule. While a full explanation must await the further studies of the solute dynamics, a few comments on this issue are in order.

It is important to note that there are two key factors that affect the location of solutes: (1) the free energy surface as a function of solute position, ΔA(d), and (2) the free volume available for the solute to occupy. The former is discussed in detail above in Sec. III B. The latter naturally depends on the pore geometry, but for cylindrical pores, there is naturally more space available near the pore interface than in the pore interior. Both the free energy and volume factors are implicitly included in the raw solute probability distribution in Cartesian coordinates, P_O(x, y), in the results shown in Fig. 2. This is not the case for the free energies as a function of, d, the distance from the nearest pore oxygen atom, given in Fig. 3, where the available volume has been divided out, as shown in Eq. (1), illustrating the intrinsic attraction of the solute to the interface independent of the free volume. In addition, the volume factor, V_O(d), is not a simple function for the silica pores considered here due to the pore roughness, though it has the generic features one would expect of a cylindrical framework. The consequence is that both the free energy and the available volume favor the solute molecule at the pore interface. Thus, the apparently small free energy bias for the solute to reside at the surface is augmented by the greater free volume at the surface.

It is also useful to consider the effect of nanoconfined solvent dynamics. We are unaware of any measurements of the diffusion constants of large solutes such as C153 in nanoconfined solvents (though such results would be a welcome addition to the literature). However, self-diffusion constants of a variety of liquids in nanoscale silica pores have been reported to be smaller by factors of 3.7-8.2 compared to the bulk liquid,⁴³ which hints at even greater slowing of large solute mobility. Similarly, molecular reorientation is dramatically slowed in silica pores, often by as much as an order-of-magnitude or more.^22,44–46 These factors suggest that solute dynamics, particularly diffusion may be slower in nanoconfined liquids due to the retarded solvent motions than indicated simply by examination of the underlying free energy surfaces.

We have shown through explicit MD simulations that the position distribution of a solute in a nanoconfined solvent changes with a shift in its charge distribution. While the system considered here is a dye molecule with a charge transfer electronic transition, this result is a general one, which is illustrated by a simple derivation. Moreover, the change in solute location is shown to be associated with the change in solute interactions with the pore interface and thus will occur not only for electronic excitation but also for any change in solute properties including, e.g., charge transfer reactions such as proton- or electron-transfer reactions or even conformational changes.

The precise nature of the change in solute position upon modification of the charge distribution will surely depend on the specific properties of the solute, solvent, and confining framework. Approaches for predicting these location changes given the system characteristics and how to tune them—for example, in a catalyzed reaction to drive the reactants toward and products away from active sites on a pore surface—would be of significant utility in the rational design of mesoporous materials. This effort should be aided by the analysis of the underlying driving forces presented here, which can be used to draw inferences about a wider range of chemical processes that involve molecular charge rearrangement. In particular, it is notable that while entropic effects are important in determining the solute position distribution, the difference between the ground- and excited-state solute locations is dominated by the internal energy, i.e., competition between pore-solute and solvent-solute interactions, including H-bonding and non-specific solvation. These energetic effects are more easily predicted and controlled than the entropic ones. Combined with the fact that the free energy variations are comparatively small, ∼1 − 3 kcal/mol, these results suggest that relatively modest changes in the confining framework or solute interactions may have a significant impact on the distribution of solute positions.

Moving forward, a number of challenges remain. Predictive models, which do not require extensive simulations, for solute positions based on the solute charge distribution along with solvent and confining framework properties represent a key goal. In addition, an improved understanding of dynamical properties that elucidates the mass transport limitations and rate processes in mesoporous materials is still needed. Experimental efforts, including additional^17–19 complementary measurements that provide insight into solute location in nanoconfined solvents will, along with theoretical and simulation work, be critical to such efforts.

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