Distribution and Dynamic Properties of Xenon ... - ACS Publications

We have investigated the structural and dynamic properties of Xe dissolved in the ionic liquid crystal (ILC) phase of 1-hexadecyl-3-methylimidazolium ...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Distribution and Dynamic Properties of Xenon Dissolved in the Ionic Smectic Phase of [C16mim][NO3]: MD Simulation and Theoretical Model Diego Frezzato*,† and Giacomo Saielli*,‡ †

Department of Chemical Sciences, University of Padova, Via Marzolo 1, 35131 Padova, Italy CNR Institute on Membrane Technology, Unit of Padova, Via Marzolo 1, 35131 Padova, Italy

‡

ABSTRACT: We have investigated the structural and dynamic properties of Xe dissolved in the ionic liquid crystal (ILC) phase of 1-hexadecyl-3-methylimidazolium nitrate using classical molecular dynamics (MD) simulations. Xe is found to be preferentially dissolved within the hydrophobic environment of the alkyl chains rather than in the ionic layers of the smectic phase. The structural parameters and the estimated local diffusion coefficients concerning the short-time motion of Xe are used to parametrize a theoretical model based on the Smoluchowski equation for the macroscopic dynamics across the smectic layers, a feature which cannot be directly obtained from the relatively short MD simulations. This protocol represents an efficient combination of computational and theoretical tools to obtain information on slow processes concerning the permeability and diffusivity of the xenon in smectic ILCs.

1-hexadecyl-3-methylimidazolium nitrate, [C16MIm][NO3],18 we present a computational and theoretical study of the structural and dynamic properties of Xe dissolved in the ionic smectic phase of the same ILC. The qualitative picture and the “raw data” are obtained from MD simulations on the nanosecond time scale which provide information about the positional probability distribution of the xenon at thermal equilibrium, and about the local diffusion coefficients, which depend on the microviscosity sensed by such a probe. We anticipate that the profile of the positional distribution features equivalent maxima and minima in correspondence, respectively, with the centers of the aliphatic and ionic layers of the ILC smectic phase. This means that the xenon preferentially resides in the aliphatic layers and the “jump” among them is an activated process. On the theoretical side, first we validate the adoption of a diffusive model for the Xe dynamics, and then we elaborate the corresponding Smoluchowski equation,19 which describes the time evolution of nonequilibrium positional probability distribution of the Xe atom. All physical dynamic observables (e.g., characteristic times and time-correlation functions or related spectral densities assessable from simulations or experiments) are ultimately determined by such a probability distribution. Here we shall focus on the estimate of two dynamic features: (1) the rate of crossing of the ILC ionic layers and (2) the “macroscopic” diffusion coefficients which refer to the migration of xenon if tracked on the long time

1. INTRODUCTION Xenon is a very important probe for the investigation of the structural properties of various kinds of matrixes since it is relatively chemically inert but at the same time, being highly polarizable, its spectroscopic signatures are generally highly sensitive to the environment. For example, when adsorbed on surfaces, photoelectron spectroscopy of Xe reveals atomic details on the surface structure and defects.1 It has also been used as a probe of cavity sites in proteins by means of X-ray scattering experiments after treatment of the proteins by highpressure xenon.2,3 Among the many techniques which exploit xenon as a probe, 129Xe NMR is probably the most widely used,4 because of the strong dependence of the 129Xe chemical shift on the surrounding environment. Therefore, 129Xe NMR has been used to study membranes,5 zeolites,6−8 liquid crystals (LCs),9 and recently ionic liquids (ILs).10,11 Given the amount of experimental work that has been carried out, it is not surprising that a large number of computational studies, mostly based on molecular dynamics (MD) simulations, have addressed the problem of the structural and dynamic properties of Xe, e.g., in membranes,12,13 LCs,14 LCs confined into nanocavities,15 and ILs.16,17 In particular, in ref 16 we have analyzed the structural and dynamic properties of Xe and its cage in 1-butyl-3-methylimidazolium chloride and hexafluorophosphate, while in ref 17 we have rationalized, by combining MD simulations and relativistic DFT calculations, the chemical shift dependence of 129Xe on the counteranion in the same two ILs. In this work, continuing the previous studies and following a recent MD investigation of the structure and void distribution in the smectic phase of the ionic liquid crystal (ILC) © XXXX American Chemical Society

Received: December 21, 2015 Revised: February 4, 2016

A

DOI: 10.1021/acs.jpcb.5b12470 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B scale; these coefficients are those which would be achieved, for example, by means of pulsed field gradient (PFG) NMR measurements20 or by other techniques sensitive to long-range displacements. The physical ingredients required for the calculations are an estimate of the local diffusion coefficient of the xenon, and the profile of its positional probability distribution at thermal equilibrium; both these ingredients are provided by the MD simulations, as stated above. The estimated crossing rate will be qualitatively compared with the observed frequency of interlayer jumps observed in the simulated trajectories. On the contrary, we shall show that the macroscopic diffusion coefficients cannot be directly attained by the simulated trajectories on the nanosecond time scale; however, the calculations based on the theoretical modeling making use of the physical ingredients provided by the MD simulations allow one to bridge the gap. In this respect, the present study depicts a methodology to exploit relatively short MD simulations to parametrize selfconsistently a diffusive model aimed at computing, with likely accuracy, the macroscopic transport coefficients which control, ultimately, the permeability of the xenon in the smectic ILCs. The same methodology can be also applied, with minor modifications, to other kinds of solutes.

Figure 1. Top: snapshot of the simulation box at 500 K (N in blue, O in red, C in cyan, Xe in yellow). Bottom: structural formula of the ILC compound [C16MIm][NO3]. Visualization of the simulation box is done with the software VMD.26

2. COMPUTATIONAL AND THEORETICAL METHODS 2.1. Computational Details of the Simulations. A simulation box was built by starting from the results obtained in the previous simulation of pure [C16MIm][NO3].18 Therefore, from a 512 ion pairs box, we have removed one ion pair selected at random and added two xenon atoms at the positions previously occupied by the terminal methyl group, C16, and the anion, NO3−. Only two xenon atoms were included in the box to avoid a perturbation of the ILC structure by the solute and to avoid Xe−Xe interactions. Then the final box contains 33217 atoms. The integration time step was set to 0.5 fs, the cutoff for the Lennard-Jones interaction was set to 14 Å, and the electrostatic interactions were treated using the Ewald sum with automatic parameter optimization to have a precision of 10−6. Two independent boxes were equilibrated at 500 and 550 K, respectively, for 5 ns, after which a production run of 10 ns was started in the NPT ensemble at 1 atm with the DL_POLY Classic software.21 All remaining details of the simulations are the same as already described in ref 18. Thus, concerning the ILC, we have used the CL&AP force field,22−24 while the Xe force field paramerters are from Bohn et al.25 Configurations were saved every 2000 time steps, that is, with a 1 ps frequency, in a trajectory file for further analysis. As an example, a snapshot of a simulation box at 500 K is shown in Figure 1. In the snapshot of Figure 1, we note that the box is composed of two smectic layers on the xy plane; that is, the director of the phase is along the z axis. One of the two ionic layers, which are composed by the nitrate anions and the imidazolium heads, is located in the middle of the box, and the other one is located at the boundaries of the box along the z axis; that is, it is split by the periodic boundary conditions, half on the top of the box and the other half at the bottom of the box. In contrast, the two hydrophobic layers are located in the middle of the two halves of the box. The smectic periodicity, L, obtained from the simulation is 37.3 Å at 500 K and 37.8 Å at 550 K. A detailed description of the ILC structure of the pure [C16MIm][NO3], in the absence of Xe solute atoms, is reported in ref 18.

2.2. Diffusive Motion of Xenon: Modeling and Computational Details. The translational motion of a xenon atom in the ILC phase is here modeled as a diffusive dynamic process. This is licit if the time interval Δt between two subsequent observations of the atom’s location is long enough that inertial effects are negligible, namely, if the time self-correlation between the velocity components is lost.19 The fulfillment of such a condition requires that Δt ≫ m/ξ̅, where m is the mass of Xe and ξ̅ is the average viscous friction which is experienced in the anisotropic medium by such a moving probe. Under the assumption of being in the diffusive regime of motion, Einstein’s relation19 states that D̅ = kBT/ξ̅, where D̅ is the average local diffusion coefficient. By combining these expressions, one finds that it has to be Δt ≫ Δt* = (m/kBT)D̅ = (MWXe/RT)D̅ , where kB is the Boltzmann constant, R is the “gas constant”, T is the absolute temperature, and MWXe is the molar weight of the xenon (131.29 × 10−3 kg/mol). By taking 10−8 m2/s as a tentative upper limit for D̅ , it follows that Δt* falls on the subpicosecond scale (on the order of 0.3 ps). This means that if the motion is tracked, as is done in the present MD simulations, at intervals Δt of 1 ps, the translational dynamics can be likely well described as a diffusive motion. On the basis of these premises, the key quantity which encodes all information on the dynamic features is the nonequilibrium distribution p(z, t) such that p(z, t) δz is the probability of finding the Xe atom in the slice between z and z + δz at time t once its location at time 0 is provided. The time evolution of p(z, t) is described by the diffusion equation in the Smoluchowski form:19 ∂p(z , t ) ∂J(z , t ) =− ∂z ∂t

(1)

where J(z, t) is the probability flux expressed by J (z , t ) = − D

,loc (z)

ρ (z )

∂[p(z , t )/ρ(z)] ∂z

(2)

In this equation, ρ(z) is the stationary profile which is reached by p(z, t) in the limit t → ∞, that is, the Boltzmann B

DOI: 10.1021/acs.jpcb.5b12470 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

By monitoring the diffusion of probes on a long time scale (hence in the long range), one extracts “macroscopic” diffusion coefficients which are related to the “local” ones and to the features of the spatial modulation of the effective potential. For smectic-like phases with periodicity L, and assuming constant local diffusion coefficients (i.e., by neglecting the spatial modulation of the microviscosity), the theory developed in ref 28 yields D ,macro = D ,loc ϵpot (5)

equilibrium distribution at thermal equilibrium. The actual form of ρ(z) is here constructed by the MD simulations and is shown in the following section. For convenience, by inverting the Boltzmann expression ρ(z) ∝ e−V(z)/kBT, such a distribution is specified by an effective potential V(z) via V (z) = const − kBT ln[ρ(z)]

(3)

where “const” stands for an immaterial additive constant. As will be shown, ρ(z) has equivalent maxima at the centers of the aliphatic layers and minima in correspondence with the ionic parts; conversely, the profile of V(z) displays equivalent “wells” which fall in the aliphatic regions and which are separated by energy barriers at the middle of the ionic regions. The other ingredient in eq 2, D∥,loc(z), is the local diffusion coefficient of the xenon in the direction parallel to the ILC director (i.e., collinear with the layers’ normal). Then the transverse local diffusion coefficients on the xy plane will be denoted as D⊥,loc(z). Note that the diffusion coefficients are taken to be z-dependent in all generality. Such a possible dependence on z is meant to arise only from a modulation of the microviscosity probed in the different regions of the ILC. In practice, for the coefficient entering eq 2, we shall adopt the form D∥,loc(z) = D0 f mod(z), with D0 constant and f mod(z) a dimensionless positive-valued function. The diffusion equation is here solved numerically by means of the finite-differences scheme described in the Appendix. The unbounded ILC phase is handled by considering a range extended from −kL to +kL, with k an integer, and employing periodic boundary conditions compatible with the translational periodicity of the phase. For a chosen k, the discretization route yields a relaxation matrix (see the “symmetrized form” given in the Appendix) whose eigenvectors specify the “relaxation modes”, and the eigenvalues are the rates of the dynamical processes which take place in the selected spatial range. In particular, the eigenvalues are split into two sets separated by a gap. The set of smaller eigenvalues corresponds to slow activated transitions among wells of the effective potential, while the set of much higher eigenvalues is related to the fast fluctuations inside the wells. The number of low eigenvalues depends on the extension of the considered spatial range (i.e., on the number of wells included in it); however, the highest eigenvalue of such a set, which is denoted as λcross in the following, is fixed and unambiguously corresponds to the transition between two adjacent wells (i.e., the fastest process among the slow ones). As the energy barrier increases, λcross decreases roughly exponentially by following an Arrhenius-like law. Moreover, the method of the “localized site functions” developed by Moro and Nordio27 allows one to establish a link between λcross and the kinetic rate kcross of the unidirectional process of barrier crossing. For equivalent minima of potential, the outcome is kcross =

λcross 2

where ϵpot < 1 is the dimensionless factor ⎛ ϵpot = L2⎜ ⎝

∫0

L

⎞−1⎛ dz e−V (z)/ kBT ⎟ ⎜ ⎠ ⎝

while D⊥,macro ≡ D⊥,loc

∫0

L

⎞−1 dz e+V (z)/ kBT ⎟ ⎠

(6)

(7)

Equation 5 expresses the slowing of the translational motion along the phase director due to the energy barriers that exert a hindering effect (ϵpot decreases, roughly exponentially, as the energy barrier increases28), while eq 7 states that the macroscopic and local diffusion coefficients coincide on the transverse directions over which the diffusion is free.

3. RESULTS AND DISCUSSION 3.1. Structural Information from the MD Simulations. Structural information concerning xenon dissolved in the ILC matrix can be obtained by inspection of the radial distribution functions (RDF) between Xe and representative atoms of the ILC ion pair. These are reported in Figure 2. We note that

Figure 2. Radial distribution function: top, 550 K; bottom, 500 K.

there is a significant correlation between Xe and the carbon atoms of the alkyl chains, especially the terminal methyl group, C16. The first peak of the RDF decreases as we move along the alkyl chain, indicating a weaker correlation as we approach the charged part of the imidazolium cation. In fact, no correlation is observed between Xe and either the imidazolium ring or the anion. These results are in perfect agreement with those reported recently by some of us17 and by Morgado et al.10 concerning Xe dissolved in imidazolium-based ionic liquids. There, it was observed that, even in the isotropic phase of ILs, the RDF of Xe with the terminal methyl of the alkyl chain was the most structured with a clear peak, while the intensity of the first peak

(4)

The inverse of this kinetic rate is then taken as a characteristic time, τcross, of transition of the xenon atom across the ionic layer. In summary, the input information required in eqs 1 and 2 is constituted by (i) the profile of the effective potential V(z), which specifies the equilibrium distribution by means of eq 3, and (ii) the profile of D∥,loc(z), which is an estimate of D0 and a likely form of f mod(z). These issues are discussed in the following section. C

DOI: 10.1021/acs.jpcb.5b12470 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the ionic layer presents a close packing of imidazolium heads and anions so that no free volume exists. Although the solubility of Xe is due to several factors, the availability of voids in the hydrophobic layer indeed increases the probability of finding Xe there rather than in the ionic layers. From eq 3, we can now calculate the effective potential in kBT units (apart from the immaterial additive constant) as −[ln ρ(z)]. The profile is shown in Figure 4 and is revealed to

decreased for the RDF of Xe with carbons of the alkyl chain closer to the imidazolium ring. This suggested that, in the isotropic phase of ILs, Xe was preferentially located in the hydrophobic domains of the nanosegregated structure. It is not surprising, therefore, that such an effect is strongly enhanced in ILCs where the nanosegregation extends to the macroscopic scale, at least in two dimensions, with a clear separation and alternation of ionic and hydrophobic layers. This result alone already suggests that Xe is strongly confined within the hydrophobic layer of the ILC, as will be apparent in the following analysis of the density distribution resolved parallel to the director of the phase. In Figure 3 we show the probability distribution ρ(z) of Xe resolved along the director of the phase. The distribution has

Figure 4. Profiles of the effective potential: top, 550 K; bottom, 500 K. Red lines connect the raw data from the equilibrium distribution constructed by histograms. The black lines are the fitting curves (see the text for details).

Figure 3. Probability distribution of Xe along the director of the phase: top, 550 K; bottom, 500 K.

been obtained by means of a histogram construction from the MD data. The z coordinate is expressed in units of the ILC periodicity length L (the simulation box extends from −1 to +1 in such units). As it appears from the snapshot in Figure 1, one of the ionic layers is at the edges of the box along the z axis (because of the periodic boundary conditions), while the other ionic layer is in the middle of the box, therefore around z = 0. As shown in Figure 3, finding a Xe atom in the middle of the hydrophobic layers has the maximum probability, while the probability of finding it in the ionic layers is strongly reduced. There is an evident lack of symmetry in the distribution at 550 K, with one layer much more populated than the other. This is due to the limited length of the simulation (10 ns production run) and to the limited number of Xe atoms in the box. As a consequence, the two Xe atoms may not spend the same amount of time in the two layers (see below the discussion of the Xe trajectories). This is, however, not a major problem since we can symmetrize the probability distribution referring it to just one smectic layer. To do this, we have averaged the profile between 0.5 and 1.0 with that between −0.5 and 0.0, and the profile between −1.0 and −0.5 with that between 0.0 and 0.5. After this procedure, we obtain a more symmetric curve. It is remarkable that the distribution of Xe within the layers appears to be strongly correlated with the distribution of voids calculated in the pure ILC in ref 18. There, it was found that voids of sufficient size to host a probe such as a xenon atom are almost exclusively formed within the hydrophobic layer, while

be broader and lower at higher temperature. The produced values have then been fitted with the lowest number of Fourier components which yield a satisfactory reproduction of the visible features. The functional form for both cases, 500 and 550 K, is 4

Vfit(z)/kBT = p0 +

∑ pk cos(2πkz/L)

(8)

k=1

where only cosines are employed due to the pairwise symmetry about z = 0. The set of best parameters is given in Table 1, and the fitting curves are also shown in Figure 4. Table 1. Fitting Parameters of the Effective Potential 550 K 500 K

p0

p1

p2

p3

p4

0.848 1.085

1.011 1.524

0.338 0.721

0.137 0.239

−0.008 −0.017

3.2. Simulated MD Trajectories of Xenon. In Figure 5, we show the trajectories (after removing the effect of the periodic boundary conditions) of the two xenon atoms during the 10 ns production at 550 and 500 K. Xenon explores a relatively large space, within the smectic layer, as indicated by the variation of the x and y coordinates. In contrast, the dynamics along the director, that is, the z coordinate, is characterized by a sort of small-step motion with interlayer D

DOI: 10.1021/acs.jpcb.5b12470 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

the xenon atom explores different local viscous environments in the course of the simulated 10 ns. We stress that eq 7 is applicable if the time τ is long enough so that the motion appears to be diffusive, but short enough so that the diffusion is locally “free” (i.e., if the effective potential is approximatively flat on the short length scale explored by the motion). Since the time step Δt of 1 ps is likely long enough to adopt the diffusive model, the coefficients Deff u,loc have been estimated by taking the smallest number of points which allow a good linear fit of Δu (τ )2 /2 vs τ. Namely, only five points have been considered. The time dependence of Δu = x, y, z (τ )2 /2 at 550 and 500 K is shown in Figure 6. On such a time scale (see the insets in the

Figure 5. Time evolution of the three Cartesian coordinates of the two Xe atoms at (top) 550 K and (bottom) 500 K. PBC effects have been removed. Dotted lines represent the values on the z axis corresponding to the middle of a hydrophobic layer.

jumps. For example, in Figure 5a, we note that, after more or less 5 ns of fluctuations roughly around z = −20 Å (that is, within one hydrophobic layer), the Xe atom jumps into the other hydrophobic layer, roughly centered about z = 20 Å, crossing the middle ionic layer (see also Figure 1 for the location of the layers along the z axis). In contrast, the second Xe atom, in Figure 5b, remains in the same layer for the whole length of the production run. At 500 K, Figures 5c,d, we note one Xe atom that crosses the middle ionic layer after about 3.5 ns, while the other atom crosses the other ionic layer at the bottom of the box after about 7.5 ns. The very limited statistics available does allow one only to say that the order of magnitude of the frequency of jumps across a ionic layer is about one every 10 ns. These two different kinds of dynamics, that is, the slow crossing of the ionic layers and the fast fluctuations within the potential wells, are clearly the result of the presence of the potential energy barrier displayed in Figure 4. As stated in section 2.2, the time scale of the former process is determined by the eigenvalue λcross of the relaxation matrix, while the time scales of the latter processes are determined by the higher eigenvalues. 3.3. Local and Macroscopic Diffusion Coefficients of Xenon from the MD Simulations. The local diffusion coefficients have been determined from the mean square displacement of the Xe atom.19 Namely, by denoting with u the unit vector which specifies a direction of monitoring (either the parallel direction or any perpendicular direction), the relation here exploited is Δu (τ )2 ≡

1 t tot

= 2Dueff,locτ

∫0

t tot

Figure 6. Time dependence of mean square displacements for the Xe atom at (top) 550 K and (bottom) 500 K along selected directions of monitoring. The insets show blowups of the same profiles on the picosecond time scale, within which the linear regression yields an estimate of the local diffusion coefficients.

figure), the three profiles are indeed linear. The resulting values of the diffusion coefficients at 500 and 550 K are given in Table 2, where D⊥,loc is the average of the outcomes from the displacements on the x and y directions, while D∥,loc refers to the z direction. Note that the values for the transverse and longitudinal diffusion coefficients are very close. This means that the

dt0 [(rXe(t0 + τ ) − rXe(t0)) ·u]2 (9)

where the ensemble average has been computed, as indicated in the equation above, via a time average numerically performed by “sliding” the initial time over the entire trajectory at disposal.29 The accuracy is further improved by averaging the outcomes for the two Xe atoms simultaneously present in the simulation box. Correspondingly, the resulting diffusion coefficients have to be interpreted as “effective” ones, since

Table 2. Local Diffusion Coefficients (10−9 m2 s−1)

550 K 500 K E

D⊥,loc

D∥,loc

D̅

9.3 ± 0.2 7.1 ± 0.2

10.0 ± 0.2 8.3 ± 0.2

9.5 7.5

DOI: 10.1021/acs.jpcb.5b12470 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

jumps, by (1) slightly increasing the energy barrier and/or (2) introducing a modulation of the local diffusion coefficient such that D∥,loc is smaller in the ionic region and higher in the aliphatic region. The effect of (1) can be quantified by considering that the barrier crossing is an activated process and the eigenvalue λcross roughly follows an Arrhrenius-like law for large enough barriers (Kramers’ theory19). An increase of the energy barrier of only 1 kBT causes τcross to roughly increase by a factor 2.7. The effect of (2) has been inspected by performing calculations of λcross with a function f mod(z) whose profile is taken to be the same as that of the effective potential. This is nothing but a subjective choice devoid of a physical basis and adopted only for illustrative calculations. Namely, we have employed the form f mod(z) = 1 + 2f [w(z) − 1/2], with w(z) = [V(0) − V(z)]/[V(0) − V(L/2)] (we recall that z = 0 corresponds to the maximum of the potential in the ionic region, while z = L/2 corresponds to a minimum in the aliphatic region) and where f sets the amplitude of the modulation: (D∥,ionic − D0)/D0 = −f and (D∥,aliphatic − D0)/D0 = +f. The calculations with f = 0.5 show that τcross increases by a factor 1.6 at 500 K and by a factor 1.5 at 550 K. Such a modulation hence plays a minor role compared to the effect of an increase of the energy barrier. In conclusion, the slightly short crossing times may be likely attributed to an underestimation of the energy barrier due to the loss of accuracy in V(z) in the ionic region where the statistical sampling is poorer. Now let us turn to the estimate of the macroscopic diffusion coefficient along the ILC director. Application of eqs 5−7 yields the results reported in Table 3 (columns 5 and 6). The effect of the energy barriers causes the diffusion along the director to be hindered (as naturally expected), and the hindrance to increase as the temperature is lowered since the barrier crossing is an activated process. Again, such a theoretical estimate is sensitive to a likely underestimation of the energy barrier. The algebraic elaboration of eq 6 under a parabolic approximation of the effective potential at the maxima and minima indicates that D∥,macro decreases roughly exponentially as the energy barrier increases.28 By following the same reasoning made above about the uncertainty on τcross, an increase of only 1 kBT causes D∥,macro to decrease by a factor of 2.7.

anisotropy of the viscous environment is small or, in any case, that it is not detected by this spherical probe. On this basis, for the calculations we shall likely take the local diffusion as isotropic from the point of view of the viscous drag; that is, we shall employ D∥,loc = D̅ , where D̅ stands for the average of the values. It is interesting to note that D̅ turns out to be very close to the diffusion coefficient in the liquid alkane n-C16H34 estimated by employing an empiric law specific for the diffusion of xenon in linear alkanes.30 This outcome supports the physical picture that the leading contribution to the effective diffusion coefficient is due to the viscous drag in the aliphatic regions (i.e., in the wells of potential). Moreover, also on a time scale of about 1 ns (after which the statistics drops dramatically), there is a linear relation between mean square displacements and time. On the other hand, it not possible to extract, from the linear regressions of these profiles, a likely estimate of the macroscopic diffusion coefficient. In fact, we see that for the x and y directions of monitoring the profiles are close to one another (as it must be because of the statistical equivalence of these directions) only at 550 K. We conclude that the statistics is too poor to make a quantitative analysis, even if the average slope of the fitted profiles at 550 K yields a value of D⊥,macro of (4.5 ± 0.2) × 10−9 m2 s−1 which is on the same order of magnitude as D⊥,loc in Table 2, in accord with eq 7. The situation is even worse for the macroscopic diffusion parallel to the ILC director. In fact, on a nanosecond time scale, the Xe atom typically remains confined within one to two vicinal layers (see the trajectories in Figure 5); hence, the character of the translational motion is still “local”, and a longrange diffusion coefficient cannot be extracted. With these remarks, we stress that the macroscopic diffusion coefficients in the ILC phase cannot be obtained even from such relatively long MD simulations. On the other hand, the sound estimate of the local diffusion coefficients, in addition to the construction of the equilibrium distribution ρ(z), both obtained from the MD simulations, allows one to calculate the macroscopic coefficients as shown in the following section. 3.4. Theoretical Estimates of the Crossing Time and of the Macroscopic Diffusion Coefficient. A FORTRAN code has been written to implement the strategy outlined in section 2.2 with the parametrization presented in sections 3.1 (effective potential) and 3.2 (local diffusion coefficient). The details of the computational setup are given in the Appendix. The reference situation is that of a constant local diffusion coefficient; namely, we set D∥,loc(z) = D0 f mod(z), with D0 = D̅ , where D̅ is reported in Table 2, and f mod(z) = 1. The outcomes of the calculations are collected in Table 3 (see columns 1−4). First, note that the characteristic times of the barrier crossing fall in the correct order of magnitude if they are compared with the number of “jumps” which happened in the course of the simulated 10 ns (see Figure 5). Notably, the crossing times can be slightly increased, to better match the observed frequency of

4. CONCLUSIONS In this work we have shown how relatively short molecular dynamics simulations (here performed on a time scale of 10 ns) can provide information about both the energetics (positional distribution) and transport coefficients (local diffusion coefficients) of xenon dissolved in the layered smectic structure of an ILC. An important outcome of the MD simulations is that Xe is mostly dissolved in the hydrophobic layers, while its efficiency in “probing” the ionic layer is much lower. This observation is important in view of the extensive use of Xe as a probe of inhomogeneous environments, such as ionic liquids, which, at the nanoscale level, are also nanosegregated in hydrophobic and ionic regions. It is also noteworthy that the xenon distribution within the ILC phase is qualitatively correlated with the void distribution found in ref 18, both reaching a maximum in the middle of the hydrophobic layer. The information obtained from the MD simulations is then used to parametrize the stochastic description of the xenon translation motion. The likelihood of the parametrization is assessed by numerically solving the diffusion equation and checking that the computed rate of crossing of the ionic layers (which separate two vicinal wells of potential) is in agreement

Table 3. Dynamic Properties from the Theoretical Model barrier crossing process

550 K 500 K

long-range diffusion

λcross (L2/D∥,loc)

kcross (s−1)

τcross (ns)

εpot

D∥,macro (10−9 m2 s−1)

2.91 1.18

9.7 × 108 3.2 × 108

1.0 3.1

0.58 0.27

5.4 2.0 F

DOI: 10.1021/acs.jpcb.5b12470 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B with the frequency of “jumps” directly observed in the MD simulations. The local diffusion coefficients and the energetics, both obtained from the MD simulations, are then used to compute the macroscopic diffusion coefficients according to eqs 5−7. Such a computational/theoretical self-consistent approach hence constitutes a tool to achieve the macroscopic transport coefficients that would be otherwise hardly obtainable by MD simulations, since the motion of the probe is hindered by the layer-to-layer energy barrier (to be overtaken by thermal fluctuations) and a longer simulation time than the 10 ns used here would be necessary to let the probe explore a sufficiently wide spatial scale. For example, by employing the value of D∥,macro in Table 3, one sees (using the Einstein relation for the root-mean-square displacement) that a spatial range even of only 10 times the periodicity L would require, to be covered at 500 K, an MD simulation of about 350 ns, that is, 35 times longer than the production run used here. In our analysis, we have pointed out that the most critical step is the accurate estimate of the energy barrier to be overtaken to move from a well of potential to a vicinal one. In fact, D⊥,macro roughly decreases exponentially as such a barrier increases. On the other hand, the comparison between computed and observed crossing rates (which also depend in a similar way on the magnitude of the energy barrier) provides a self-consistent check of reliability on the outcomes. We stress that the same methodology can be easily adapted and applied to other probes different from xenon. The approach here depicted to obtain D⊥,macro and D∥,macro can be adopted to predict the diffusivity in a novel ILC. Moreover, we recall that the macroscopic diffusion coefficients can be experimentally obtained, for example, by means of PFG-NMR measurements. Thus, the present approach can be interfaced with measurements to assess, by a comparison of the outcomes, the likelihood of the physical information at the microscopic level.

MN ,1 = − ρ(z1)−1 ρ(zN+) D

MN , N − 1 = − ρ(zN − 1)−1 ρ(zN−) D MN , N = ρ(zN )−1[ ρ(zN+) D

M1,2 = − ρ(z 2)−1 ρ(z1+) D

Mn , m = δn , mρ(zn)−1[ρ(zn+) D

−1

M1, N = − ρ(zN )

ρ(z1−)

D

(Δz)−2

− ,loc(z1 )

−2

+ ρ(zN−) D

(Δz)−2 − −2 ,loc(zN )](Δz)

+ ,loc(zn )

+ ρ(zn−) D

− −2 ,loc(zn )](Δz)

− δm , n + 1ρ(zn + 1)−1 ρ(zn+) D

+ ,loc(zn )

(Δz)−2

− δm , n − 1ρ(zn − 1)−1 ρ(zn−) D

− ,loc(zn )

(Δz)−2 (A3)

In the expressions above, δ denotes the Kronecker delta function and zn± = zn ± Δz/2, with z−1 = zmin and z+N = zmax. The eigenvalues of M are the rates of relaxation of p(z, t) toward the stationary profile ρ(z), while the eigenvectors can be interpreted as the related “modes”. One eigenvalue is equal to 0, and the related eigenvector corresponds to (the discretized form of) ρ(z). The other eigenvalues are all realvalued and non-negative. For computational convenience, by turning to p̃n(t) = . ̃ ̃ , where the p = −Mp p(zn, t)/ρ(z)1/2, one achieves the format ∼ new matrix M̃ , with elements M̃ n,m = Mn,m[ρ(zn)/ρ(zm)]1/2, is symmetric and has the same eigenvalues as M. Calculations have been performed on scaled quantities by setting the length unit equal to the ILC periodicity L and the time unit equal to L2/D0. The diagonalization of the matrix M̃ has been performed with Jacobi’s pivoting route. The results presented in Table 3 have been obtained by performing accurate calculations with k = 1, that is, by taking into account only two vicinal wells of potential. Accordingly, only fluctuations inside the wells and the barrier crossing event are described. In this case, a single eigenvalue of small magnitude, directly identified as λcross, appears to be well separated by the higher ones. The results presented refer to N = 300 discretization intervals, after the check that passing from 200 to 300 intervals affects the value of λcross by at most 3.5 × 10−3%. We have checked that calculations with any k > 1 yield the same value of λcross identified as the largest eigenvalue of the “slow set” (see section 2.2). Moreover, we have also performed calculations with very extended spatial ranges (namely, up to k = 100) by applying reflecting boundary conditions at the extrema, i.e., J(zmin, t) = J(zmax, t) = 0. This mimics the situation of diffusion into an extended slab of ILC phase. The outcome is that, as k increases, the largest eigenvalue of the “slow set” asymptotically tends to the λcross determined by applying periodic boundary conditions.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.

− −2 ,loc(z1 )](Δz)

+ ,loc(z1 )

− ,loc(zN )

and the elements for rows 2 ≤ n ≤ N − 1 are

APPENDIX The diffusion equation is solved numerically by means of a finite-differences scheme based on the approximation of the spatial derivatives by incremental ratios with a partition of the z axis into N intervals of equal width Δz. The spatial range is limited within zmin = −kL and zmax = +kL, where L is the spatial periodicity of the ILC and k ≥ 1 is an integer value. The coordinate z = 0 is set to correspond to a maximum of the effective potential, i.e., to the middle of an ionic layer which separates two vicinal hydrophobic layers (potential wells). The constraint J(zmin, t) = J(zmax, t), that is, persistent equality of the probability flux at the extrema of the spatial range, corresponds to the application of periodic boundary conditions. Such a choice has been pursued to handle the unboundedness of the ILC phase (see the remarks below). The discretization procedure turns eq 1 into the matrix format ṗ = −Mp, where the column vector p has elements pn(t) = p(zn, t), with zn the central point of the nth interval and M the N × N matrix whose first-row elements are + − ,loc(z1 ) + ρ(z1 ) D

+ ,loc(zN )

(Δz)−2

(A2)

■

M1,1 = ρ(z1)−1[ρ(z1+) D

+ ,loc(zN )

■

ACKNOWLEDGMENTS Financial support from MIUR (Grants PRIN 2010N3T9M4 and FIRB RBAP11C58Y) is gratefully acknowledged. We thank CNR-CAS bilateral agreement 2014-2016 for support. We thank the LICC computational facility of the Department of

(Δz)

(A1)

while the last-row elements are G

DOI: 10.1021/acs.jpcb.5b12470 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(16) Frezzato, D.; Bagno, A.; Castiglione, F.; Mele, A.; Saielli, G. MD Simulation of Xenon in Ionic Liquids: Disentangling the Cationic and Anionic Cage Effects on the Structural and Dynamic Properties. J. Mol. Liq. 2015, 210 (Part B), 272−278 10.1016/j.molliq.2015.03.039. (17) Saielli, G.; Bagno, A.; Castiglione, F.; Simonutti, R.; Mauri, M.; Mele, A. Understanding Cage Effects in Imidazolium Ionic Liquids by 129 Xe NMR: MD Simulations and Relativistic DFT Calculations. J. Phys. Chem. B 2014, 118, 13963−13968. (18) Saielli, G. Fully-Atomistic Simulations of the Ionic Liquid Crystal [C16MIm][NO3]: Orientational Order Parameters and Voids Distribution. J. Phys. Chem. B 2016, DOI: 10.1021/acs.jpcb.5b12469. (19) Zwanzig, R. Nonequilibrium Statistical Mechanics; Oxford University Press: New York, 2001. (20) Johnson, C. S., Jr. Diffusion Ordered Nuclear Magnetic Resonance Spectroscopy: Principles and Applications. Prog. Nucl. Magn. Reson. Spectrosc. 1999, 34, 203−256. (21) Smith, W.; Forester, T. R.; Todorov, I. T. DL_POLY Classic, version 2010; Daresbury Laboratory: Daresbury, U.K. http://www. ccp5.ac.uk/DL_POLY_CLASSIC (accessed Feb 5, 2016). (22) Canongia Lopes, J. N.; Deschamps, J.; Padua, A. A. H. Modeling Ionic Liquids using a Systematic all-Atom Force Field. J. Phys. Chem. B 2004, 108, 2038−2047. (23) Canongia Lopes, J. N.; Deschamps, J.; Padua, A. A. H. Modeling Ionic Liquids using a Systematic all-Atom Force Field. J. Phys. Chem. B 2004, 108, 11250−11250. (24) Canongia Lopes, J. N.; Padua, A. A. H. Molecular Force Field for Ionic Liquids III: Imidazolium, Pyridinium, and Phosphonium Cations; Chloride, Bromide, and Dicyanamide Anions. J. Phys. Chem. B 2006, 110, 19586−19592. (25) Bohn, M.; Fischer, J.; Kohler, F. Prediction of Excess Properties for Liquid Mixtures: Results from Perturbation Theory for Mixtures with Linear Molecules. Fluid Phase Equilib. 1986, 31, 233−252. (26) Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (27) Moro, G.; Nordio, P. L. Diffusive and Jump Description of Hindered Motions. Mol. Phys. 1985, 56, 255−269. (28) Moro, G.; Nordio, P. L.; Segre, U. Translational Diffusion in Smectic-A Phases. Chem. Phys. Lett. 1984, 105, 440−443. (29) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, U.K., 1991. (30) A parametric law which relates the diffusion coefficient D of xenon to the viscosity η of alkanes has been proposed by Pollack et. al.: Pollack, G. L.; Kennan, R. P.; Himm, J. F.; Stump, D. R. Diffusion of Xenon in Liquid Alkanes: Temperature Dependence Measurements with a New Method. Stokes-Einstein and Hard Sphere Theories. J. Chem. Phys. 1990, 92, 625−630. The relation is Dηp = AT, where D is expressed in square centimeters per second, η is expressed in centipoise, T is the temperature in kelvin, p = 0.708, and A = 9.80 × 10−8. As a model of the aliphatic region of the ILC, we have taken the liquid n-C16H34. The viscosity has been calculated by applying the following parametric form given by Yaws (Yaws, C. L. Handbook of Viscosity; Gulp Publishing Co.: Houston, TX, 1995; Vol. 3): log η = A + B/T + CT + DT2, with A = −8.1894, B = 1.5571 × 103, C = 1.5270 × 10−2, and D = −1.2371 × 10−5. This relation is applicable in the temperature range between 291 and 721 K. The results are η(500 K) = 0.293 cP and η(550 K) = 0.199 cP, which lead to the values D(500 K) = 1.17 × 10−8 m2/s and D(550 K) = 1.69 × 10−8 m2/s, very close to those given in Table 2.

Chemistry of the University of Padova and the CINECA consortium (Bologna, Italy) for granting CPU time through ISCRA Projects HP10BDI66L and HP10C1V9ME.

■

ABBREVIATIONS IL, ionic liquid; ILC, ionic liquid crystal; LC, liquid crystal; NMR, nuclear magnetic resonance; MD, molecular dynamics; DFT, density functional theory; RDF, radial distribution function; PBC, periodic boundary conditions; PFG, pulsed field gradient

■

REFERENCES

(1) Guo, H.; Zaera, F. Xenon as a Probe for Minority Sites on Solid Surfaces. Nat. Mater. 2006, 5, 489−493. (2) Duff, A. P.; Trambaiolo, D. M.; Cohen, A. E.; Ellis, P. J.; Juda, G. A.; Shepard, E. M.; Langley, D. B.; Dooley, D. M.; Freeman, H. C.; Guss, J. M. Using Xenon as a Probe for Dioxygen-Binding Sites in Copper Amine Oxidases. J. Mol. Biol. 2004, 344, 599−607. (3) de Sanctis, D.; Dewilde, S.; Pesce, A.; Moens, L.; Ascenzi, P.; Hankeln, T.; Burmester, T.; Bolognesi, M. Mapping Protein Matrix Cavities in Human Cytoglobin through Xe Atom Binding. Biochem. Biophys. Res. Commun. 2004, 316, 1217−1221. (4) Fraissard, J. Xenon as a Probe Atom: Introduction, Characteristics, Investigation of Microporous Solids. Hyperpolarized Xenon-129 Magnetic Resonance: Concepts, Production, Techniques and Applications; The Royal Society of Chemistry: Cambridge, U.K., 2015; Chapter 1, pp 1−15. (5) Amor, N.; Hamilton, K.; Küppers, M.; Steinseifer, U.; Appelt, S.; Blümich, B.; Schmitz-Rode, T. NMR and MRI of Blood-Dissolved Hyperpolarized Xe-129 in Different Hollow-Fiber Membranes. ChemPhysChem 2011, 12, 2941−2947. (6) Huang, S.; Zhao, Q.; Chen, W.-H.; Han, X.; Bao, X.; Lo, P.-S.; Lee, H.-K.; Liu, S.-B. Structure and Acidity of Mo/H-MCM-22 Catalysts Studied by NMR Spectroscopy. Catal. Today 2004, 97, 25− 34. (7) Liu, S. B.; Lee, C. S.; Shiu, P. F.; Fung, B. M. Interpretaion of Xenon Adsorption Isotherms and Xe-129 NMR Chemical Shifts on Ion-Exchanged Nay Zeolites. Stud. Surf. Sci. Catal. 1994, 83, 233−241. (8) Smith, M. L.; Corbin, D. R.; Abrams, L.; Dybowski, C. Flexibility of Zeolite Rho: Xe-129 NMR Studies of Entrapped Xenon in Cd-Rho. J. Phys. Chem. 1993, 97, 7793−7795. (9) Lounila, J.; Muenster, O.; Jokisaari, J.; Diehl, P. Temperature Dependence of Nuclear Shielding and Quadrupolar Coupling of Noble Gases in Liquid Crystals. J. Chem. Phys. 1992, 97, 8977−8985. (10) Morgado, P.; Shimizu, K.; Esperanca, J. M. S. S.; Reis, P. M.; Rebelo, L. P. N.; Canongia Lopes, J. N.; Filipe, E. J. M. Using 129Xe NMR to Probe the Structure of Ionic Liquids. J. Phys. Chem. Lett. 2013, 4, 2758−2762. (11) Castiglione, F.; Simonutti, R.; Mauri, M.; Mele, A. Cage-Like Local Structure of Ionic Liquids Revealed by a 129Xe Chemical Shift. J. Phys. Chem. Lett. 2013, 4, 1608−1612. (12) Yamamoto, E.; Akimoto, T.; Shimizu, H.; Hirano, Y.; Yasui, M.; Yasuoka, K. Diffusive Nature of Xenon Anesthetic Changes Properties of a Lipid Bilayer: Molecular Dynamics Simulations. J. Phys. Chem. B 2012, 116, 8989−8995. (13) Jansen, J. C.; Macchione, M.; Tocci, E.; De Lorenzo, L.; Yampolskii, Y. P.; Sanfirova, O.; Shantarovich, V. P.; Heuchel, M.; Hofmann, D.; Drioli, E. Comparative Study of Different Probing Techniques for the Analysis of the Free Volume Distribution in Amorphous Glassy Perfluoropolymers. Macromolecules 2009, 42, 7589−7604. (14) Lintuvuori, J.; Straka, M.; Vaara, J. Nuclear Magnetic Resonance Parameters of Atomic Xenon Dissolved in Gay-Berne Model Liquid Crystal. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 2007, 75, 031707. (15) Karjalainen, J.; Vaara, J.; Straka, M.; Lantto, P. Xenon NMR of Liquid Crystals Confined to Cylindrical Nanocavities: A Simulation Study. Phys. Chem. Chem. Phys. 2015, 17, 7158−7171. H

DOI: 10.1021/acs.jpcb.5b12470 J. Phys. Chem. B XXXX, XXX, XXX−XXX