Molecular Dynamics Simulation of Highly Confined Glassy Ionic

Dec 29, 2015 - News Ed., Am. Chem. ... Chemical Engineering Program, Education City, Texas A&M .... were examined, focusing mainly on the detailed ana...
0 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCC

Molecular Dynamics Simulation of Highly Confined Glassy Ionic Liquids Georgios Kritikos,*,† Niki Vergadou,*,† and Ioannis G. Economou†,‡ †

Institute of Nanoscience and Nanotechnology, Molecular Thermodynamics and Modelling of Materials Laboratory, National Center for Scientific Research “Demokritos”, GR-153 10, Aghia Paraskevi Attikis, Greece ‡ Chemical Engineering Program, Education City, Texas A&M University at Qatar, P.O. Box 23874, Doha, Qatar S Supporting Information *

ABSTRACT: We present a molecular dynamics simulation study of 1octyl-3-methylimidazolium tricyanomethanide ([Omim+][TCM−]) ionic liquid capped by two silica planar surfaces. The study extends over a wide temperature range and various interwall distances. Our results indicate that the structure and dynamics of the confined system is significantly affected by the width of the film. At the shortest interwall distance of 25 Å, which is comparable to the ion pair dimensions, the bulk structure is breached. The dynamics of the cation in the adsorbed layer is accelerated for the time scale of 1 ns and decelerates for longer time scales. In the most confined film, we observe a suppression of the cooperative characteristics in the diffusion. The whole phenomenon seems to be related to an Arrhenius behavior. Our proposed model suggests a stable, static liquid path in the center of the pore that facilitates the diffusion. The simulations results are consistent with a recent experimental study on the same confined system.

1. INTRODUCTION Ionic Liquids (ILs) are salts with a melting point (Tm) or a glass transition temperature (Tg) below 100 °C.1 In case they retain liquid characteristics like diffusion at room temperature, they are called room temperature ionic liquids (RTILs).1 There is a variety of cation and anion pairs that can form an IL. Generally, one of the counterions is organic and charge delocalization is present in the IL systems, preventing the formation of a crystal state and therefore preserving a wide liquid range. In many cases, RTILs exhibit dynamics of complex systems (cooperative diffusion) and at a lower temperature (Tg) they undergo glass transition (glassy materials). Among their common properties are the very low vapor pressure, the high solubility for a wide range of compounds, and their thermal stability.1 One of the most studied families of RTILs is the imidazolium-based ILs. For a specific anion, the structure of the liquid depends on the length of the alkyl tail of the cation. Systems with cations of long alkyl chains can form nanoscale polar and nonpolar domains.2 X-ray diffraction3,4 confirms the existence of nanosize heterogeneities that depend on the length of the alkyl chain. Longer alkyl groups in most cases lead to increased CO2 absorption capacity due to an enhancement in the free volume.5 In the case of confinement of ILs into pore channels, the produced membranes are called supported ionic liquid membranes (SILMs).6 It is a new class of materials also characterized as “solid liquid composites”.7 Inside the pore, a middle layer with close to bulk behavior exists together with a © XXXX American Chemical Society

surface layer where the structure and the dynamics are governed by the wall−fluid interaction.7 These layers depict different mobility from each other.6−10 In the case of highly confined SILMs, with a pore diameter close to double the size of the adsorbed layer, no middle-bulk layer is present. The influence of the wall and the suppression of the dynamics in the direction perpendicular to the wall become dominant.7 For the cations, that are not immobilized near the pore walls, an enhanced diffusivity compared to the bulk has been observed.8,9 Even a phase transition for a part of the film may take place. The existence of polymorphism in crystal structure or a transition to one-dimensional phase with different diffusivity have been recorded in the literature.6,7,11−13 The Tg was observed either to increase14 or to decrease.6 In cases of phase transitions captured or not by the differential scanning calorimetry (DSC) experiment,15 the remaining liquid phase is expected to depict a weak glass transition signature in the respective thermogram.6,14 The phenomenon is related to a reduced heat capacity step. This might explain the uncertainty in the evaluation of the Tg. More systematic studies on other confined organic liquids have shown suppression in Tm and Tg that is inversely proportional to the diameter of the nanopore.7 For the alkyl-methyl-imidazolium tricyanomethanide ([Rmim+][TCM−]) SILMs capped by silica, the experimentally evaluated diffusivity/permeability of CO2 is enhanced.6 MoreReceived: October 11, 2015 Revised: December 26, 2015

A

DOI: 10.1021/acs.jpcc.5b09947 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 1. (a) Single ions of the [Omim+][TCM-] IL and (b) the capped system.

same arguments of thermodynamic origin, for highly confined glassy materials (pore size comparable to the CRR) an immobilization on the solid surface may provide a stable, static liquid path. The suppression of the caging effect in the center of the pore should be related to an Arrhenius behavior in the activation energy.10,15,16,28−32 Atomistic simulations on confined ILs are relatively limited. Molecular dynamics (MD) simulations on such systems can predict the properties of SILMs and test the accuracy of theoretical models. In this case, a reliable force-field to describe accurately the bulk behavior of the specific family of ILs is needed.33−40 The accuracy of such a work is usually assessed in comparison to experimental data for specific properties (such as diffusivity, density, etc.). Furthermore, the simulation of the confinement conditions is a challenging task. Since the use of a grand canonical ensemble for these large, highly charged ions is difficult, other approaches have been employed. In most cases, simulations in the canonical ensemble (NVT) are preferred. According to the methodology followed,41−43 the mean density in the simulation box of the confined system is assumed independent from the wall size, having a constant value usually equal to the density of the bulk. In this way the excluded volume of the surface is not taken into account accurately.44,45 Hence, a potential introduction of additional stresses would affect the dynamics of the confined ILs. Also, it is not ensured that the confined film is in equilibrium with the bulk outside the pore. Another ensemble used is the isothermal - isobaric ensemble (NPT).46,47 In this case an atomistic description of the wall is more convenient. All diagonal components of the stress tensor (Pxx, Pyy, and Pzz) relax independently and a change in the box shape is allowed.46 A third approach is to test different loadings of the pore that may affect the dynamics of the ILs. By reducing the concentration of the ILs,48,49 a stiffening of the adsorbed layer has been observed, while upon increasing the ILs loading a decrease in the characteristic residence time of the ions on the

over, some dielectric experiments on highly confined ILs and other glassy materials indicated an Arrhenius thermally activated process.10,16 In this case, the activation energy in the most confined system depicted a minimum.10 The studied family of [Rmim+][TCM−]5,6 RTILs belongs to a wider category of materials that undergo glass transition. According to the Adam−Gibbs theory17−19 for the glass transition, at low temperatures the molecules of glassy materials move in cooperative way. Local fluctuations in the enthalpy and the creation of cooperative rearranging regions (CRRs) allow diffusion.17−19 On the other hand, the basic idea of the “free volume theory”18−20 claims that the molecules can achieve liquid mobility only when a minimum of space is accessible. The kinetic arguments consider the glassy materials as if they are below the Tm (supercooled liquids).18,19,21 Based on this model, the difficulties in diffusion (reflected in the nonArrhenius behavior) and the dramatic increase in viscosity are interpreted by the assumption that in equilibrium the material is in the solid phase. The existence of free energy barriers explains the experimentally observed state. Semicrystallinity in polymers is a characteristic example of such a case. Traditional kinetic arguments do not predict a facilitation in the diffusion and the dynamics in confined glassy materials.14 Confinement of complex systems introduces increased relaxation times. At an attractive surface, slow adsorption/ desorption kinetics explain an additional time scale separation.15,22,23 Nevertheless in some cases the high confinement of complex systems has been related to a non-Einstein decrease in viscosity.24,25 By modifying the Sanchez−Lacombe equation of state26,27 (SL EoS) one of the authors presented a model that explains the creation of CRRs with free volume arguments.15 The fluctuations in the pairwise energy embedded in a new EoS15 may explain the non-Arrhenius behavior of glassy bulk materials in the glass transition region. Since the dynamic CRRs are formed in order to facilitate the diffusion15 their size should be comparable to the size of the molecule. According to the B

DOI: 10.1021/acs.jpcc.5b09947 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

error procedure) equal to 2.6 Å. This value satisfies our assumptions of equilibrium with the bulk for the capped40 system at the temperature of 398.15 K. Then, we have kept fixed σw for each capped system, for all other temperatures. This approach results to a larger volume of the simulation box in the case of the capped25 system than the volume of the capped40 system. All simulations were performed using the LAMMPS software.52 The initial configurations for the capped, supported and bulk systems were obtained by the amorphous builder53,54 as implemented in MAPS software55,56 and underwent energy minimization using the conjugate gradient method. Then for the bulk systems a 50 ns run in the NPT ensemble using the Nosé-Hoover thermostat and barostat followed in order to obtain the average density at 398.15 K. The pressure damping parameter was 1000 fs while the temperature damping parameter was 100 fs. In addition we have taken advantage of the rRESPA multitime scale integrator57 as implemented in LAMMPS.52 The innermost loop (bond, angle, dihedral and improper forces) had a time step equal to 0.5 fs, the intermediate loop (pair forces) had a time step of 1 fs and the outermost loops (long-range forces) had a time step of 2 fs. The choice of these multiple time steps enabled the conduction of fully stable simulations at all cases studied. The cooling rate of the bulk systems was 50 K per 50 ns in the NPT ensemble. In the capped systems, the cooling rate in the NVT ensemble remained the same and the volume in the target temperature was predefined (eq 1). The equilibrium was identified by the energy equilibration criterion and the symmetry in the density profile for the capped systems. Production runs of the order of 100 ns were performed in the NVT ensemble using the Nosé-Hoover thermostat. Additional runs in the microcanonical ensemble (NVE) were performed at 298.15 K, as a means of verification that the dynamics of the system was properly reproduced in the NVT ensemble. The same procedure was followed for the supported systems as well. 2.2. Force Field. We have used an all-atom force field for the representation of the ILs.34,35,37,58 The expression for the intramolecular interactions (bond stretching, bond bending, dihedral angle distortion and improper potentials) is of the form

surface and an increase in their self-diffusivity (close to bulk properties) has been detected.48 Most molecular simulations of confined ILs42,43,47,49 revealed an interfacial layer consisting mainly of cations and extending up to the first density minimum. The structure and the dynamics of the adsorbed layer have been also confirmed in the case of a one solid surface film.50 In the present work, we study the confinement effects on the structure and the dynamics of 1-methyl-3-octylimidazolium tricyanomethanide [Omim+][TCM−] IL using MD simulations. In section 2, we describe a new methodology that allows us to follow the temperature dependence of the density in a confined system. We give the simulation parameters that reproduce the studied IL/silica interface. Detailed information about the force field parameters of the bulk IL are given elsewhere.51 In section 3, we present our results and propose a model for the interpretation of the observed behavior. The examination of our simulation outcomes is performed in comparison with a previous experimental study6 on the same system. Finally, we end with the conclusions of the study.

2. SIMULATION DETAILS 2.1. Methodology. The systems examined consisted of 100 pairs of [Omim+][TCM−] ions, confined by two solid silica surfaces (labeled as capped systems) as shown in Figure 1. The capped systems were compared to the respective cases in which one of the two silica surfaces was substituted by vacuum (labeled as supported systems). Interwall distances equal to 25, 40, 60, and 80 Å (labeled as capped25, capped40, capped60, and capped80, respectively) were examined, focusing mainly on the detailed analysis of the two systems of highest confinement (capped25). The density profile (density oscillations) near the surface obtained from the supported systems was used as guide for the adjustment of the system’s dimensions in the xy plane (parallel to the surface). Additionally, in the case of the capped40, capped60, and capped80 systems the requirement of equilibration with the bulk (outside the pore) was also satisfied by the assumption of a middle layer with the same density with the one of the bulk. The bulk pressure at all temperatures was set equal to 0.1 MPa. Since the shortest interwall distance (25 Å) used is about twice the distance in which the oscillations in the cation−anion radial distribution function vanish, the simulated system can be considered close to equilibrium with the bulk. The films were studied at 253.15, 298.15, 353.15, and 398.15 K. In all cases, the confinement effects on the calculated properties were examined in direct comparison with bulk simulations at the same temperature and pressure conditions. Based on the assumptions already mentioned, a relation between the volume of the bulk and the volume of the confined system should take into account the excluded volume due to the presence of the opposite silica walls.44,45 Assuming that a linear relationship is valid for a wide range of interwall distances and temperatures (much higher than Tg), a good approximation is given by the equation xcappedycapped =

Vbulk (zcapped − 2 σw )

Uintra =

∑

kb(b − b0)2 +

bonds

∑

kθ(θ − θ0)2

angles n=1

+

∑ ∑ kχ[1 + cos(nχ − δ)] dihedral

+

∑

4

k ψ (ψ − ψ0)2 (2)

impropers

where b, θ, χ, and ψ denote bond, angle, dihedral, and improper angle respectively, the subscript “0” refers to the equilibrium values, n is the multiplicity parameter of the dihedral angle, and δ is the phase shift of the dihedral potential over the full range of rotation. The intermolecular interactions (electrostatic and LennardJones (LJ)) are given as

(1)

where xcapped, ycapped are the x, y dimensions of the simulation box of the capped system (wall dimensions), zcapped is the interwall distance, Vbulk is the volume of the bulk system as calculated by NPT simulations and σw is the average excluded volume parameter. The parameter σw was set (by a trial and

⎧ ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ qiqj ⎫ σij ⎥ ⎪ ⎪ ⎢ σij ⎬ = ∑ ∑ ⎨4εij⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥ + r ⎝ rij ⎠ ⎦ 4ε0rij ⎪ i=1 j>1 ⎪ ⎭ ⎩ ⎣⎝ ij ⎠ n−1 n

Uinter

C

(3)

DOI: 10.1021/acs.jpcc.5b09947 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C where qi stands for the partial charges, ε0 is the vacuum permittivity, and ε and σ are the LJ parameters. The pairwise interactions between atoms separated by less than three covalent bonds were excluded. Also, the 1−4 Coulombic interactions for atoms separated by three bonds were scaled by a factor of 0.5, while the 1−4 LJ interactions were fully taken into account. A cutoff distance of 12 Å was used for both electrostatic and van der Waals interactions. Long range interactions in both cases were treated according to particle− particle particle-mesh solver as embedded in LAMMPS52 for both bulk and slab/film cases. The set of force field parameters that was used is presented in ref 51. Charge distributions for both ions were calculated quantum mechanically (QM) using various schemes and the anion parameters were determined.51 A total ionic charge of ±0.75e was introduced in this force field as dictated by the QM calculations. Bonded and LJ parameters for the cations were obtained from the literature58,59 with some slight modifications. The accuracy of the optimized force field has been validated for the whole [Cnmim+][TCM−] (n = 2, 4, 6, 8) IL family, in both NVT and NVE ensembles. We have employed both LAMMPS52 and NAMD60 software. The silica surface was described by an integrated LJ potential of the form42 Vwall(z) =

2πσw 3εw ρ ⎡ 2 ⎛ σw ⎞9 ⎛ σw ⎞3⎤ ⎢ ⎜ ⎟ −⎜ ⎟⎥ ⎝z⎠⎦ 3 ⎣ 15 ⎝ z ⎠

capped40 system). Similar behavior is observed at all temperatures examined (253.15, 298.15, 353.15, and 398.15 K) and are not shown here. The trends were similar for the capped60 and the capped80 systems. In Figure 2, the densities of each counterion are also shown. Based on the minima in the density profile we may define layers that enable the study of the structure and the dynamics as a function of the distance from the surface. As shown by the vertical dot lines, in the capped40 system we distinguish the layerI which corresponds to distances up to 16 Å from the wall and the layerII which corresponds to distances larger than that. Additionally, for both capped40 and capped25 systems we introduce a second discretization, according to which distances up to 5.8 Å from the wall belong to the adsorbed layer (layerAds) and greater distances belong to the free layer (layerFre). The density profiles in both systems are obtained from a 100 ns MD run and are not symmetrized. A symmetrization is applied only in specific cases in which the centers of mass are involved, in order to refine the profiles and gain better statistics. Considering the symmetry of the density profile as indication of equilibration, it is interesting to note that equilibration is reached much faster in the capped25 system compared to the other capped systems with larger interwall distances. In Figure 3, the radial distribution function (g(r)) of the center of masses of the cations is presented. In the bilayer

(4)

where σw = 3.0 Å, εw = 0.8 kJ/mol, and 4πρ = 0.5 sites/Å3. In all interactions we have used the Lorentz−Berthelot mixing rule.

3. RESULTS AND DISCUSSION 3.1. Structure. The mass density profiles of the capped25 and capped40 systems at 298.15 K are depicted in Figure 2. The values in bulk are shown with dashed horizontal lines. In the case of the capped40 system, the middle layer has a mean density close to the bulk. In addition, the density oscillations in the capped40 and capped25 systems overlay the density profile of the supported system (the same x,y box dimensions with the

Figure 3. Radial distribution function of the cation−cation centers of mass. We show results for the two layers in the capped40 system, the capped25 system and the bulk at 298.15 K. In all cases, the normalization was done based on the mean density in bulk.

(layerI and layerII) approach of the capped40 system, we assume that the central ions are in the middle (±3 Å in the zaxis) of each layer. Their neighbors are anywhere in the film. Although the immediate ion neighborhood in not uniform (especially for the layerI), the g(r) presents a mean picture of the neighboring ions, which can be used for the interpretation of the observed stresses and dynamics in the adsorbed and the middle layers of the silica pore. In the curve of the layerI of the capped40 system, we observe a first peak at 5.7 Å which does not exist in the bulk. It reflects the increased cation density near the surface.42,43,47,49 Most of the neighboring cations are located in the x, y directions. Moreover, at greater distances the same curve depicts values lower than unity. In layerII of the capped40 system, the immediate neighborhood resembles the

Figure 2. Density profiles of the capped40 (solid red line), capped25 (solid blue line) and the supported (solid black line) systems, at 298.15 K. z = 0 corresponds to the wall. Also, we depict separately the densities of the each ion (squares for the cations and circles for the anions). In all cases we compare with the respective values in bulk (horizontal dashed black lines). The vertical dot lines separate the systems in layers. D

DOI: 10.1021/acs.jpcc.5b09947 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

reaches values close to the ones in the bulk. In this way, due to bulk behavior the polar and nonpolar domains are not distinguished in a pronounced way. In case of the capped25 system, as the dimensions of the gap are comparable with the dimensions of the formed anion−cation ordering, we observe a depletion of tails in the middle layer. The concentration of the cation rings and the anions is higher compared to the bulk in the center of the film. The first peak in the anion density in the capped25 system is slightly shifted toward the core of the pore compared to the capped40 system. Moreover, the structure in the center of the pore of the capped25 system resembles the structure of the layerI of the capped40 system. The overall picture is in agreement with the results from the g(r). It explains the first peak in the cation−cation g(r), in the capped40_layerI curve. Also, the depletion of the cation centers of mass at ∼13 Å from the surface in the capped40 system justifies the reduced values of the g(r) at short distances in the curve capped40_layerII, in Figure 3. We mention that the concentration profile is imposed by the wall, which we have assumed carries no charges. Nevertheless, even a wall with sites of electrostatic interaction with the IL is expected to undergo a partial shielding by a monolayer of adsorbed ions, possibly cations rings, forming hydrogen bonds with the silica.49,61 Especially in case of a phase transition of the adsorbed layer,6 the depicted picture (Figure 4) could be representative of the structure next to the immobilized/adsorbed layer. Moreover, the aggregation of the tails on the wall was also confirmed for the case of a wall with less adsorption strength49 that exerts the same force in all atoms (not shown). The increased concentration of cations on the surface is in agreement with other computational studies in various confined ILs.42,47,50 It supports the model proposed6 for the experimentally observed facilitation of CO2 diffusion in the [Rmim+][TCM−]/silica SILMs under extreme confinement. According to the same model the anions are in the center of the film and the CO2 diffuses through a “hoping” mechanism from anion to anion.6 In Figure 5a, the stress profiles of the capped40 and capped80 systems are shown. The total stresses (all contributions) describe the formation of regions, defined by peaks, where the pressure is positive. In order to explain this behavior, we have plotted the total stresses of a supported film with the same xy box area as in the respective capped40 system. We observe that the stress oscillations in the supported system depict the same behavior. Consequently, we may assume that it is the few nm width of the film that is responsible for the layering. It reflects the spatial organization of the ILs (size ∼12 Å) in the z-axis. In our simulation of the confined ILs, no additional stresses (negative or positive) are induced. Moreover the IL/IL contributions of the stresses (after subtracting the wall contributions) are in agreement with the density profile. The IL/IL contributions, at ∼5 Å from each surface depict (not shown) positive pressure.62 In the capped25 system (Figure 5b), the gap is equal to twice the size of the clustering. Diffusion in this case does not promote the formation of regions in the middle layer. The low negative pressure values are in agreement with the behavior in the respective supported system (Figure 5b). It is worth noting that the standard deviation in the capped25 system is lower compared to the case of capped40 system, although the length of both runs is exactly the same. The orientation of the confined [Omim+][TCM−] at 298.15 K is examined in Figure 6. A second order Legendre polynomial has been used for the calculation of the order

one of the bulk. However, at distances greater than 10 Å we observe values higher than unity. The g(r) is affected by the intense density oscillations in the capped40 system. In the case of the capped25 system, we assume there is no middle-bulk layer (Figure 2). The capped25 curve represents a scan over cation−cation neighbors in the whole film. It depicts an intermediate picture between the capped40_layerI and capped40_layerII curves. Concerning the cation−anion g(r), Figure S1 in the Supporting Information (SI) shows that in all cases, in the capped systems the first peak is higher compared to the bulk. Even in the shortest interwall gap (extreme confinement), the first peak in the cation−anion g(r) appears at the same distance.46 Moreover, in the capped25 system the same peak depicts the highest value. The anion−anion g(r) is plotted in Figure S2 in the SI. It is shown that the capped systems preserve the two peaks. The influence of the surface, as illustrated in curve capped40_layerI, becomes evident in the middle region of the capped25 system. The first peak in the capped25 curve has the highest value. It describes an increased density of neighbor anions caused by the organization due to confinement. On the other hand in layerII of the capped40 system the same peak depicts a value close to the bulk. The fluctuations in the anion−anion g(r) extend to distances greater than the 12 Å. The spatial distribution function having as origin in z-axis the lower silica wall (Figure 1b) is presented in Figure 4. In all cases

Figure 4. Spatial distribution of the anion, the cation and the cation “ring” centers of mass at 298.15 K. Also, we present the distribution of the terminal carbon atoms in the cation tail (CT3, “free-end”). The results for capped40 system are depicted with red, while the results for the capped25 system are drawn with blue. The walls are placed at ±20 and ±12.5 Å on the z-axis for the capped40 and capped25 systems, respectively.

the width of the layer for the discretization in z-axis is 1 Å. We observe that the centers of mass for both ions have a peak at almost 8 Å from the surface, for the capped40 and capped25 systems. The nonpolar tails tend to aggregate on the surface. At further distances in layerI of the capped40 system, we can detect the formation of polar domains.2−4 The rings of the cations appear at the distances of almost 12 Å from the wall. At the same distance we also observe a second peak in the concentration of the anions. In the layerII of the capped40 system the concentration of the centers of mass of both ions E

DOI: 10.1021/acs.jpcc.5b09947 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

parameter considering an angle θ between the z-axis and a specific vector. A value equal to 1 describes an alignment parallel to the z-axis, while a value equal to −0.5 describes an alignment parallel to the surface. A value of zero corresponds to a no order case. In our analysis, three vectors were examined: (i) the vector perpendicular to the cation’s ring plane, (ii) the vector along the alkyl tail, defined by the first and the terminal carbon atom, and (iii) the vector perpendicular to the anion plane. In both systems (capped40 and capped25) on the surface (layerAds), we observe that the small fraction of ions located in this layer exhibits increased order and conformations parallel (cation ring plane and anion plane) to the silica wall.43 The order parameter in layerAds could give us a sense about the experimentally6 observed phase transition on the silica surface. In the layerI of the capped40 system, at distance ∼8 Å from the silica wall where the concentration of both ions is increased, the cation ring plane is slightly tilted perpendicular to the surface.46 Moreover at the same distance the tail vector shows a maximum value in the order parameter (∼0.6).46 Our predictions for the orientation of the cation with the use of a wall potential are in agreement with other studies in which the silica surface is described atomistically.46 It points out the entropic character of the aggregation of the cation on the surface. Also, the results are in general agreement with the results from previous experimental studies.61,63,64 The anion plane, at ∼8 Å from the wall, has still a tendency for an orientation parallel to the surface but in a less pronounced way. In the middle layer (layerII) of the capped40 system there is no order in anions. For the cations, the influence of the surface in layerII also tends to zero. However, for the capped25 system, in which no middle layer with bulk behavior exists, the orientation observed in the layerI (not in layerAds) of the capped40 system is extended in the whole capped25 system. The ring plane in the middle layer is slightly tilted perpendicular to the wall, while the tail vector is also slightly oriented parallel to the z-axis. The anions obtain less ordered conformations. Using the eigenvalues of the gyration tensor,65 a picture about the shape of a molecule can be generated. In Figure 7, we show the dimensions of the more deformable ion which is the cation. The radius of gyration for both capped40 and capped25 systems indicates extended conformations on the surface. This

Figure 5. (a) Total stress profiles of the capped40 and the respective supported system and also the capped80 system, at 298.15 K. (b) Total stress profiles of the capped25 and the respective supported system, at 298.15 K. The supported systems have the same xy box area with the respective capped systems.

Figure 6. Orientational-order parameter (P2(cos θ)) of the ring plane (perpendicular vector), the tail vector (first and terminal tail carbons) and the anion plane (perpendicular vector) at 298.15 K. The results for the capped40 system are drawn with red, while the results for the capped25 system are drawn with blue.

Figure 7. Radius of gyration (eq s1 in the SI) of the cation for the capped40 (red squares) and the capped25 (blue squares) systems, at 298.15 K. With dash line we depict the bulk mean value. F

DOI: 10.1021/acs.jpcc.5b09947 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

capped40 and capped25 systems. In the layerI, the decorrelation time is higher compared to the layerII (Table 1). After 10 ns, the curve capped40_layerI depicts a plateau,

is consistent with the order parameter depicted in Figure 6. Although in the capped40 system a middle layer with bulk cation dimensions exists, for the capped25 system a slight elongation of the cation is indicated. The anisotropy, asphericity and acylindricity defined according to eqs s2−s4 in the SI and presented in Figure S3 reveal that on the surface, the tendency for anisotropy and asphericity becomes more pronounced for both capped25 and capped40 systems. At distances just beyond the layerAds (∼6 Å from the wall), in both systems the cations take configurations less anisotropic, more spherical and less cylindrical than in the bulk case. As expected, the planar structure of the anion remains unchanged in the whole film (⟨Rg2⟩ = 3.65 Å2). 3.2. Dynamics. The description of the silica wall was implemented with the use of an integrated LJ potential (eq 4). This representation neglects the wall roughness. Nevertheless, our simulation introduces42,66 in an adequate way the effect of confinement in one dimension. The diffusion in the capped systems is not isotropic as in the bulk case. We have chosen to examine the mobility in the z-axis separately from the mobility in x,y axis (2D displacement). The restriction in the dynamics due to the wall potential is studied by means of the displacement of the centers of mass in the z-axis, according to the incoherent scattering function22 Fs(qz , t ) =

1 M

Table 1. Relaxation Times for the Autocorrelation Functions Presented in Figures 8−11 and S4 of the SIa τc [ps]

a

k=1

Fs(qz,t) for cation

Fs(qxy,t) for cation

P2(t) for CN7B-NR1

P2(t) for CT2-CT3

P2(t) for CCM-CN

layer I layer II capped25 bulk

16 180 11 953 13 095 3000

2974 1978 1683 3317

2479 675 831 1327

22 29 19 38

231 221 135 322

The β parameter is given in Table S1 of the SI.

while the curve capped40_layerII continues to drop. However, at short times (