Structure and Dynamics of Octamethylcyclotetrasiloxane Confined

Feb 25, 2016 - Influence of surface commensurability on the structure and relaxation dynamics of a confined monatomic fluid. Vadhana Varadarajan ...
0 downloads 0 Views 5MB Size
Article pubs.acs.org/JPCB

Structure and Dynamics of Octamethylcyclotetrasiloxane Confined between Mica Surfaces V. Vadhana† and K. G. Ayappa*,†,‡ †

Department of Chemical Engineering, Indian Institute of Science, Bangalore - 560012, India Centre for Biosystems Science and Engineering, Indian Institute of Science, Bangalore - 560012, India



ABSTRACT: Using a molecular model for octamethylcyclotetrasiloxane (OMCTS), molecular dynamics simulations are carried out to probe the phase state of OMCTS confined between two mica surfaces in equilibrium with a reservoir. Molecular dynamics simulations are carried out for elevations ranging from 5 to 35 K above the melting point for the OMCTS model used in this study. The Helmholtz free energy is computed for a specific confinement using the two-phase thermodynamic (2PT) method. Analysis of the in-plane pair correlation functions did not reveal signatures of freezing even under an extreme confinement of two layers. OMCTS is found to orient with a wide distribution of orientations with respect to the mica surface, with a distinct preference for the surface parallel configuration in the contact layers. The self-intermediate scattering function is found to decay with increasing relaxation times as the surface separation is decreased, and the two-step relaxation in the scattering function, a signature of glassy dynamics, distinctly evolves as the temperature is lowered. However, even at 5 K above the melting point, we did not observe a freezing transition and the self-intermediate scattering functions relax within 200 ps for the seven-layered confined system. The selfdiffusivity and relaxation times obtained from the Kohlrausch−Williams−Watts stretched exponential fits to the late α-relaxation exhibit power law scalings with the packing fraction as predicted by mode coupling theory. A distinct discontinuity in the Helmholtz free energy, potential energy, and a sharp change in the local bond order parameter, Q4, was observed at 230 K for a five-layered system upon cooling, indicative of a first-order transition. A freezing point depression of about 30 K was observed for this five-layered confined system, and at the lower temperatures, contact layers were found to be disordered with long-range order present only in the inner layers. These dynamical signatures indicate that confined OMCTS undergoes a slowdown akin to a fluid approaching a glass transition upon increasing confinement, and freezing under confinement would require substantial subcooling below the bulk melting point of OMCTS.



INTRODUCTION Over the last two decades, the physics of fluids confined to nanoscale dimensions have been investigated using a variety of experimental methods, molecular simulations, and theory.1−11 It is now well accepted that a fluid confined between two infinite surfaces can remain a liquid or freeze and undergo solid−solid transitions.9−18 The state of the confined fluid is naturally dependent on the fluid density and interactions of the fluid with the confining surface. Recent simulations show that a confined fluid can also exist as a glass11,19−21 and varying degrees of confinement can induce re-entrant glass transitions.22 Surface force experiments have been used to study the static and dynamic response of fluids confined to molecular separations between two mica surfaces. Confinement effects on octamethylcyclotetrasiloxane (OMCTS), a model nonpolar fluid, have been extensively studied in surface force experiments.5,6,23−25 Since mica surfaces are atomically smooth, solvation force oscillations as a consequence of fluid layering have been observed both in experiments6,23,25−27 as well as in simulations.28−30 These experiments also reveal that the viscosity increases by several orders of magnitude for a © XXXX American Chemical Society

confined system when compared with their corresponding bulk values. Solvation forces have also been observed in atomic force experiments, indicating that layering can also occur when a fluid is confined between a tip and an infinite surface.28,31 Surface force experiments with OMCTS are carried out at a temperature of 298 K which is about 8 K above the melting point of bulk OMCTS.32 Under these conditions, the confined fluid is in equilibrium with a liquid reservoir close to the liquid−solid transition of OMCTS. The question then arises as to the state of OMCTS under confinement. The interpretation of static and/or oscillatory surface force measurements has resulted in two views for the state of confined OMCTS. Does the fluid undergo a liquid to solid transition upon confinement,6,23 or is the transition more gradual, akin to a dynamic slowdown observed in a glass transition?25,33 Interpretation of the measurement was further complicated by the presence of Pt nanoparticles deposited on the mica surface during sample Received: December 31, 2015 Revised: February 14, 2016

A

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

Article

The Journal of Physical Chemistry B

obtained from the bath simulations, NVT molecular dynamics simulations are carried out in periodic systems to study the structure and dynamics as a function of surface separation. In addition to density, pair, and orientational distribution functions, we compute the mean squared displacement and the self-intermediate scattering functions as a function of surface separation. The latter quantity provides a signature of the effects of caging and relaxation associated with the dynamic slowing down of a liquid approaching the glass transition. Using a mode coupling analysis, we show that both the diffusivity and late α relaxation of the scattering function possess a power law divergence as a function of the packing fraction of the confined fluid. To further investigate the possibility of freezing under confinement, we compute the Helmholtz free energy using the two-phase thermodynamic (2PT) method. In contrast to other methods, the entropy is computed from the density of state (DoS) in the 2PT method by a decomposition into translational, rotational, and vibrational components. In the 2PT method, as the name suggests, the DoS for each of these components is obtained by a decomposition into gas and solid contributions. The method had been validated for both atomic and molecular fluids as well as for confined molecular fluids such as water.39−47 In contrast to umbrella sampling and other methods used to compute the free energy,48 the 2PT method makes use of the molecular dynamics trajectory of the fluid at a given state point. This makes it particularly attractive for sampling confined molecular fluids without the complexity of locating a reversible intergration path. Alternate methods for monatomic confined fluids have been based on a Landau free energy approach.13 Using the 2PT method, we obtain the free energy for a five-layer confined OMCTS system which is gradually cooled and show that freezing occurred at temperatures about 30−35 K below the bulk melting point of OMCTS.

preparation, and several experiments have been carried out to resolve this aspect.34−36 The normal force is measured as a function of mica surface separation during static experiments, or changes in phase and amplitude are recorded during an oscillatory rheological measurement. The state of the confined fluid is then inferred from a suitable model for the rheological response of the fluid. Since the confined fluid structure is not directly probed in a conventional surface force apparatus, a direct structure property correlation is not possible. Variations such as the X-SFA and e-SFA attempt to directly correlate the fluid structure with the measured force response.37 Molecular simulations provide an atomic level perspective of the fluid arrangement which cannot always be probed using experiments, especially under conditions of strong confinement. Monte Carlo and molecular dynamics (MD) simulations bridge this gap by a direct sampling of the atomic coordinates given an appropriate intermolecular force field. Simulations for OMCTS, until recently, have been based on a spherical monatomic representation which predict the effects of layering and oscillatory solvation forces observed in experiments.10 However, the spherical model assumption ignores shape and packing effects associated with a fully atomistic model for OMCTS which accounts for the somewhat flattened structure of the molecule. Shape anisotropy effects have been illustrated by computing the intermolecular potential as a function of different orientations by Matsubara et al.38 Molecular models for OMCTS have recently been used to study the structure and dynamics of OMCTS between mica surfaces,30,38 as well as between an AFM gold tip and a mica surface28 using molecular dynamics simulations. These studies reveal a distinct decrease in diffusivity and increase in rotational relaxation times upon increasing confinement;28 however, these studies do not show evidence of freezing except at extreme confinement, when two OMCTS layers are present.28 Since these molecular models of OMCTS have been parametrized with the liquid density and diffusivity of OMCTS, the accuracy of the melting point predicted by these models has not been investigated. The melting point is an important consideration, as this has a direct bearing on conclusions drawn about the state of a fluid under confinement. This assumes special significance for OMCTS, since surface force experiments are carried out at an elevation of 8−10 K from the melting point of OMCTS which is 290 K.32 In this manuscript, we examine structural and dynamical quantities of OMCTS confined between mica surfaces, with the primary focus on determining whether OMCTS behaves like a fluid approaching a glass transition or undergoes a freezing transition upon confinement. Molecular dynamics simulations are carried out using the rigid molecular model developed by Matsubara et al.38 Using the liquid−solid interface method, we estimate the melting point for this molecular model to lie between 260 and 265 K. In order to accurately reflect the influence of confinement with respect to the bulk melting point, we carry out simulations at temperatures ranging from 265 to 300 K, thereby encompassing the range of temperature elevations in surface force experiments. Unlike previous studies with the molecular model for OMCTS, we perform molecular dynamics simulations with two mica surfaces immersed in a liquid bath of OMCTS to mimic the surface force experiments. This method is expected to yield accurate equilibrium densities for the confined fluid as a function of surface separation. Further, this is expected to overcome the drawback of limited sampling at high densities typically encountered in GCMC simulations for molecular systems. Once the number density is



THEORY AND SIMULATION PROCEDURE Molecular dynamics simulations are carried out with OMCTS modeled as an eight-site neutral molecule with 12−6 LennardJones interactions between the eight methyl groups treated as united atoms. We used the model OMCTS parameters reported as model B in the parametrization by Matsubara et al.,38 since the bulk density and diffusivity obtained with this parametrization at 300 K and 1 atm match with experimental data for OMCTS. The Lennard-Jones parameters for the methyl groups are σ = 3.54 Å and ε = 0.39 kcal/mol. Molecular dynamics simulations were carried out using GROMACS 4.5.6, and initial configurations were obtained using PACKMOL.49 In order to carry out MD simulations within the rigid body framework, we used restraint algorithms where the bonded, nonbonded distances and dihedral angles are maintained within a specified tolerance from the initial molecular structure. The restraint algorithm is based on applying potentials to the internal degrees of freedom to prevent deviations from the molecular structure. The initial molecular structure was adapted from the structure reported by Matsubara et al.38 It is important to note that these additional potentials are not part of the force field. Two kinds of restraints were used for simulating a rigid body. • Distance Restraints: A penalty to the potential is added when the distance between any pairs of atoms including both bonded as well as nonbonded distances exceeds a B

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

Article

The Journal of Physical Chemistry B threshold value. The force constant used was 50 000 kJ mol−1 nm−2. • Dihedral Restraints: These are used to restrain the dihedral angle formed by four atoms or groups of the molecule. The force constant used in the simulation is 500 kJ mol−1 rad−2. Bulk Fluid Simulations. 216 OMCTS molecules were simulated at 300 K with a time step of 2 fs in a cubic box with an edge length of 4.83 nm using a Berendsen thermostat for 2 ns. The time constant for the Berendsen thermostat was 0.2 ps. NPT simulations were carried out for 50 ns from the output trajectories of the NVT simulations using a Berendsen thermostat and a Berendsen barostat with time constants of 0.2 and 0.3 ps, respectively. The final density obtained from NPT simulations was 957.6 kg m−3 which was within 1% of the experimentally reported density of 948 kg m−3 for OMCTS.38 The radial distribution function (Figure 1a) based on the intermolecular methyl−methyl distance is in good agreement with the data reported by Matsubara et al.38 Figure 1b illustrates the mean squared displacement based on center of mass positions of the OMCTS molecules. The mean squared displacement shows a transition from a ballistic regime at short times to a diffusive regime at long times. The self-diffusion coefficient is obtained using the Einstein relation D = lim

t →∞

⟨|r(t ) − r(0)|2 ⟩ 6t

(1)

where r = exrx + eyry + ezrz and the angular brackets denote averages over both shifted time orgins as well as the number of particles. The self-diffusivity is 3.22 × 10−6 cm2 s−1 which compares well with the experimental value of 4.03 × 10−6 cm2 s−1.50 The self-intermediate scattering function, Fs(k, t), computed with the center of mass positions using N

Fs(k, t ) =

1 ⟨∑ exp(i k· (rj(t ) − rj(0)))⟩ N j=1

(2)

−1

with the wavenumber k = 2.618 nm is illustrated in Figure 1c which shows a single relaxation time typical of simple liquids. In the long wavelength and low frequency limit, the selfintermediate scattering function can be expressed as

Fs(k , t ) = exp( −Dk 2t )

(3)

Figure 1. Structural and dynamical quantities from bulk liquid simulations of octamethylcyclotetrasiloxane at 300 K. (a) Pair correlation functions based on intermolecular methyl groups in OMCTS at 300 K. The solid line represents data from this study, and the triangles are data extracted from Matsubara et al.38 (b) Mean squared displacement for OMCTS at 300 K. (c) Self-intermediate scattering function for OMCTS at 300 K. The dotted line corresponds to the exponential fit, Fs(k, t) = 0.822 exp(−0.00667t).

where D is the self-diffusion coefficient. For the smallest wavenumber k = 2.607 nm−1 based on half the box size, we obtain D = 3.24 × 10−5 cm2 s−1, which compares well with the self-diffusion coefficient obtained from the mean squared displacement. We point out that the value of D = 2.666 × 10−5 cm2 s−1 reported by Matsubara et al.38 is slightly lower than the value obtained in our studies. This could be attributed to the restraint algorithm in our study to keep the rigid body assumption for the model. Nevertheless, the diffusivity obtained by us is in closer agreement with the experimental data. As a further confirmation, we also ran a Monte Carlo simulation with an inhouse code at 300 K and 1 bar for the rigid OMCTS model. The density was in excellent agreement with that reported by Matsubara et al.38 Melting Point of OMCTS. Since the surface experiments are carried out with the liquid OMCTS at about 8° above the melting point which is about 290 K,32 it is important to assess the melting point predicted by the molecular model for

OMCTS. Since the melting point was not one of the criterion used to parametrize OMCTS,38 we determined the melting point of OMCTS by a direct simulation of the solid−liquid interface51−53 which is expected to yield the melting point to within an accurcacy of 3−5 K. NPT simulations were carried out for a system consisting of 960 molecules of crystal OMCTS in contact with 960 molecules of liquid OMCTS at 1 atm with a time step of 1 fs. The crystal structure was built from the lattice parameters reported by Steifink et al.54 The initial liquid structure was created by placing molecules at random orientations in a volume similar to that used for the solid C

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

Article

The Journal of Physical Chemistry B phase. The temperature was fixed using a Nosé−Hoover thermostat, and pressure was controlled using Parrinello− Rahman barostat. The time constants were obtained by validating the method for the liquid−solid coexistence of a Lennard-Jones fluid. Using a time constant of 5 ps for the thermostat and 2 ps for the barostat, the melting point obtained from the simulation was 81.5 K, which is in good agreement with the melting point of 84 K55 for argon, obtained from the Lennard-Jones phase diagram. Simulations were run for 100 ns with the above-mentioned time constants. Since the Nosé− Hoover thermostat can produce oscillatory temperature fluctuations, the system was initially equilibrated using a Berendsen thermostat for 2.5 ns before using the Nosé− Hoover thermostat. Melting or freezing can be monitored from the total energy, which is the sum of the potential and kinetic energies (neglecting the internal degrees of freedom). If the system temperature is maintained above the melting point, the crystal will melt and the thermostat supplies heat to the system. If the system temperature is maintained below the melting point, the liquid region will freeze and the thermostat will extract heat from the system. Figure 3a shows the evolution of

Figure 3. Determination of melting point using the solid−liquid coexistence method. (a) Potential energy of OMCTS crystal at different temperatures. The initial crystal structure corresponds to the lattice parameters at 223 K.54 (b) Variation of crystal and liquid thicknesses with time at different temperatures. Solid lines represent the variation of crystal thickness with time, and dotted lines represent variation of liquid thickness with time.

consisting of two surfaces held at a fixed separation was placed inside a bath of 48670 OMCTS molecules. The bath is a cube with an edge of 29.5 nm with periodic boundary conditions in x, y, and z directions. The system temperature is set to 300 K using a Berendsen thermostat with a time constant equal to 0.2 ps, and the time step used in the simulation is 2 fs. The simulations were run for 4 ns and data sampled every 0.4 ps. The size of the bath used in the simulation was decided on the basis of ensuring that the average OMCTS density away from the mica surfaces is 957.6 kg m−3 at 300 K. In all cases, we found that the number of particles between the mica surfaces equilibrated within about 2 ns and the equilibrium number density of confined OMCTS was obtained by averaging the particle numbers over the last 1 ns from a 4 ns simulation after removing edge effects. The influence of the edge was estimated by evaluating the pore density as a function of the distance from the mica edges. The densities reported here were obtained by counting the molecules within 2 nm (4 lattice parameters) away from each edge along the x direction and 2.7 nm (3 lattice parameters) away from each edge along the y direction. The resulting number densities for different surface separations investigated in this study are given in Table 1. Figure 2 illustrates a small density enhancement with respect to the bulk density at all surface separations. The density is also enhanced at lower temperatures for H = 5.3 nm. The equilibrium number densities obtained from the bath simulations were used in molecular dynamics simulations to

Figure 2. Time averaged number of molecules between the mica surfaces as a function of surface separations after removing edge effects for surface area ≈45.06 nm2 at 300 K. Data for simulations at lower temperatures for H = 5.3 nm is also illustrated. The solid line represents the number of molecules based on the bulk density of OMCTS at 300 K.

the total energy for the crystal OMCTS during the initial equilibration runs. The total energy at 260 K initially decreases and then reaches a plateau, indicating that heat is removed from the system. Energies at 270, 275, and 280 K increase with time, suggesting melting of the crystal. From the energy evolution, the melting point was estimated to lie between 260 and 270 K. In order to sharpen this estimate, we extract the thickness of the solid and liquid regions as a function of simulation time (Figure 3b) from the density distributions normal to the liquid−solid interface. From the variation of thickness in the solid and liquid regions, we estimate the melting point to lie between 260 and 265 K. Confined Fluid Simulations. The structured mica wall was constructed from a unit cell of muscovite mica as described by Malani and Ayappa,56 and similar potential parameters were used in this study. The mica walls consisted of 20 × 10 unit cells which make up a periodic cell of lateral dimensions Lx = 104.088 Å and Ly = 90.18 Å. The surface separation, H, between two mica walls is based on the distance between the oxygen layers of the top and bottom walls. An empty mica pore D

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

Article

The Journal of Physical Chemistry B

2D simulation consisted of 16 unit cells along x and 6 unit cells along the y direction. Initially, canonical ensemble, NVT simulations are carried out using Berendsen scaling with a time constant of 0.2 ps for 2 ns. After this initial equilibration, simulations were continued using a Nosé−Hoover thermostat with a time constant of 0.5 ps and a temperature of T = 300 K for a duration of 50−100 ns. We also carried out a few additional longer simulations for 140−150 ns for H = 2.11 and 2.83 nm and did not observe any differences in either the structural or dynamical quantities reported in this study. The density normal to the mica surface is evaluated using

Table 1. Pore Density Obtained for a Confined OMCTS− Mica System for Various Mica Surface Separations, H, and Temperatures H (nm)

T (K)

ρ (nm−3)

5.30 5.30 5.30 4.12 3.30 2.83 2.48 2.11 1.77

300 275 265 300 300 300 300 300 300

1.987 2.073 2.150 1.996 1.996 2.070 2.110 2.185 2.223

(

N z− ρ (z ) =

study the dynamics of OMCTS confined between surfaces periodically replicated in the x−y plane. The mica walls for the

Δz , 2

z+

AΔz

Δz 2

) (4)

where the angular brackets in the above equation represent the time averages, Δz is the bin thickness, and A is the cross section

Figure 4. Density distributions of confined OMCTS in the direction normal to the mica surface for various surface separations: (a) H = 5.3 nm, (b) H = 4.12 nm, (c) H = 3.3 nm, (d) H = 2.83 nm, (e) H = 2.11 nm, (f) H = 1.77 nm. E

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

Article

The Journal of Physical Chemistry B area of the mica surface. In all cases, 200 bins were used to obtain the density distribution.



RESULTS AND DISCUSSION The density distributions for different surface separations (Figure 4) clearly show the formation of layers. However, the nonzero density between layers in the central region for H = 5.3 nm (n = 7) layers to H = 3.3 nm (n = 4) layers implies that the layered structures are not compact. Similar features have been observed by Xu and Leng28 for the fully flexibile OMCTS model. The density peaks for the contact layers for n = 7 to n = 3 are approximately similar, while the maximum peak intensity was observed for the contact layers at H = 2.11 nm. We observe a reduced peak height and a broader density distribution for n = 2, corresponding to H = 1.77 nm, indicating that the OMCTS molecules do not pack as efficiently at this surface separation, despite the increase in pore density. In cases where distinct layers are formed, the density distribution resembles those obtained with a spherical OMCTS model10 with an effective diameter of 8 Å. Since the OMCTS molecular model yields a bulk melting point between 260 and 265 K, we carried out bath simulations at lower temperatures, in order to investigate changes in pore density as the temperatures was lowered (Table 1). Simulations carried out for n = 7 (H = 5.3 nm) reveal a marginal increase in density as the temperature is lowered, resulting in a slight increase in intensity of the density peaks, as illustrated in Figure 5. The diffuse nature of the density distributions suggests that

Figure 6. Interaction energies as a function of surface separation for OMCTS confined between mica surfaces at 300 K. Continuous slope change is observed at H = 2.11 nm (n = 2 layers) which indicates the absence of a first-order transition.

functions (PCFs) for different layers based on center of mass positions using g (r ) =

1 Nl

Nl

∑ i=1

Ni(r + Δr /2, r − Δr /2) ρl (z)2πr Δr

(5)

where Ni(r + Δr/2, r − Δr/2) is the number of particles from particle i at a distance r, located in a cylindrical shell of thickness Δr. ρl(z) is the areal density of the layer containing Nl particles. The bounds of each layer are obtained from the density distributions. For some cases, we have computed the PCF using the standard definition with a spherical cutoff. Figure 7a illustrates the PCF for the contact layers as a function of surface separation. The overall features of the PCF are qualitatively similar to bulk liquid OMCTS despite the contact layers having the highest local density. We did not observe any signatures of long-range periodicity characteristic of crystalline structure even under extreme confinement. The first peak in the PCF which corresponds to the coordination in the first neighbor shell and a measure of the local packing appears at 8.7 Å, for all the surface separations, except for H = 1.77 nm where the peak is sharper and shifted to 9 Å. In contrast, the first peak occurs at 8.88 Å for bulk OMCTS at 300 K. Additionally, a weak second peak occurs for the first time at 1.83 nm for H = 1.77 nm, indicating that the in-plane order is a maximum for this two-layered situation, despite the density peak having a lower maximum when compared with H = 2.11 nm. Hence, there is a competition between layering and local packing effects at these small surface separations. With the exception of H = 1.77, the magnitude of the first peak height is only marginally higher than that observed for the bulk fluid. The in-plane PCFs for the contact layers are illustrated in Figure 7b for pore densities obtained at lower temperatures, for H = 5.23 nm (n = 7). The PCFs of the contact layers do not show any significant increase in the in-plane structure even at a temperature of 265 K which is about 5 K above the freezing point predicted by the molecular model for OMCTS. Orientational Distribution Functions. The orientation of molecules is an important characteristic of polyatomic molecules, as it is related to fluid wall interactions and local packing effects. We determine the probability distribution of the angle between the vector joining the two diagonally

Figure 5. Density distribution obtained at different temperatures for H = 5.3 nm. A weak enhancement in peak density is observed at lower temperatures of 265 K.

OMCTS remains liquid-like even at 265 K and no evidence of distinct layering, the first signatures of freezing, were observed. The fluid−fluid, fluid−wall, and total interaction energies as a function of surface separation for OMCTS confined between mica surfaces at 300 K are illustrated in Figure 6. The fluid− wall interaction increases as H is decreased, while fluid−fluid interactions increase with decreasing H. The strength of the total interaction energy increases with decreasing surface separation. The total internal energy changes continuously as the degree of confinement (H) is decreased, indicating the absence of of a first-order transition. Pair Correlation Functions. The structure of the confined fluid is examined by analyzing in-plane pair correlation F

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

Article

The Journal of Physical Chemistry B

corresponds to an angle of 45°, and the latter corresponds to a perpendicular orientation with the surface normal as observed in the snapshots for H = 1.77 nm (n = 2) in Figure 9b. At H = 1.77 nm, the distribution distinctly favors molecular orientations that are parallel to the mica surface. For H = 2.11 nm (n = 2), the distribution distinctly shifts toward molecules oriented at about 45° to the mica surface. At larger separations, n = 3 (H = 2.83 nm) and n = 7 (H = 5.3 nm), a more uniform sampling of angles between 60 and 90° is observed with a distinct peak at 45° to the surface normal. Our data shows that contact layers orient parallel to the mica surface only under extreme confinement, as seen by a distinct peak corresponding to 0° inclination with the surface normal for H = 1.77 nm, which is also revealed in the snapshot plotted in Figure 9b. The orientational distribution function was analyzed for the inner layers and contrasted with the contact layers for separations corresponding to H = 5.3 nm (7 layers) and H = 2.83 nm (3 layers), as illustrated in Figure 9c and d, respectively, at 300 K. For n = 7 layers, the inner layers uniformly sample 45−90° orientations with the surface normal and a distinct preference for the tilted configuration of 45° is observed only in the contact layers. However, for n = 3 layers, the inner layer showed similar orientational distributions as the contact layers, indicative of a templating effect at these smaller surface separations. The results indicate that molecules adopt an optimal relative configuration with the orientational degrees of freedom at a given surface separation. These structural rearrangements have a strong impact on the particle dynamics which we investigate next. Mean Squared Displacement. The dynamics of the confined liquid is studied by evaluating the mean squared displacement (MSD) in the plane parallel to the mica surface. Since the fluid is confined between two extended mica surfaces, the displacement perpendicular to the confining surface is limited by the extent of separation between the two surfaces. Hence, one cannot extract a diffusion coefficient from the long time limit of the mean squared displacment in this direction.11 Figure 10a illustrates the mean squared displacement calculated for OMCTS at 300 K at different surface separations for all the fluid layers. Two distinct regimes are captured in the MSD versus time plots. At short times, corresponding to the ballistic regime, the MSD is proportional to t2. At longer times, the mean squared displacement is proportional to t. We did not observe any significant caging effects which would result in flattening of the MSD between the ballistic and diffusive regimes.11,19 However, we observe a distinct region spanning over 3 orders in the time scale, between the ballistic and diffusive times where the MSD scales as tα with 1 < α < 2. This intermediate regime is observed for both the bulk liquid as well as confined OMCTS with a distinct gradual dynamical slowdown for the confined fluid. A similar transition behavior between ballistic and diffusive regimes has been observed for low density confined alkanes.57 The in-plane diffusion coefficient for the confined liquid is evaluated in the linear regime using the Einstein relation

Figure 7. Structural properties of OMCTS confined between mica surfaces for various surface separations and temperatures. (a) In-plane pair correlation functions for contact layers at various surface separations. The radial distribution function for the bulk liquid corresponds to 300 K. (b) In-plane pair correlation function for contact layers corresponding to H = 5.3 nm at various temperatures.

opposite Si−Si atoms in OMCTS and the unit normal to the mica surface using P(cos θ ) =

1 Nl

Nl

∑ δ(cos θi) i=1

(6)

where Nl is the number of molecules in a given layer; r = exrx + eyry; cos θi = (r·ez)/|r|; θi is the angle between the Si−Si vector and the surface normal of the ith molecule, ez, as illustrated in Figure 8. Figure 9a illustrates the orientational probability distribution for n = 2, 7, and 3 layers. We observed two dominant peaks at cos θ = 0.66 and cos θ = 0. The former

D = lim

t →∞

⟨|r(t ) − r(0)|2 ⟩ 4t

(7)

where r = exrx + eyry and the angular brackets in the above equation represent the time averages. The diffusion coefficients for various surface separations are listed in Table 2. From the values of the self-diffusion coefficient, confinement induced

Figure 8. Molecular representation of OMCTS, indicating the vector r used to define the orientational probability distribution of OMCTS confined between mica surfaces. Yellow, silicon atoms; red, oxygen atoms; blue, methyl groups. G

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

Article

The Journal of Physical Chemistry B

Figure 9. (a) Orientational probability distribution function for OMCTS confined between mica surfaces for various surface separations at 300 K. (b) Snapshot of a contact layer for H = 1.77 nm (2 layers) at 300 K. Triangles represent molecules inclined at about 45° with the surface normal, and circles represent molecules inclined at 90° to the surface normal. (c) Layerwise analysis of the probability distribution function for H = 5.3 nm (n = 7 layers) at 300 K. (d) Layerwise analysis of the probability distribution function for H = 2.83 nm (n = 3 layers) at 300 K.

Table 2. Diffusion Coefficients for Various Surface Separations of the Confined System Estimated for All the Fluid Layers as Well as for the Inner Layers at 300 K D × 10−5 cm2 s−1 H (nm)

T (K)

all layers

inner layers

5.3 4.12 3.3 2.83 2.11 1.77 5.3 5.3

300 300 300 300 300 300 275 265

0.299 0.248 0.190 0.114 0.067 0.025 0.122 0.04

0.350 0.419 0.423 0.698

diffusion slowdown is observed for lower surface separations and the diffusion coefficient steadily decreases with decreasing H. The diffusion coefficient for H = 5.3 nm (n = 7 layers) is comparable to that of bulk OMCTS liquid. The diffusion coefficient reported by Xu and Leng28 for a three-layered system is 0.147 × 10−5 cm2 s−1, which is comparable with the diffusion coefficient obtained in this study for the same number of layers, for H = 2.83 nm, as seen in Table 2. Confinement induced diffusion slowdown, but liquid-like behavior was also demonstrated for the OMCTS−mica system by Matsubara et al.29 where the diffusion slowdown was explained using the concept of activated diffusion theory. Additionally, we observe a nonmonotonic dependence of self-diffusion coefficient on surface separations, a phenomenon well documented for confined monatomic fluids.4,58 At this juncture, it is useful to point out that we do not observe any signatures of freezing.

Figure 10. (a) Mean squared displacements for various surface separations, H, at 300 K. (b) Self-diffusion coefficient as a function of packing fraction. The solid line represents a power law fit of the form D(ϕc − ϕ)γ, where ϕc = 0.92 and γ = 1.852. H

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

Article

The Journal of Physical Chemistry B Even at the smallest surface separation of H = 1.77 nm where n = 2, we observe that the fluid is able to sample the diffusive regime. While investigating dynamical arrest of a fluid approaching the glass transition temperature using mode coupling analysis, one typically studies the divergence in diffusivity as the temperature is lowered.11,19 However, in our study, since we are interested in the diffusion slowdown with increasing confinement at a fixed temperature, we investigate the divergence in diffusivity with the fluid density as reflected in the packing fraction.22,59 Since each OMCTS molecule is made up of eight united atom methyl groups, four oxygen and silicon atoms, the volume of the molecule (0.318 nm3) was obtained by summing the volume of the individual atoms. The packing fraction is evaluated as a ratio between the volume occupied by the molecules to the volume of the accessible space between the two surfaces. In order to compute the accessible volume, for a given area of the surfaces, the layer thickness is obtained from the distance between two minima that define the layers in the density distributions for a given surface separation H. Figure 10b illustrates the dependence of diffusion coefficient on the packing fraction based on all the fluid layers. A power law of the form D ∼ (ϕc − ϕ)γ as predicted by MCT scalings59 is seen to fit the data reasonably well. The critical packing fraction is ϕc = 0.92, and the exponent γ is 1.852. We also analyzed the inner layers separately and obtained a similar value of ϕc = 0.92 with γ = 1.921. We next investigate the trends in the MSD at lower temperatures, which correspond to temperatures in the range 5−35 K above the determined melting point of the model OMCTS used in this study. The mean squared displacements for H = 5.3 nm (n = 7 layers) at various temperatures are plotted in Figure 11. As the temperature is lowered, a mild

or density−density correlation, for various degrees of confinement at 300 K are plotted in Figure 12 evaluated at the

Figure 12. (a) Self-intermediate scattering function for all the fluid layers at 300 K, with wavenumber, k = 0.722 Å−1. The corresponding Kohlrausch−Williams−Watts (KWW) parametrization is represented by the fine dashed line. (b) Self-intermediate scattering function for the inner fluid layers at 300 K, corresponding to the wavenumber, k = 0.722 Å−1. The corresponding Kohlrausch−Williams−Watts (KWW) parametrization is represented by the fine dashed line.

wavenumber corresponding to the first peak in the in-plane PCF for the contact layers. At short times, the curves collapse onto a single curve, and at longer times, we observe a one step relaxation following a stretched exponential decay. This trend is observed until n = 3 layers. However, we see a weak two-step relaxation at H = 2.11 nm (n = 2 layers) typical of glassy dynamics. The mean squared displacement for the corresponding surface separations (Figure 10a) has a weak and extended subdiffusive regime. As the surface separation is decreased, we observe a progressive delay in the onset of the α relaxation of the Fs(k, t), resulting in a distinct increase in the relaxation time. In order to understand the effect of contact layers on the dynamics, the self-intermediate scattering function was estimated separately for the inner layers, as illustrated in Figure 12b. At larger surface separations, the relaxation of the inner layers is distinctly faster when compared with the relaxation obtained when the contact layers are included (Figure 12a). A layerwise analysis of the Fs(k, t) for H = 2.83 nm (n = 3) is illustrated in Figure 13. We see that the inner layers relax about 10 times faster than the contact layers, indicating that the contact layers play an important role in determining the overall relaxation characteristics of the confined fluid. The different relaxation regimes in the Fs(k, t) can be classified as ballistic (for time 230 K, the confined crystal evolves to its final disordered state within about 0.5 ns from the beginning of the MD simulations. The disordering is more rapid when the surface separation is decreased. For example, at H = 3.95 nm, the initial crystalline state rapidly disorders within the first 50 ps. In-Plane Pair Correlation Function. To further confirm the evolution of the confined fluid structure from the initial crystalline configuration, the in-plane pair correlation functions based on the OMCTS center of mass for the contact and inner layers were evaluated (Figure 19a and b, respectively). There is a distinct absence of long-range ordering in the contact layers for all the temperatures investigated. We observe shouldering in the first peak of the PCF as a consequence of preferential orientations induced by the confining surfaces. For T > 230 K, the disordering in the contact layer percolates to the inner layers which also remain disordered. However, we do not observe shoulders in the peaks for the innermost layer. At 220 K, we clearly observe long-range translational ordering indicative of a crystalline innermost layer. The first inner layers also exhibited long-range translational ordering (not illustrated). This suggests that the crystallization phenomenon for the confined fluid is nonuniform and the contact layers are likely to crystallize at a much lower temperature than observed for the innermost layers. Such heterogenity in crystallization for confined fluids has been observed earlier.10,15 Experimental investigations for confined solids also suggest that the contact layers can be liquid-like while the inner layers can remain frozen.65 These studies reaffirm that confinement tends to reorient the molecules breaking the crystal symmetry even at temperatures below the melting point of the bulk crystal. We point out that, as the system was cooled, voids began to appear in the final configuration at 230 K, indicating that a higher density is required to stabilize and pack the crystal under confinement at the lower temperatures. These freezing transitions also depend on the surface separations which alter the density of the confined fluid. Free Energy Analysis: Two-Phase Thermodynamic Method (2PT). The equilibrated configurations obtained from the 5 ns NVT simulations for the confined five-layer OMCTS crystal were further simulated for 100 ps using the Nosé− Hoover thermostat with a coupling constant equal to 0.5 ps and a time step equal to 2 fs. Data from the last 40 ps sampled every 4 fs were used for the 2PT analysis. The MDWiZ platform was used to convert the structure file from GROMACS to LAMMPS format66 for use with the 2PT program. We do not present the details of the underlying theory behind the 2PT

2⎥

|q lm(i)|

T (K)

(13)

Since the bulk crystal is a tetragonal crystal with space group symmetry P42/n,54 Q4 = 1. The variation of Q4 with temperature is illustrated in Figure 18 and Table 4. For temperatures above 230 K, disordered structure is observed, as reflected in the low values of Q4. As the system is cooled below 230 K, a sharp increase in Q4 is observed, indicative of an ordering transition in OMCTS. This transition which occurs

Figure 18. Temperature dependence of the bond order parameter. Data from MD simulations carried out with a confined crystal as the initial configuration. An order−disorder transition occurs at 230 K for the confined system. Error bars are small (Table 4) and hence not indicated in the figure. L

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

Article

The Journal of Physical Chemistry B

of 1 indicates that the system behaves like a gas without any solid component and a value of zero indicates that the system is a solid without a gas component. Low values of fluidicity were observed for temperatures less than and equal to 230 K which suggests a nondiffusive behavior akin to a solid phase. Above 230 K, a sharp rise in fluidicity is observed which indicates melting of the confined crystal and the values remain almost constant, independent of temperature. The dependence of free energy on temperature is illustrated in Figure 20a. We observe a discontinuity in the slope at about 230 K which is consistent with the sharp change observed in the bond order parameter with temperature (Figure 18) signaling an order−disorder transition when the temperature is lowered. We further investigated the variation of potential energy (fluid−fluid and fluid−wall) with temperature for the confined crystal, illustrated in Figure 20b. We observe a discontinuity in slope at 230 K which is a signature of a firstorder transition without accounting for possible finite size effects. Our simulations reveal that OMCTS undergoes a freezing point depression by about 30 K under confinement. The phenomenon of freezing point depression is well documented in the confined fluids literature, and several experiments and simulations have shed light on this phenomenon. The melting point depression of about 30 K observed for OMCTS crystal when confined between mica surfaces (H = 4.12 nm, n = 5 layers) compares well with previously reported differential scanning calorimetry data67 where a depression in melting point of about 25 K was observed for small organic molecules such as benzene and cdecalin and about 45 K for cyclohexane when confined in controlled pore glasses with a pore diameter of 4 nm. The elevation or depression in freezing point is governed by several factors such as fluid−wall strength, fluid−fluid interactions, wall density, surface to volume ratio of the confined solid, and separation distance between the confining surfaces. Further, the frozen phase can have several topological variants such as the inner layer frozen situation as observed in our simulations. These aspects as well as the limitations in experimentaly detecting these variations have been reviewed by AlbaSimionesco et al.68 The ratio of fluid−wall to fluid−fluid interaction was less than 1 for this system, and the hypothesis that the freezing point gets suppressed for low ratios of fluid− wall to fluid−fluid interaction is in agreement with the confined fluid phase diagram proposed by Radhakrishnan etal.15 The influence of all of these factors on the melting point of the confined crystal was not investigated and is beyond the scope of this study. From the free energy analysis, confined OMCTS with the parameters used in our study undergoes a depression in freezing point by about 30 K from the bulk OMCTS crystal. The symmetry of the confined frozen phase observed at low

Figure 19. In-plane pair correlation functions for the confined fluid where simulations were carried out with a five-layered confined crystal as the initial configuration (H = 4.32 nm). The pair correlation functions depict the structure for the equilibrated 5 ns simulations. (a) Contact layers and (b) contact layers and first inner layers excluded.

free energy method, and the reader is referred to several papers with applications of the 2PT method to both bulk39,40,43 and confined fluids.41,42,44−47 The entropy S, fluidicity factor f obtained from the 2PT method, potential energy from simulations (PE), and hence the Helmholtz free energy (A) calculated on a per molecule basis are tabulated in Table 5. The fluidicity factor f defined as the fraction of gas-like component in the overall system is calculated from the ratio of the selfdiffusion coefficient of the confined fluid to that of the hard sphere diffusivity under the corresponding thermodynamic conditions. Since the self-diffusion has both translational and rotational components, the fluidicity from each of these contributions can be evaluated in the 2PT method. A value

Table 5. Thermodynamic Quantities Obtained from Confined Crystal Simulations at Various Temperatures T (K)

S (kJ mol−1)

200 220 230 240 250 260 270 280

0.264 0.289 0.302 0.319 0.331 0.339 0.348 0.357

PE (kJ mol−1) 37.874 47.596 52.583 58.734 63.153 67.389 71.173 75.289

± ± ± ± ± ± ± ±

A (kJ mol−1) −14.926 −15.984 −16.877 −17.826 −19.597 −20.751 −22.787 −24.671

0.129 0.290 0.150 0.568 0.106 0.019 0.048 0.078 M

± ± ± ± ± ± ± ±

0.158 0.064 0.083 0.301 0.063 0.123 0.148 0.157

f trans

f rot

0.0282 0.0391 0.0549 0.0726 0.0917 0.0909 0.0881 0.0904

0.0241 0.0457 0.0615 0.0895 0.0998 0.1025 0.1023 0.107 DOI: 10.1021/acs.jpcb.5b12759 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

findings, we computed the free energy of confined OMCTS using the two-phase thermodynamic (2PT) model as a function of temperature and monitored the order parameter, free energy, and internal energy as a function of temperature for a specific surface separation. We carry out simulations at temperatures corresponding to an elevation ranging from 5 to 40° above the melting point of 260 K predicted by the molecular model of OMCTS. Our simulations reveal that OMCTS undergoes increasing stuctural arrest akin to a fluid approaching a glass transition as the degree of confinement increases from n = 7 layers to n = 2 layers. Even under extreme confinement when two fluid layers are present, the self-intermediate scattering function shows a weak two-step relaxation characteristic of glass-like dynamics of a fluid above the glass transition temperature. At this extreme confinement, the scattering function is seen to relax completely within 10 ns. At n = 7, a distinct two-step relaxation in the selfintermediate scattering function is observed with a well-defined β relaxation regime as the temperature is reduced. Even at 5° above the melting point, we did not observe freezing of the confined fluid and a distinct glass-like relaxation signature is observed at this surface separation. The lack of freezing can be partly explained due to the inherent anisotropy in the OMCTS molecules which pre-empts a freezing transition upon confinement. Hence, long-range peridocity in the in-plane pair distribution functions is distinctly absent even for the fluid layers directly in contact with the mica surface. Under extreme confinement, the molecules tend to preferentially orient in a direction parallel to the confining surface. Our observation that confined OMCTS undergoes a transition similar to vitrification or a glass transition is borne out by several studies of small molecules such as benzene and phenol which do not crystallize upon confinement.69,70 Thus, crystallization under confinement is prevented in these small molecules due to the presence of a surface which disrupts the preferred packing symmetry to form a crystalline phase. These trends are consistent with our observations in this study for OMCTS which shows a gradual slowing down upon confinement. We characterize the divergence in the diffusion coefficient and relaxation time as a function of the packing fraction which is an appropriate reflection of the density enhancement with increasing confinement. Both the self-diffusion coefficient and relaxation time in the late α relaxation regime show a power law divergence with packing fraction. A critical value of packing fraction has also been extracted to define the glass transition line based on which a phase diagram was constructed in this study. On the basis of different interpretations from the experimental surface force experminents and our results from the dynamical characterization of confined OMCTS, we conclude that OMCTS undergoes a gradual glass-like transition upon confinement between mica surfaces. We digress briefly to compare our results with previous simulation studies which are related to freezing of confined fluids in general and studies specifically related to OMCTS confined between mica surfaces. Monte Carlo simulations and Landau free energy calculations with monatomic Lennard-Jones fluids13 show that the freezing point can either be elevated or depressed with respect to the bulk freezing point. This shift in the freezing transition temperature is correlated with the strength of the fluid−wall interaction. Other studies with monatomic systems have shown that confined fluids not only freeze but undergo solid−solid transitions similar to the scenarios observed with confined colloids.10,71,72 More recently,

Figure 20. (a) Dependence of Helmholtz free energy with temperature. The continuous line represents a linear fit, and circles represent data points. Error bars are indicated by black solid lines. The first-order transition is characterized by the equality of Helmholtz free energy which is observed at about 235 K from the intersection of the linear fits. (b) Variation of potential energy (fluid−fluid and fluid−wall) with temperature. The continuous line represents linear fits, and data points are marked by squares. Error bars indicated by black solid lines lie within the data points except at the transition. Discontinuity in slope is observed between 230 and 240 K which suggests the presence of a first-order transition.

temperatures is complicated with a surface layer oriented preferentially toward the mica surface with long-range order observed only in the inner layers at lower temperatures.



SUMMARY AND CONCLUSIONS We have carried out molecular dynamics simulations in order to understand the fluid structure and dynamics of OMCTS confined between mica surfaces. We use a recently developed eight-site molecular model for OMCTS which provides a realistic model for studying layering and orientational effects in contrast to the spherical model used in the past. In addition to examining the density distributions, pair-correlation functions, and orientational distributions which yield structural information about the confined fluid, we also compute the mean squared displacements and self-intermediate scattering functions to study the dynamics and relaxation behavior of the confined fluid. Using these appropriate dynamical signatures, we comment on whether the confined fluid undergoes a freezing transition or whether the transition is akin to that of a fluid approaching a glass transition. To further support our N

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

Article

The Journal of Physical Chemistry B

probe the state of the confined fluid without the application of an external shear and hence describe the fluid state in the low frequency or long time limit.

we have illustrated that even a monatomic fluid is capable of depicting a glass-like transition under appropriate conditions.11 Our earlier study with a spherical model for OMCTS confined between mica surfaces revealed a distinct tendency of the confined fluid to freeze upon increasing confinement. This is consistent with freezing scenarios encountered for confined monatomic fluids. Recently, several molecular dynamics simulation studies for OMCTS have been based on a molecular model for OMCTS, a departure from the spherical model assumption. Given the significant shape anisotropy of OMCTS,38 it is clear that a spherical or monatomic model used previously is an unrealistic representation for OMCTS and at best represents freezing trends typically encountered in monatomic systems. While using any molecular model, the interpretation of results is then dependent on the accuracy of the molecular model itself. In this regard, we point out that the present molecular model, which is parametrized to predict bulk liquid properties of OMCTS, underestimates the freezing point by 30−35 K, as illustrated in this work. A molecular model parametrized to predict the freezing point of OMCTS correctly would be an improvement over the existing model. However, it is unlikely to alter the main conclusions. We carried out a detailed free energy analysis using the 2PT method for an n = 5 layered system, from 280 K to temperatures about 60 K below the melting point of the OMCTS crystal. A distinct discontinuity in the Helmholtz free energy and potential energy was observed at 230 K, indicative of a first-order transition. Consistent with this transition, a sharp rise in the local bond order parameter, Q4, was observed for temperatures below 230 K. At temperatures below 230 K where the confined fluid begins to order, the in-plane PCF for inner layers revealed the onset of long-range translational ordering. However, the contact layers were found to remain disordered at all temperatures with the molecules preferring a predominantly surface parallel orientation at lower temperatures. The pair correlation functions and the orientational bond order parameter are related to the Landau free energy,64 leading us to conclude that confinement of a molecular crystal such as OMCTS disrupts long-range ordering even at temperatures below the melting point of bulk OMCTS. Hence, the crystal is unstable under confinement up to temperatures about 30 K below the melting point of the bulk crystal. Thus, a distinct freezing point depression was observed for OMCTS under confinement. A similar observation of freezing point depression and glassy dynamics have been documented for moelcular fluids such as benzene and toluene under confinement.60 Although we have computed the free energy for a specific surface separation, we expect to see similar trends at smaller surface separations. On the basis of both the dynamic and thermodynamic analysis carried out in this study, our findings support the interpretation that OMCTS undergoes a glass-like transition when confined between mica surfaces. Additionally, free energy calculations reveal a freezing point depression of about 30 K for a specific confinement of five layers (H = 4.12 nm). The structurally arrested low temperature solid phase consists of frozen inner layers with the OMCTS molecules having a surface parallel orientation at the mica surface. Freezing point depression in nanopores has been well documented using experiments as well as simulations for small organic molecules68 confined in controlled pore glasses,67 and our study shows that OMCTS follows a similar behavior. We finally mention that equilibrium molecular dynamics simulations used in this work



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +91 (0)80 22932769. Fax: +91 (0)80 23608121. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The authors thank Ateeque Malani for providing the coordinates of the mica surface used in this study and Uzi Landmann for helpful discussions related to the melting point determination. We thank Sudeep Punnathanam and V. Kumaran for several insightful discussions related to this work. We would like to thank Prabal K. Maiti for sharing the 2PT code as well as for the assistance in implementing the free energy calculations reported in this manuscript. We also acknowledge Shiang-Tai Lin for discussions on the 2PT method. We thank Riccardo Baron for providing access to the MDWiZ platform for converting GROMACS structure files to LAMMPS format. The authors also thank the Supercomputer Education and Research Centre (SERC) for computational resources as well as Department of Science and Technnology (DST) India for support.



REFERENCES

(1) Rhykerd, C. L.; Schoen, M.; Diestler, D. J.; Cushman, J. H. Epitaxy in Simple Classical Fluids in Micropores and Near-solid Surfaces. Nature 1987, 330, 461. (2) Pozhar, L.; Gubbins, K. Transport Theory of Dense, Strongly Inhomogeneous Fluids. J. Chem. Phys. 1993, 99, 8970−8996. (3) Somers, S. A.; Davis, H. T. Microscopic Dynamics of Fluids Confined Between Smooth and Atomically Structured Solid Surfaces. J. Chem. Phys. 1992, 96, 5389−5407. (4) Gao, J.; Luedtke, W. D.; Landman, U. Layering Transitions and Dynamics of Confined Liquid Films. Phys. Rev. Lett. 1997, 79, 705− 708. (5) Granick, S. Motions and Relaxations of Confined Liquids. Science 1991, 253, 1374−1379. (6) Klein, J.; Kumacheva, E. Confinement Induced Phase Transitions in Simple Liquids. Science 1995, 269, 816−819. (7) Klein, J.; Kumacheva, E. Simple Liquids Confined to Molecularly Thin Layers.I.Confinement-Induced Liquid to Solid Phase Transitions. J. Chem. Phys. 1998, 108, 6996. (8) Kumacheva, E.; Klein, J. Simple Liquids Confined to Molecularly Thin Layers.II. Shear and Frictional Behavior of Solidified Films. J. Chem. Phys. 1998, 108, 7010. (9) Ghatak, C.; Ayappa, K. G. Solid Solid Transformations in a Confined Soft Sphere Fluid. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 64, 051507. (10) Ayappa, K. G.; Mishra, R. K. Freezing of Fluids Confined Between Mica Surfaces. J. Phys. Chem. B 2007, 111, 14299−14310. (11) Krishnan, S. H.; Ayappa, K. G. Glassy Dynamics in a Confined Monatomic Fluid. Phys. Rev. E 2012, 86, 011504. (12) Bock, H.; Gubbins, K. E.; Ayappa, K. G. Solid/Solid Phase Transitions in Confined Thin Films: A Zero Temperature Approach. J. Chem. Phys. 2005, 122, 094709. (13) Radhakrishnan, R.; Gubbins, K. E. Free Energy Studies of Freezing in Slit Pores: An Order-Parameter Approach Using Monte Carlo Simulation. Mol. Phys. 1999, 96, 1249−1267. (14) Sliwinska-Bartkowiak, M.; Gras, J.; Sikorski, R.; Radhakrishnan, R.; Gelb, L.; Gubbins, K. E. Phase Transitions in Pores: Experimental O

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

Article

The Journal of Physical Chemistry B and Simulation Studies of Melting and Freezing. Langmuir 1999, 15, 6060−6069. (15) Radhakrishnan, R.; Gubbins, K. E.; Sliwinska-Bartkowiak, M. Effect of the Fluid-Wall Interaction on Freezing of Confined Fluids: Toward the Development of a Global Phase Diagram. J. Chem. Phys. 2000, 112, 11048−11057. (16) Miyahara, M.; Gubbins, K. E. Freezing/Melting Phenomena for Lennard-Jones Methane in Slit pores: A Monte Carlo Study. J. Chem. Phys. 1997, 106, 2865−2880. (17) Camara, L.; Bresme, F. Molecular Dynamics Simulations of Crystallization Under Confinement at Triple Point Conditions. J. Chem. Phys. 2003, 119, 2792−2800. (18) Sałamacha, L.; Patrykiejew, A.; Sokołowski, S.; Binder, K. The Structure of Fluids Confined in Crystalline Slitlike Nanoscopic Pores. J. Chem. Phys. 2005, 122, 074703. (19) Gallo, P.; Rovere, M.; Spohr, E. Glass Transition and Layering Effects in Confined Water: A Computer Simulation Study. J. Chem. Phys. 2000, 113, 11324−11335. (20) Binder, K.; Baschnagel, J.; Paul, W. Glass Transition of Polymer melts: Test of Theoretical Concepts by Computer Simulation. Prog. Polym. Sci. 2003, 28, 115−172. (21) Khare, R.; de Pablo, J.; Yethiraj, A. Molecular Simulation and Continuum Mechanics Investigation of Viscoelastic Properties of Fluids Confined to Molecularly Thin Films. J. Chem. Phys. 2001, 114, 7593−7601. (22) Mandal, S.; Lang, S.; Gross, M.; Oettel, M.; Raabe, D.; Franosch, T.; Varnik, F. Multiple Reentrant Glass Transitions in confined HardSphere Glasses. Nat. Commun. 2014, 5, 1−5. (23) Mizukami, M.; Kusakabe, K.; Kurihara, K. Progress in Colloid and Polymer Science; Springer: Berlin Heidelberg, 2004; pp 105−108. (24) Becker, T.; Mugele, F. Nanofluidics: Molecularly Thin Lubricant Layers Under Confinement. Mol. Simul. 2005, 31, 489−494. (25) Bureau, L. Nonlinear Rheology of a Nanoconfined Simple Fluid. Phys. Rev. Lett. 2010, 104, 218302−218305. (26) Klein, J.; Kumacheva, E. Simple Liquids Confined to Molecularly Thin Layers. I. Confinement Induced Liquid to Solid Phase Transitions. J. Chem. Phys. 1998, 108, 6996−7009. (27) Zhu, Y.; Granick, S. No-Slip Boundary Condition Switches to Partial Slip When Fluid Contains Surfactant. Langmuir 2002, 18, 10058−10063. (28) Xu, R.-G.; Leng, Y. Solvation Force Simulations in Atomic Force Microscopy. J. Chem. Phys. 2014, 140, 214702. (29) Matsubara, H.; Pichierri, F.; Kurihara, K. Mechanism of Diffusion Slowdown in Confined Liquids. Phys. Rev. Lett. 2012, 109, 197801−5. (30) Matsubara, H.; Pichierri, F.; Kurihara, K. Unraveling the Properties of Octamethylcyclotetrasiloxane Under Nanoscale Confinement: Atomistic View of the Liquidlike State from Molecular Dynamics Simulation. J. Chem. Phys. 2011, 134, 044536. (31) O’Shea, S.; Welland, M.; Pethica, J. Atomic Force Microscopy of Local Compliance at Solid−Liquid Interfaces. Chem. Phys. Lett. 1994, 223, 336−340. ́ (32) Sterczyńska, A.; Deryło-Marczewska, A.; Sliwiń ska-Bartkowiak, M.; Piotrowska, J.; Jarek, M.; Domin, K. Phase Transitions of Octamethylcyclotetrasiloxane Confined Inside Aluminosilicate and Silicate Nanoporous Matrices. J. Therm. Anal. Calorim. 2014, 118, 263−276. (33) Demirel, L. A.; Granick, S. Glasslike transition of a confined simple fluid. Phys. Rev. Lett. 1996, 77, 2261. (34) Lin, Z.; Granick, S. Platinum Nanoparticles at Mica Surfaces. Langmuir 2003, 19, 7061−7070. (35) Israelachvili, J. N.; Alcantar, N. A.; Maeda, N.; Mates, T. E.; Ruths, M. Preparing Contamination-free Mica Substrates for Surface Characterization, Force Measurements, and Imaging. Langmuir 2004, 20, 3616−3622. PMID: 15875391. (36) Perkin, S.; Chai, L.; Kampf, N.; Raviv, U.; Briscoe, W.; Dunlop, I.; Titmuss, S.; Seo, M.; Kumacheva, E.; Klein, J. Forces Between Mica Surfaces, Prepared in Different Ways, Across Aqueous and Non-

aqueous Liquids Confined to Molecularly Thin Films. Langmuir 2006, 22, 6142−6152. (37) Heuberger, M. The Extended Surface Forces Apparatus. Part I. Fast Spectral Correlation Interferometry. Rev. Sci. Instrum. 2001, 72, 1700−1707. (38) Matsubara, H.; Pichierri, F.; Kurihara, K. Design of a Versatile Force Field for the Large Scale Molecular Simulation of Solid and Liquid OMCTS. J. J. Chem. Theory Comput. 2010, 6, 1334−1340. (39) Lin, S.-T.; Blanco, M.; Goddard, W. A., III The Two-Phase Model for Calculating Thermodynamic Properties of Liquids from Molecular Dynamics: Validation for the Phase Diagram of LennardJones Fluids. J. Chem. Phys. 2003, 119, 11792−11805. (40) Lin, S.-T.; Maiti, P. K.; Goddard, W. A., III Two-phase Thermodynamic Model for Efficient and Accurate Absolute Entropy of Water from Molecular Dynamics Simulations. J. Phys. Chem. B 2010, 114, 8191−8198. (41) Lin, S.-T.; Maiti, P. K.; Goddard, W. A. Dynamics and Thermodynamics of Water in PAMAM Dendrimers at Subnanosecond Time Scales. J. Phys. Chem. B 2005, 109, 8663−8672. (42) Kumar, H.; Dasgupta, C.; Maiti, P. K. Driving Force of Water entry into Hydrophobic Channels of Carbon Nanotubes: Entropy or Energy? Mol. Simul. 2015, 41, 504−511. (43) Huang, S.-N.; Pascal, T. A.; William, A.; Goddard, I.; Maiti, P. K.; Lin, S.-T. Absolute Entropy and Energy of Carbon Dioxide Using the Two-Phase Thermodynamic Model. J. Chem. Theory Comput. 2011, 7, 1893−1901. (44) Mukherjee, B.; Maiti, P. K.; Dasgupta, C.; Sood, A. Strong correlations and Fickian Water Diffusion in Narrow Carbon Nanotubes. J. Chem. Phys. 2007, 126, 124704. (45) Kumar, H.; Dasgupta, C.; Maiti, P. K. Structure, Dynamics and Thermodynamics of Single-file Water Under Confinement: Effects of Polarizability of Water Molecules. RSC Adv. 2015, 5, 1893−1901. (46) Debnath, A.; Ayappa, K.; Maiti, P. K. Simulation of Influence of Bilayer Melting on Dynamics and Thermodynamics of Interfacial Water. Phys. Rev. Lett. 2013, 110, 018303. (47) Jana, B.; Pal, S.; Maiti, P. K.; Lin, S.-T.; Hynes, J. T.; Bagchi, B. Entropy of Water in the Hydration Layer of Major and Minor Grooves of DNA. J. Phys. Chem. B 2006, 110, 19611−19618. (48) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: California, USA, 2002. (49) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157− 2164. (50) Fischer, J.; Weiss, A. Transport properties of liquids. V. Self Diffusion, Viscosity, and Mass Density of Ellipsoidal Shaped Molecules in the Pure Liquid Phase. Ber. Bunsenges. Phys. Chem. 1986, 90, 896− 905. (51) García Fernández, R.; Abascal, J. L.; Vega, C. The Melting Point of Ice Ih for Common Water Models Calculated from Direct Coexistence of the Solid-Liquid Interface. J. Chem. Phys. 2006, 124, 1− 11. (52) Landman, U.; Luedtke, W. D.; Ribarsky, M. W.; Barnett, R. N.; Cleveland, C. L. Molecular Dynamics Simulations of Epitaxial Crystal Growth from the Melt. I. Si(100). Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 4637−4646. (53) Luedtke, W. D.; Landman, U.; Ribarsky, M. W.; Barnett, R. N.; Cleveland, C. L. Molecular Dynamics Simulations of Epitaxial Crystal Growth from the Melt. II. Si(111). Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 4647−4655. (54) Steinfink, H.; Post, B.; Fankuchen, I. The Crystal Structure of Octamethyl Cyclotetrasiloxane. Acta Crystallogr. 1955, 8, 420−424. (55) Hansen, J.-P.; Verlet, L. Phase Transitions of the Lennard-Jones System. Phys. Rev. 1969, 184, 151−161. (56) Malani, A.; Ayappa, K. G. Adsorption Isotherms of Water on Mica: Redistribution and Film Growth. J. Phys. Chem. B 2009, 113, 1058−1067. P

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

Article

The Journal of Physical Chemistry B (57) Cruz, F.; Müller, E. Behavior of Ethylene and Ethane Within Single-Walled Carbon Nanotubes, 2: Dynamical properties. Adsorption 2009, 15, 13−22. (58) Magda, J. J.; Tirrell, M.; Davis, H. T. Molecular Dynamics of Narrow, Liquid Filled Pores. J. Chem. Phys. 1985, 83, 1888. (59) Gotze, W. Complex Dynamics of Glass-Forming Liquids: A ModeCoupling Theory; Oxford University Press: Oxford, 2008. (60) Alba-Simionesco, C.; Dosseh, G.; Dumont, E.; Frick, B.; Geil, B.; Morineau, D.; Teboul, V.; Xia, Y. Confinement of Molecular Liquids: Consequences on Thermodynamic, Static and Dynamical Properties of Benzene and Toluene. Eur. Phys. J. E: Soft Matter Biol. Phys. 2003, 12, 19−28. (61) Kob, W.; Andersen, H. C. Testing Mode-Coupling Theory for a Supercooled Binary Lennard-Jones Mixture I: The Van Hove Correlation Function. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1995, 51, 4626−4641. (62) Kob, W.; Andersen, H. C. Scaling Behavior in the β-Relaxation Regime of a Supercooled Lennard-Jones Mixture. Phys. Rev. Lett. 1994, 73, 1376−1379. (63) Donth, E. The Glass Transition: Relaxation Dynamics in Liquids and Disordered Materials; Springer Series in Materials Science; Springer: Berlin, Heidelberg, 2010. (64) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. BondOrientational Order in Liquids and Glasses. Phys. Rev. B: Condens. Matter Mater. Phys. 1983, 28, 784−805. (65) Dosseh, G.; Xia, Y.; Alba-Simionesco, C. Cyclohexane and Benzene Confined in MCM-41 and SBA-15: Confinement Effects on Freezing and Melting. J. Phys. Chem. B 2003, 107, 6445−6453. (66) Rusu, V. H.; Horta, V. A.; Horta, B. A.; Lins, R. D.; Baron, R. MDWiZ: A Platform for the Automated Translation of Molecular Dynamics Simulations. J. Mol. Graphics Modell. 2014, 48, 80−86. (67) Jackson, C. L.; McKenna, G. B. The Melting Behavior of Organic Materials Confined in Porous Solids. J. Chem. Phys. 1990, 93, 9002−9011. (68) Alba-Simionesco, C.; Coasne, B.; Dosseh, G.; Dudziak, G.; Gubbins, K. E.; Radhakrishnan, R.; Sliwinska-Bartkowiak, M. Effects of Confinement on Freezing and Melting. J. Phys.: Condens. Matter 2006, 18, R15−R68. (69) Johnston, K.; Harmandaris, V. Properties of Benzene Confined between Two Au(111) Surfaces Using a Combined Density Functional Theory and Classical Molecular Dynamics Approach. J. Phys. Chem. C 2011, 115, 14707−14717. (70) Audonnet, F.; Brodie-Linder, N.; Morineau, D.; Frick, B.; AlbaSimionesco, C. J. Non-Cryst. Solids 2015, 407, 262−269. Seventh IDMRCS: Relaxation in Complex Systems. (71) Ghatak, C.; Ayappa, K. G. Solid-Solid Transitions in Slit Nanopores. Colloids Surf., A 2002, 205, 111−117. (72) Vishnyakov, A.; Neimark, A. V. Specifics of Freezing of LennardJones Fluid Confined to Molecularly Thin Layers. J. Chem. Phys. 2003, 118, 7585−7598.

Q

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