Structure and Mass Transport Characteristics at the Intrinsic Liquid

Article pubs.acs.org/JPCB

Structure and Mass Transport Characteristics at the Intrinsic Liquid− Vapor Interfaces of Alkanes Hari Krishna Chilukoti,* Gota Kikugawa, and Taku Ohara Institute of Fluid Science, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan ABSTRACT: In this paper, an instantaneous interface definition has been used to study the intrinsic structure and self-diffusion coefficient in the vicinity of the liquid−vapor interfaces of decane and tetracosane at three different temperatures using molecular dynamics simulations, and the results have been compared with those obtained on the basis of the conventional Gibbs dividing surface (time- and space-averaged interface). The alkane molecules were modeled using the united atom NERD force field. Partial layered structures of alkane molecules at the liquid− vapor interface are observed as a pinned structure of alkane liquids based on the intrinsic interface. This kind of characteristic has not been observed in the density profiles obtained based on the Gibbs dividing surface. By examining the orientation order parameter and radius of gyration of the alkane molecules, it was observed that the alkane molecules were preferentially oriented to be more parallel to the intrinsic interface than to the Gibbs dividing surface, and the shape of the alkane molecules is slightly changed in the vicinity of the liquid−vapor interfaces. The self-diffusion coefficient parallel to the intrinsic interface was examined using the Green− Kubo relation, where the projection of the velocity in the parallel direction to the local intrinsic interface is used in the velocity correlation function. It was found that the self-diffusion coefficient in the direction parallel to the intrinsic interface changes as the position approaches the interface in a more obvious manner as compared with the self-diffusion coefficient obtained with respect to the Gibbs dividing surface. These results suggest that the use of an instantaneous interface definition allowed us to capture sharp variations in transport properties which are originating due to steeper structure at the liquid−vapor interfaces. influenced by the averaging effects of thermal fluctuations, they will be surely different from those obtained based on the timeand space-averaged Gibbs dividing surface. This kind of molecular-level understanding of physical and structural quantities gives a good insight into physics and chemistry taking place at the liquid−vapor interfaces. For example, the chemical reactivity of a liquid−vapor interface on a local scale is mainly controlled by the functional groups offered to the vapor phase.14 Knowledge of equilibrium and nonequilibrium properties at the interfaces which are obtained by considering the capillary fluctuations have applications in capillary drying and the kinetics of self-assembly.15 To consider the molecular-scale surface fluctuations while calculating properties, the main difficulty is to define the interface position which is changing with time and space by considering the movement of molecules between liquid and vapor phases. Recently, a few researchers have defined the intrinsic interface and examined the intrinsic structure and physical properties at the liquid−liquid and liquid−vapor interfaces.16−24 Due to its importance in chemical and biological processes, mass transfer phenomena in liquids have been widely examined by many researchers using MD simulations.25−28 In MD simulations, the self-diffusion coefficient is obtained using the Einstein relation or the Green−Kubo relation. Direct

1. INTRODUCTION A molecular level understanding of structure and transport properties at liquid interfaces has great importance in the study of industrial and biological processes. Experimental techniques such as nonlinear spectroscopic methods,1,2 X-ray scattering methods,3 and ellipsometry4 and computer simulations such as classical5−8 and ab initio methods9 have enhanced our molecular level understanding of liquid interfaces. However, the presence of thermal fluctuations and capillary waves at liquid−vapor interfaces makes accurate measurement of structural properties, and their analysis at the interfaces is a challenging task from the standpoint of both computational and experimental methods. The molecular dynamics (MD) simulation technique is the most commonly used computer simulation method to investigate structure and transport properties at the liquid− vapor interfaces.10−13 In the past studies,10−13 quantities were obtained on the basis of the time- and space-averaged Gibbs dividing surface, where the simulation system was divided into rectangular bins to examine the variation of the desired quantities in the interface region. In these results, most of the instantaneous and local structures at the interfaces are lost because the quantities are averaged under the influence of molecular-scale surface fluctuations at the free surfaces. To obtain clear information about the physical quantities at the interfaces, one has to consider the interface definition by accounting for the molecular-scale fluctuations. If we can calculate the physical and structural quantities that are not © XXXX American Chemical Society

Received: May 26, 2016 Revised: July 3, 2016

A

DOI: 10.1021/acs.jpcb.6b05332 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B implementation of these methods is only suitable for homogeneous fluids.29 Recently, new methods have been developed to accurately determine the diffusion coefficient in inhomogeneous regions such as interfaces.29,30 It is reported that the calculation of the self-diffusion coefficient in the parallel direction to the interface in the interface region is essentially straightforward because the inhomogeneity is only considered to exist in the direction perpendicular to the interface.30 If one wishes to find the diffusion in the perpendicular direction to the interface, then the abovementioned new methods have to be used.29,30 To obtain the local self-diffusion coefficient parallel to the interface in the interface region, one has to measure it in a thin control volume region where the main difficulty is that the molecules are likely to leave the region of interest during the calculation time. If the interface definition is not considered carefully, an attempt to define the control volume region may result in an averaging effect of two phases in the interface region due to the capillary waves. So, the interface definition has an important role in measuring the local diffusion coefficient accurately. By considering the space- and time-averaged Gibbs dividing surface, the self-diffusion coefficient of water and alkane in the liquid−vapor interface region was examined by a few researchers,11,30 where the results in the interface region were affected by the averaging effect of liquid and vapor phases. To find the diffusion coefficient more accurately in the interface region, an intrinsic interface definition proposed by our group16 was used in this work. To the best of our knowledge, there are no reports available in the literature concerning the self-diffusion of realistic fluids consisting of molecules with internal molecular motion, such as n-alkane liquids, at the liquid−vapor interfaces based on the local and instantaneous (intrinsic) interface definition. In this paper, the molecular-scale intrinsic structure and self-diffusion in the direction parallel to the intrinsic interface in the vicinity of liquid−vapor interfaces of liquids is investigated by studying liquid decane (C10H22) and tetracosane (C24H50) at various temperatures. The results were also compared with the quantities obtained with respect to the time- and spaceaveraged Gibbs dividing surface.

Figure 1. Snapshot of the decane liquid−vapor computational system with intrinsic interfaces at T = 350 K.

where rij is the distance between united atoms i and j. The values of constants kr and beq are 1.3323 × 10−18 J/Å2 and 1.54 Å, respectively. The bond bending is defined by 1 U (θ ) = kθ(θ − θ0)2 (2) 2 where θ is the bond angle. The equilibration angle θ0 is 114.0°, and the constant kθ is 8.6291 × 10−19 J/rad2. The torsional term is modeled using the Ryckaert−Bellemans (RB) potential form as follows U (φ) = V0 + V1(1 + cos φ) + V2(1 − cos 2φ) + V3(1 + cos 3φ)

where φ is the dihedral angle and the constants V0, V1, V2, and V3 are 0, 4.9018 × 10−21, − 9.4146 × 10−22, and 1.0925 × 10−20 J, respectively. The nonbonded interaction between two pseudoatoms located on different molecules and pseudoatoms located on the same molecule separated by at least four bonds is given by the Lennard-Jones (LJ) potential

2. MOLECULAR MODEL AND COMPUTATIONAL DETAILS

⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij U (rij) = 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠ LJ

To study the intrinsic structure and diffusion coefficient of liquid alkanes at the liquid−vapor interfaces, the simulation system constructed in this study is a filmwise liquid n-alkane with two liquid−vapor interfaces as depicted in Figure 1. Decane and tetracosane liquids were considered as an example of n-alkanes in this work. The united atom model was used to describe alkane molecules, where methyl (CH3−) and methylene (−CH2−) groups were treated as pseudoatoms located at carbon atoms. The NERD force field31 was used to define molecular interaction between pseudoatoms. In this force field, the pseudoatoms interact via nonbonded and bonded forces. Bonded interactions consist of bond bending, bond stretching, and torsional terms. The bond stretching is modeled in terms of a harmonic potential as follows U (rij) =

kr (rij − beq )2 2

(3)

(4)

where rij is the distance between two atoms. The energy and size parameters for CH3 and CH2 groups are εCH3 = 1.4358 × 10−21 J, εCH2 = 6.3233 × 10−22 J, and σCH3 = 3.91 Å, σCH2 = 3.93 Å, respectively. Interaction parameters between dissimilar interaction sites, like CH3 and CH2 groups, were attained using Lorentz−Berthelot mixing rules. In all simulations, a cutoff radius of 14 Å was used for the LJ interactions. The rRESPA method,32 which involves the time integration scheme for the equations of motion with multiple time scales (MTS), was used. The integration time steps used for intermolecular and intramolecular motions are 1 and 0.2 fs, respectively. Periodic boundary conditions were applied in three spatial directions. The lengths of the rectangular simulation system are Lx = Ly = 40 Å and Lz = 160 Å for decane and Lx = Ly = 60 Å and Lz = 240 Å for tetracosane. The numbers of molecules used

(1) B


Article

The Journal of Physical Chemistry B

3. DEFINITION OF THE INTRINSIC INTERFACE

in the simulation of decane and tetracosane liquids were 250 and 300, respectively. To create the liquid−vapor simulation system, initially molecules were randomly positioned at the center of the rectangular computational system in such a way to form a liquid film surrounded by vacuum. Next, the system temperature was slowly increased to the desired temperature by the velocity scaling method. Subsequently, the system was run in the NVT ensemble for 6 ns with a Nosé−Hoover chain thermostat33 for equilibration. When the system reached an equilibrium state, the thickness of the film was altered from its initial thickness as the number of molecules in the computational system was kept constant so as to have the equivalent bulk liquid density in the film at that temperature and the corresponding saturated pressure. After the equilibration, for the analysis, the data were collected for 2 ns. For the investigation of structural quantities, data were saved every 100 time steps, and for the evaluation of the self-diffusion coefficient, data were collected every 10 time steps. The simulations were performed using in-house code.

In this work, we have used the local and instantaneous interface definition developed by Kikugawa et al.16 for obtaining the structure and self-diffusion coefficient in the vicinity of liquid− vapor interfaces. The details of the interface definition are briefly summarized as follows. In this method, an instantaneous single particle density of alkane pseudoatoms is expressed as the summation of Dirac delta functions N

ρ(r) =

∑ δ(r − ri)

where N is number of pseudoatoms in the system, r is the field position, and ri is the position of the pseudoatom i. This discrete density distribution was expressed into a continuous field quantity by using the following smoothed delta function on every pseudoatom in the system

3 ⎧ ⎪(2k Δx)−3 ∏ 1 + cos π (rα − rα(par)) k Δx B(r − r(par)) = ⎨ α=1 ⎪ 0, ⎩

{

} if |r − r

where r(par) is the position of the atom. Here, rα indicates the coordinates of the grid points in the x, y ,and z directions, and Δx is the grid spacing (1 Å in this work) generated in the MD simulation box. The variable k specifies broadening of the smoothed delta function. Therefore, k defines the ability of this method to capture the molecular scale fluctuations of the interface. In this study, the value of k was chosen in such a way that the smoothing of the delta function is approximately equal to the mean intermolecular distance of the atoms. The smoothed delta function in eq 6 was used in place of Dirac delta function in eq 5 as

∑ B(r − ri) i=1

(par) α |

α

< k Δx

otherwise

(6)

4. MEASUREMENT OF IN-PLANE SELF-DIFFUSION COEFFICIENT AT THE INTRINSIC LIQUID−VAPOR INTERFACE The shape of the intrinsic interface changes with space and time. It is a challenging task to consider these spatial and temporal changes of the interface into account while counting the number of molecules in a local volume in the interface region because of the interface movement over the time. In this work, we used the slab definition changing with time based on the instantaneous intrinsic interface. During calculation of the self-diffusion coefficient based on migration of molecules in the parallel direction to the instantaneous interface, i.e., the inplane diffusion coefficient, the topological changes of the intrinsic interface over the calculation time mentioned below were not considered. This means that for counting the number of molecules in a region of interest from a particular time step we used the shape of the intrinsic interface at that time step until the desired time interval mentioned below for calculating the local diffusion coefficient in a slab. To determine the inplane self-diffusion coefficient parallel to the intrinsic interface in the vicinity of interface, we have obtained the velocity autocorrelation function (VACF) profiles because of their small relaxation time over MSD profiles. The VACF can be calculated as follows

N

ρD (r) =

(5)

i=1

(7)

The simulation box was divided into rectangular grids of equal spacing in x, y, and z directions to calculate the density value and other quantities on each grid point. To determine the position of the interface, the probability distribution of the density on every grid point evaluated by eq 7 was utilized. It showed two peaks corresponding to liquid and vapor phases. The central density value between the two density peaks gives the density value of the instantaneous interface. The instantaneous surface is defined as the isosurface which has the density value of the instantaneous interface. After obtaining the interface position, the reinitialization procedure of the level set method34,35 is adopted to determine the distance function from the intrinsic surface on each grid point. After that, the distance between all atoms and the intrinsic interface is evaluated by the interpolation technique. The distance function, ϕ, from the instantaneous interface is defined in such a way that positive values means liquid side and negative values means vapor side and zero means instantaneous interface. The gray color surfaces in Figure 1 are examples of the intrinsic interfaces obtained using the above-mentioned method in the decane liquid−vapor simulation system at 350 K.

Cv(t ) = ⟨V (0)V (t )⟩

(8)

where V∥ is a velocity vector of the center of mass of a molecule in the parallel direction to the intrinsic interface, which was obtained by projecting the velocity vector, V, onto the plane parallel to the intrinsic interface as follows V =V−

(V n )n |n|2

(9)

where n is the local normal vector to the intrinsic surface, which was evaluated from the distance function ϕ as35 C


Article


Figure 2. Density distributions based on the intrinsic interface and the Gibbs dividing surface for (a) decane at T = 300, 350, and 400 K and (b) tetracosane at T = 400, 450, and 500 K.

n=

∇ϕ |∇ϕ|

same quantities evaluated based on the time- and spaceaveraged Gibbs dividing surface are also shown in the figures, where the results are plotted by taking the Gibbs dividing surface as the origin, whose position was determined by fitting the time-averaged number density profile along the z-axis to a hyperbolic tangent function.37 Even though the shapes of the instantaneous interfaces on the top and bottom sides of the liquid film are different in Figure 1, the characteristics on both sides are expected to be similar like the averaged density profiles on both sides, which were found to be equivalent. Thus, the results of the structural quantities and self-diffusion coefficient obtained for both interfaces were averaged for better statistical accuracy. Figure 2 shows the intrinsic density profiles for decane at T = 300, 350, 400 K and for tetracosane at T = 400, 450, 500 K. For the case of decane, the obtained values of bulk liquid densities at T = 300, 350, 400 K are 705, 665, 621 kg/m3. For tetracosane, the obtained values of bulk liquid densities at T = 400, 450, 500 K are 721, 683, 643 kg/m3. For reference, the critical point temperature of NERD decane and NERD tetracosane are assumed to be 617 and 805 K, respectively.31 When defining the intrinsic interface, the values of k in eq 6 were selected for decane at T = 300, 350, and 400 K to be 7, 8, and 9 and for tetracosane at T = 400, 450, 500 K to be 8, 9, 10, respectively. Figure 2 illustrates that the intrinsic density profiles are characterized by oscillations at the free surface,

(10)

To obtain the variation of the self-diffusion coefficient in the interface region, separate measuring regions having a width of 3 Å were assigned with respect to the distance from the intrinsic surface. During calculation of the ensemble average of the VACF in a selected local region of interest, only the molecules continuously present in the region of interest over the entire time interval τ were considered. By applying the Green−Kubo relation, the self-diffusion coefficient in the parallel direction is calculated from the VACF as follows:36 Di =

1 lim 2 τ →∞

∫0

τ

⟨V , i(0) ·V , i(t )⟩ dt

(11)

5. RESULTS AND DISCUSSION 5.1. Intrinsic Density Profiles. Before studying the variation of the self-diffusion coefficient in the vicinity of the liquid−vapor interfaces, we first discuss the molecular structure. Distributions of the self-diffusion coefficient and structural quantities such as density, orientation of molecules, and their shape, which in this case were of the liquids in the vicinity of the liquid−vapor interface obtained by tracking the molecularscale fluctuations, are presented with respect to the distance from the intrinsic interface. For the sake of comparison, the D


Article


Figure 3. Orientation order parameter distributions based on the intrinsic interface and the Gibbs dividing surface for (a) decane at T = 300, 350, and 400 K, and (b) tetracosane at T = 400, 450, and 500 K. Density profiles are shown with dashed lines.

distance function from the instantaneous interface. ⟨ ⟩ indicates that the average is taken over all vectors located in a specified region and time steps. The range of this parameter is [−0.5, 1]. A negative value indicates preferential alignment parallel to the intrinsic interface, whereas a positive value suggests preferential alignment perpendicular to the intrinsic interface. The orientation order parameter value of −0.5 and 1 suggests that molecules are completely parallel and perpendicular to the interface, respectively. Zero indicates completely random arrangement of molecular segments. To obtain the variation of the orientation order parameter as a function of distance to the intrinsic interface, the measuring slabs, away from each other, having a width of 1 Å were assigned with respect to the distance from the intrinsic interface. A vector linking between two atoms was allotted to the region where the central point between the two atoms was positioned. Figure 3 illustrates the distributions of the orientation order parameter obtained based on the intrinsic interface and the Gibbs dividing surface for decane and tetracosane at three temperatures. Density distributions are also given in Figure 3. The interface region can be divided into three regions with respect to the orientation order parameter obtained based on the instantaneous interface, which is similar to the discussion in our previous paper11 based on the time- and space-averaged Gibbs dividing interface. In the innermost layer, which is close to the bulk liquid side, the orientation order parameter

which clearly indicates the partial layering of molecules in the vicinity of the interface. The oscillations extend into the bulk liquid and disappear ∼10 Å away from the interface. It was reported that the oscillations in the intrinsic density profile disappear at ∼11 Å from the liquid−vapor interface for dodecane at 400 K,15 which is consistent with our results. These oscillations are not found in the density profiles obtained based on the Gibbs dividing surface due to the time- and spaceaveraging caused by the thermal fluctuations of the interface. The oscillatory nature is more pronounced at lower temperatures and almost disappears at higher temperatures. This suggests that the degree of layering of molecules at the liquid− vapor interfaces decreases with an increase in temperature. A similar kind of oscillatory nature is observed for decane and tetracosane in the examined temperature range. 5.2. Molecular Alignment and Chain Conformation at the Interfaces. In order to determine the average molecular alignment of molecules at the free surfaces, the following orientation order parameter was investigated. It is defined as follows H(ϕ) =

1 ⟨3 cos2 λ − 1⟩ 2

(12)

where λ is defined as the angle between a vector joining a pair of united atoms, which are two atoms away in a molecule, and the local normal vector to the instantaneous interface. ϕ is a E


Article


Figure 4. For tetracosane at T = 400, 450, and 500 K: (a) total radius of gyration distributions based on the intrinsic interface and the Gibbs dividing surface and (b) components of radius of gyration distributions based on the intrinsic interface and the Gibbs dividng surface. Density profiles are shown with dashed lines.

obtained based on the intrinsic interface has a positive value and shows oscillations in the vicinity of the interface according to the oscillations in the intrinsic density profile. This indicates that some segments of molecules are bent in the perpendicular direction to the interface in between the layers and are preferentially parallel to the interface in the subsequent layer. In the middle layer, adjacent to the innermost layer, the orientation order parameter obtained based on the intrinsic interface exhibits a negative value, which suggests that molecules are preferentially parallel to the intrinsic interface. In the third layer, close to the vapor side, the orientation parameter approached a positive value and is expected to become zero in the vapor phase since the molecules are randomly oriented in the vapor phase. Because of less sampling in the vapor phase, the orientation order parameters have not yet approached to zero in this study. In the innermost layer, oscillatory characteristics seen in the orientation order parameter obtained based on the intrinsic interface are not clearly visible in the profiles obtained based on the Gibbs dividing surface, due to the surface fluctuations. In the middle layer, the values of orientation order parameter obtained based on the intrinsic interface exhibit a significantly smaller value than those obtained based on the time- and space-averaged Gibbs dividing surface. This means that the preference of molecules to be parallel to the interface is observed more

clearly when studied from a viewpoint of the intrinsic interface. In all three layers, preferences of molecular orientation in each layer are less obvious at higher temperatures. It is also observed that the ordering parameter close to the intrinsic interface for tetracosane at T = 400 K is approximately equal to −0.5, which indicates that molecules are almost parallel to the intrinsic interface, whereas the decane molecules are significantly less parallel to the intrinsic surface at the same temperature, presumably because the studied temperature is relatively close to the critical point of decane where the thermal agitation is relatively higher than that for tetracosane. To acquire information about molecular shape and its variation at the interfaces, the distributions of the radius of gyration (ROG) of molecules were investigated. The ROG gives the squared average distance between the center of mass of a molecule to all the other atoms in the molecule. The meansquared total ROG and mean-squared ROG in the z- and xdirections as a function of distance from the instantaneous interface are shown in Figure 4 for tetracosane at three temperatures. Results obtained based on the Gibbs dividing surface are also presented in Figure 4. As expected, the meansquared total ROG is constant in the bulk liquid region as the molecular shape is uniform. It can be observed that the meansquared total ROG obtained based on the intrinsic interface decreases slightly in the partial layered region and then F


Article


Figure 5. Gauche fraction distributions based on the intrinsic interface and the Gibbs dividing surface for (a) decane at T = 300, 350, and 400 K and (b) tetracosane at T = 400, 450 and 500 K. Density profiles are shown with dashed lines.

For the case of decane at 300 K, the accuracy of the gauche fraction in the vapor side of the interface region is poor because of the small sampling number. The gauche fraction distributions together with total ROG distributions and orientation order parameter distributions suggest that molecules have a slightly globular shape in the partial layered region and slightly more linear shape in the vapor side of the interface region. These kinds of tendencies were not observed when the quantities were measured with respect to the Gibbs diving surface. The present gauche distributions obtained based on the Gibbs dividing surface are consistent with the results reported for tetracosane in the literature.38 The molecular conformations at interface regions were primarily influenced by the molecular interaction between the molecules at the interfaces with the dense liquid phase. 5.3. Diffusion Coefficient Variation at the Liquid− Vapor Interfaces. The characteristic VACF versus time curves obtained based on the intrinsic interface and the Gibbs dividing surface for tetracosane at 400 K are shown in Figure 6. An example of the positions of slabs considered in the study is shown on the density profiles in Figure 6a. It should be noted that the shape of the slab is not rectangular when the quantities are obtained based on the intrinsic interface. The VACF distributions obtained in the bulk liquid region and in the vicinity of interfaces are illustrated in Figure 6b. It is observed that the VACF distributions obtained based on the intrinsic

marginally increases in the adjacent region to the layered region toward the vapor phase. This implies that the molecule shape is marginally altered two times in the interface region, which was not observed in the distributions obtained with respect to the Gibbs dividing surface. In the interface region the difference in the components of the mean-squared ROG in the z- and xdirections is much more obvious in the distributions based on the intrinsic interface than those that are obtained based on the Gibbs dividing surface. This also suggests that the chain molecules are more noticeably flattened in the z-direction than in the other two directions. In addition to the ROG distributions, gauche fraction profiles were examined to provide more details about the molecular shape change in the vicinity of interfaces. The gauche fraction distributions obtained based on the intrinsic interface and the Gibbs dividing surface for decane and tetracosane at three temperatures are shown in Figure 5. In order to obtain the gauche fraction in a region of interest, the total number of gauche angles was divided by the total number of dihedral angles in that region over the sampling time. A gauche state is defined as having a dihedral angle between −120 and 120°. As expected, the gauche fraction is constant in the bulk liquid region for both decane and tetracosane at all temperatures. Gauche fraction distributions obtained based on the intrinsic interface increase marginally in the partial layered region and then slightly decrease in the vapor side of the interface region. G


Article


Figure 6. For tetracosane at T = 400 K: (a) division of slab regions on density profiles based on the intrinsic interface and the Gibbs dividing surface; (b) VACF distributions based on the intrinsic interface and the Gibbs dividing surface.

interface have similar characteristics to that of VACF distributions obtained based on the Gibbs dividing surface in the bulk liquid region. The VACF distributions obtained based on the intrinsic interface decay noticeably more slowly compared to the VACF distributions calculated with respect to the Gibbs dividing surface in slab no. 1, which suggests that the self-diffusion coefficient obtained with respect to the intrinsic interface is different from that of the diffusion coefficient measured based on the Gibbs dividing surface in slab no. 1. Figure 7 shows the variation of the self-diffusion coefficient according to molecular migration parallel to the instantaneous interface in the vicinity of the liquid−vapor interface as a function of distance from the intrinsic interface for decane and tetracosane at three temperatures. The distributions of the self-diffusion coefficient parallel to the Gibbs dividing surface are also shown in Figure 7. The selfdiffusion coefficients in the bulk liquid region from both approaches show good agreement. It is also noted that the selfdiffusion coefficients obtained from both approaches increases in the interface region for both decane and tetracosane. It can also be seen that the self-diffusion coefficient increases monotonically as the density decreases monotonically at the interface when the quantities were obtained with respect to the Gibbs dividing surface. Contrary to this, the self-diffusion coefficient increases monotonically as the density changes in a fluctuating manner at the interfaces when the quantities were

obtained with respect to the intrinsic interface. From these results and the results reported in the literature38,39 it can be understood that the correlation between the density and the self-diffusion coefficient is not sufficient to understand the diffusion characteristics in the vicinity of interfaces, as well as in the bulk liquid region. For example, the self-diffusion coefficient of hexane in the binary liquid mixture of hexane-tetracosane increases with an increase in molar fraction of hexane in the bulk liquid region, which means both the density and selfdiffusion coefficient of hexane increase simultaneously in the mixture.38 It has been suggested that the increase in the selfdiffusion coefficient of hexane in the hexane-tetracosane mixture is due to the larger free volume available in the mixture with the increase in molar fraction of hexane. From these observations it can be interpreted that even though the density increases in a local region, the self-diffusion coefficient might increase when the local free volume around the molecules is higher. At the same time, it has been observed that the nature of the self-diffusion characteristics of both hexane and tetracosane in their binary mixture in the liquid− vapor interface region are mostly governed by the changes in local free volume, which is influenced by the local structure of liquid at the interfaces.38 From the results in Figures 2−5, it is known that the local structure of liquid alkane in the vicinity of liquid−vapor interface region is significantly different from that in the bulk liquid region, a conclusion which is obtained from H


Article


Figure 7. Self-diffusion coefficient distributions based on the intrinsic interface and the Gibbs dividing surface for (a) decane at T = 300, 350, and 400 K and (b) tetracosane at T = 400, 450, and 500 K. Density profiles are shown with dashed lines.

relation, the self-diffusion coefficient parallel to the intrinsic interface was obtained in the vicinity of liquid−vapor interfaces. It was found that this self-diffusion coefficient varies steeply in the liquid−vapor interface region in a manner similar to that for the intrinsic density. Therefore, we can conclude that the use of a local and instantaneous definition of interface can capture more local information regarding not only structural properties but also molecular dynamics and transport properties that are intrinsically present at the interface.

both types of interface definitions. Therefore, we believe that the increase in the self-diffusion coefficient at the liquid−vapor interfaces is due to the increase in the local free volume around the molecules in the vicinity of interfaces. The main finding of this work is that the variation of the self-diffusion coefficient in the interface region is captured by the present intrinsic interface analysis much more obviously, in a steep manner, than that found using the conventional time- and space-averaged Gibbs dividing surface. This is because the local and instantaneous definition of interface enables us to capture the interfacial structures and quantities which are intrinsically present at the interface. The increase in self-diffusion coefficient is lower in the vicinity of interfaces when the quantities are calculated with respect to the Gibbs dividing surface due to the averaging over the fluctuations of the instantaneous interface.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail:[email protected]. Tel: +81-22-2175872. Fax: +81-22-217-5872. Notes

The authors declare no competing financial interest.

6. CONCLUSIONS Molecular structures and self-diffusion coefficients at the intrinsic liquid−vapor interfaces of decane and tetracosane at three temperatures were examined using MD simulations. The partial layering structure was clearly observed at the intrinsic liquid−vapor interfaces of alkanes. In the interface region, alkane molecules are preferentially oriented more parallel to the intrinsic interface as compared with those observed at the space-fixed Gibbs dividing surface. This tendency is less obvious at higher temperatures. By using the Green−Kubo

■

ACKNOWLEDGMENTS The work reported in this paper was partially supported by a Grant-in-Aid for Scientific Research and the Global COE Program “International COE of Flow Dynamics” by the Japan Society for the Promotion of Science (JSPS). Numerical simulations were performed on the SGI Altix UV1000 and UV2000 at the Advanced Fluid Information Research Center, Institute of Fluid Science, Tohoku University. I


Article


■

(24) Bresme, F.; Chacon, E.; Tarazona, P. Force Field Dependence on the Interfacial Structure of Oil-Water Interfaces. Mol. Phys. 2010, 108, 1887−1898. (25) Dysthe, D. K.; Fuchs, A. H.; Rousseau, B. Fluid Transport Properties by Equilibrium Molecular Dynamics. III. Evaluation of United Atom Interaction Potential Models for Pure Alkanes. J. Chem. Phys. 2000, 112, 7581−7590. (26) Padilla, P.; Toxvaerd, S. Self-diffusion in n-Alkane Fluid Models. J. Chem. Phys. 1991, 94, 5650−5654. (27) Mondello, M.; Grest, G. S. Molecular Dynamics of Linear and Branched Alkanes. J. Chem. Phys. 1995, 103, 7156−7165. (28) Nieto-Draghi, C. N.; Ungerer, P.; Rousseau, B. Optimization of the Anisotropic United Atoms Intermolecular Potential for n-Alkanes: Improvement of Transport Properties. J. Chem. Phys. 2006, 125, 044517. (29) Liu, P.; Harder, E.; Berne, B. J. On the Calculation of Diffusion Coefficients in Confined Fluids and Interfaces with an Application to the Liquid-Vapor Interface of Water. J. Phys. Chem. B 2004, 108, 6595−6602. (30) Wick, C. D.; Dang, L. X. Diffusion at the Liquid-Vapor Interface of an Aqueous Ionic Solution Utilizing a Dual Simulation Technique. J. Phys. Chem. B 2005, 109, 15574−15579. (31) Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. On the Simulation of Vapor-Liquid Equilibria for Alkanes. J. Chem. Phys. 1998, 108, 9905− 9911. (32) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible Multiple Time Scale Molecular Dynamics. J. Chem. Phys. 1992, 97, 1990−2001. (33) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (34) Sussman, M.; Smereka, P.; Osher, S. A Level Set Approach for Computing Solutions to Incompressible Two-phase Flow. J. Comput. Phys. 1994, 114, 146−159. (35) Osher, S.; Fedkiw, R. Level Set Methods and Dynamic Implicit Surfaces; Springer: New York, 2003. (36) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press, Inc.: New York, 1987. (37) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Dover: New York, 1982. (38) Chilukoti, H. K.; Kikugawa, G.; Ohara, T. Self-diffusion Coefficient and Structure of binary n-Alkane Mixtures at the liquidVapor Interfaces. J. Phys. Chem. B 2015, 119, 13177−13184. (39) Harmandaris, V. A.; Angelopoulou, D.; Mavrantzas, V. G.; Theodorou, D. N. Atomistic Molecular Dynamics Simulation of Diffusion in Binary Liquid n-Alkane Mixtures. J. Chem. Phys. 2002, 116, 7656−7665.

REFERENCES

(1) Lu, R.; Gan, W.; Wu, B.-H.; Chen, H.; Wang, H.-F. Vibrational Polarization Spectroscopy of CH Stretching Modes of the Methylene Groups at the Vapor/Liquid Interfaces with Sum Frequency Generation. J. Phys. Chem. B 2004, 108, 7297−7306. (2) Miranda, P. B.; Shen, Y. R. Liquid Interfaces: A Study by SumFrequency Vibrational Spectroscopy. J. Phys. Chem. B 1999, 103, 3292−3307. (3) Fradin, C.; Braslau, A.; Luzet, D.; Smilgies, D.; Alba, M.; Boudet, N.; Mecke, K.; Daillant, J. Reduction in the Surface Energy of Liquid Interfaces at Short Length Scales. Nature 2000, 403, 871−874. (4) Cho, J-H. J.; Law, B. M. Surface Orientational Order at LiquidVapor Interfaces Induced by Dipole-Image-Dipole Interactions. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 67, 031605. (5) Kikugawa, G.; Takagi, S.; Matsumoto, Y.; Ohara, T. A Molecular Dynamics Study on the Local Structure of Liquid-Vapor Interface of Water and L-J Fluid. J. Therm. Sci. Technol. 2008, 3, 234−240. (6) Chacón, E.; Tarazona, P.; Alejandre, J. The Intrinsic Structure of the Water Surface. J. Chem. Phys. 2006, 125, 014709. (7) Chacón, E.; Tarazona, P. Intrinsic Profiles beyond the Capillary Wave Theory: a Monte Carlo Study. Phys. Rev. Lett. 2003, 91, 166103. (8) Chowdhary, J.; Ladanyi, B. M. Water-Hydrocarbon Interfaces: Effect of Hydrocarbon Branching on Interfacial Structure. J. Phys. Chem. B 2006, 110, 15442−15453. (9) Hernández-Cobes, J. H.; Saint-Martin, H. S.; Mackie, A. D.; Vega, L. F.; Ortega-Blake, I. O. Water Liquid-Vapor Equilibria Predicted by Refined ab initio Derived Potentials. J. Chem. Phys. 2005, 123, 044506. (10) Harris, J. G. Liquid-Vapor Interfaces of Alkane Oligomers. Structure and Tthermodynamics from Molecular Dynamics Simulations of Chemically Realistic Models. J. Phys. Chem. 1992, 96, 5077− 5086. (11) Chilukoti, H. K.; Kikugawa, G.; Ohara, T. A Molecular Dynamics Study on Transport Properties and Structure at the LiquidVapor Interfaces of Alkanes. Int. J. Heat Mass Transfer 2013, 59, 144− 154. (12) Tsige, M.; Grest, G. S. Surface Tension and Surface Orientation of Perfluorinated Alkanes. J. Phys. Chem. C 2008, 112, 5029−5035. (13) Cao, B.-Y.; Xie, J.-F.; Sazhin, S. S. Molecular Dynamics Study on Evaporation and Condensation of n-Dodecane at Liquid-Vapor Phase Equilibria. J. Chem. Phys. 2011, 134, 164309. (14) Esenturk, O.; Walker, R. A. Surface Vibrational Structure at Alkane Liquid/Vapor Interfaces. J. Chem. Phys. 2006, 125, 174701. (15) Bresme, F.; Chacón, E.; Tarazona, P.; Tay, K. Intrinsic Structure of Hydrophobic Surface: The Oil-Water Interface. Phys. Rev. Lett. 2008, 101, 056102. (16) Kikugawa, G.; Takagi, S.; Matsumoto, Y. A Molecular Dynamics Study on Liquid-Vapor Interface Adsorbed by Impurities. Comput. Fluids 2007, 36, 69−76. (17) Partay, L. B.; Hantal, G.; Jedlovszky, P.; Vincze, A.; Horvai, G. A New Method for Determining the Interfacial Molecules and Characterizing the Surface Roughness in Computer Simulations. Application to the Liquid-Vapor Interface of Water. J. Comput. Chem. 2008, 29, 945−956. (18) Willard, A. P.; Chandler, D. Instantaneous Liquid Interfaces. J. Phys. Chem. B 2010, 114, 1954−1958. (19) Chacon, E.; Tarazona, P. Intrinsic Profiles beyond the Capillary Wave Theory: A Monte Carlo Study. Phys. Rev. Lett. 2003, 91, 166103. (20) Geysermans, P.; Pontikis, V. Interfacial Layering and Capillary Roughness in Immiscible Liquids. J. Chem. Phys. 2010, 133, 074706. (21) Duque, D.; Tarazona, P.; Chacón, E. Diffusion at the LiquidVapor Interface. J. Chem. Phys. 2008, 128, 134704. (22) Bresme, F.; Chacón, E.; Tarazona, P. Molecular Dynamics Investigation of the Intrinsic Structure of Water-Fluid Interfaces via the Intrinsic Sampling Method. Phys. Chem. Chem. Phys. 2008, 10, 4704−4715. (23) Jorge, M.; Natalia, M.; Cordeiro, D. S. Intrinsic Structure and Dynamics of the Water /Nitrobenzene Interafec. J. Phys. Chem. C 2007, 111, 17612−17626. J


Structure and Mass Transport Characteristics at the Intrinsic Liquid

Recommend Documents