Dielectric properties of liquid ethanol: a computer simulation study

Dielectric properties of liquid ethanol. A computer simulation study Leonor Saiz, Elvira Guàrdia, and Joan-Àngel Padró Citation: The Journal of Chemical Physics 113, 2814 (2000); doi: 10.1063/1.1305883 View online: http://dx.doi.org/10.1063/1.1305883 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/113/7?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Liquid 1-propanol studied by neutron scattering, near-infrared, and dielectric spectroscopy J. Chem. Phys. 140, 124501 (2014); 10.1063/1.4868556 Are dipolar liquids ferroelectric? Simulation studies J. Chem. Phys. 126, 104506 (2007); 10.1063/1.2672734 Free energy of liquid water from a computer simulation via cell theory J. Chem. Phys. 126, 064504 (2007); 10.1063/1.2434964 Dielectric relaxation of supercritical water: Computer simulations J. Chem. Phys. 113, 3499 (2000); 10.1063/1.1289919 Spectroscopic and dielectric properties of liquid water: A molecular dynamics simulation study J. Chem. Phys. 109, 4911 (1998); 10.1063/1.477102 This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.83.95.18 On: Wed, 14 Jan 2015 14:53:57 Dielectric properties of liquid ethanol. A computer simulation study Leonor Saiz Departament de Fı ´sica Fonamental, Universitat de Barcelona, Diagonal 647, E-08028 Barcelona, Spain Elvira Gua `rdia Departament de Fı ´sica i Enginyeria Nuclear, Universitat Polite `cnica de Catalunya, Sor Eula `lia D’Anzizu, Campus Nord, Mo `dul B4-B5, E-08034 Barcelona, Spain Joan-A `ngel Padro ´ Departament de Fı ´sica Fonamental, Universitat de Barcelona, Diagonal 647, E-08028 Barcelona, Spain 共Received 17 April 2000; accepted 18 May 2000兲 Static and dynamic dielectric properties of liquid ethanol have been studied as a function of the wave-vector number by computer simulation. Molecular dynamics simulations at room temperature have been performed using the optimized potentials for liquid simulations 共OPLS兲potential model proposed by Jorgensen 关J. Phys. Chem. 90, 1276 共1986兲兴. The time dependent correlation functions of the longitudinal and transverse components of the dipole density as well as the individual and total dipole moment autocorrelation functions have been calculated. The infrared spectra and the dielectric relaxation of the liquid have been also analyzed. Results have been compared with the available experimental data. Special attention has been dedicated to investigate the molecular origin of the different analyzed properties. © 2000 American Institute of Physics. 关S0021-9606共00兲51531-1兴 I. INTRODUCTION The study of dielectric properties is fundamental for a deep understanding of polar liquids. Computer simulations provide a suitable tool to analyze these properties at molecular level, providing a direct route from microscopic details to macroscopic properties of experimental interest. The molecular dynamics 共MD兲simulation method is appropriate for this purpose since it enables one to obtain both static and dynamic properties. However, difficulties arise due to the slow relaxation of the time correlation functions of interest and to the sensitivity to the long-range intermolecular interactions of the quantities involved. To overcome the first difficulty, extremely long simulations should be performed whereas fairly large systems are required to overcome the latter. Moreover, statistically accurate calculations of these collective properties will also require very long simulations since only one value of a given collective quantity is obtained at each time step. In the case of polar protic systems, the study of dielectric properties is especially relevant due to the important role played by these fluids as solvents in chemistry 共for a review on computer simulation of hydrogen bonded liquids, see Ref. 1兲. Among these liquids, water has focused the interest of most of the computer simulation studies. In the case of alcohols, only the dielectric properties of liquid methanol2–4 and methanol–water mixtures5have been analyzed by computer simulation. Many experimental studies of both static and dynamic dielectric properties of hydrogen bonded liquids have been carried out, but in the case of ethanol, the latter are restricted to the microwave region of the spectra6and to the infrared region7up to about 34 cm⫺1. To achieve a better understanding of the dielectric properties of monohydric alcohols, we show in this article the results obtained for liquid ethanol by means of MD simulations. We have computed the static dielectric constant and the transverse and longitudinal components of the dielectric permittivity at different wave-vector numbers paying special attention to the low wave-vector limit, which provides an alternative route to calculate the dielectric permittivity. We have also determined the time correlation functions of the longitudinal and transverse components of the dipole density as a function of the wave-vector number, as well as the total dipole moment autocorrelation function and its self- and distinct components. In addition, the dielectric relaxation and the far infrared spectra have been investigated. It is worth emphasizing that previous computer simulations of liquid ethanol were restricted to the computation of individual properties at different thermodynamic states.8–13 The article is organized as follows. In Sec. II, a theoretical background is provided for the necessary definitions. Details of the simulations are given in Sec. III. Static dielectric properties are reported in Sec. IV. Section V is devoted to the results concerning dynamical dielectric properties, including time correlation functions and frequency spectra. The dynamics underlying the principal dielectric band in the microwave frequency region is analyzed. The infrared spectra is investigated in Sec. VI. Comparison between experimental and simulation data is presented throughout the article. The main results and conclusions are summarized in Sec. VII. II. THEORETICAL FRAMEWORK For an infinite system with cubic or spherical symmetry, the dielectric permittivity and susceptibility tensors are related through JOURNAL OF CHEMICAL PHYSICS VOLUME 113, NUMBER 7 15 AUGUST 2000 28140021-9606/2000/113(7)/2814/9/$17.00 © 2000 American Institute of Physics This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.83.95.18 On: Wed, 14 Jan 2015 14:53:57 ⑀ L共k, ␻ 兲⫺ ⑀ ⬁ ⑀ L共k, ␻ 兲 ⑀ ⬁ ⫽ ␹ L 0共k, ␻ 兲共1兲 and ⑀ T共k, ␻ 兲⫺ ⑀ ⬁⫽ ␹ T 0共k, ␻ 兲共2兲 for the longitudinal and transverse components, respectively. ⑀ ⬁is the dielectric constant at optical frequencies. 共For a model system with nonpolarizable molecules, ⑀ ⬁⫽1.兲Linear response theory relates the susceptibility tensor ␹ 0(k ជ , ␻ )to correlation functions of the Fourier components of the dipole density of the system M ជ (k ជ ,t). The two components of M ជ (k ជ ,t), M ជ 共k ជ ,t兲⫽M ជ L共k ជ ,t兲⫹M ជ T共k ជ ,t兲,共3兲 can be expressed as follows: M ជ L共k ជ ,t兲⫽兺 j⫽1 N k ˆk ˆ• ␮ ជ j共t兲exp共ik ជ r ជ j CM共t兲兲 共4兲 and M ជ T共k ជ ,t兲⫽兺 j⫽1 N 共1⫺k ˆk ˆ兲• ␮ ជ j共t兲exp共ik ជ r ជ j CM共t兲兲,共5兲 where ␮ ជ j(t) is the dipole moment of the jth molecule and r ជ j CM(t) the position of its center-of-mass. The components of the susceptibility tensor are given by14,15 ␹ A 0共k, ␻ 兲⫽ 具 兩 M ជ A共k ជ ,0兲 兩 2 典 ␯ AVkBT ⑀ 0关1⫹i ␻ ⌽A共k, ␻ 兲兴,共6兲 where A⫽Lor T共with ␯ L⫽1 and ␯ T⫽2兲,Vis the volume of the sample, ⑀ 0is the vacuum permittivity, and ⌽A(k, ␻ )is the Fourier transform, ⌽A共k, ␻ 兲⫽ 冕 0 ⬁dt⌽A共k,t兲exp共i ␻ t兲,共7兲 of the normalized dipole density correlation function ⌽A(k,t) defined by ⌽A共k,t兲⫽ 具 M ជ A共k ជ ,t兲•M ជ A共⫺k ជ ,0兲 典 具 兩 M ជ A共k ជ ,0兲 兩 2 典 .共8兲 Because of the short-ranged nature of the correlations which determine the dielectric permittivity, ⑀ (k, ␻ ) must become independent of the wave-vector number when kis sufficiently small, that is lim k→0 ⑀ T共k, ␻ 兲⫽lim k→0 ⑀ L共k, ␻ 兲⫽ ⑀ 共 ␻ 兲.共9兲 The expressions for the components of the static wave-vector dependent dielectric constant can be deduced from Eqs. 共1兲 and 共2兲. For the static case, i.e., ␻ ⫽0, they reduce to ⑀ L共k兲⫺ ⑀ ⬁ ⑀ L共k兲 ⑀ ⬁ ⫽y 具 兩 M ជ L共k ជ ,0兲 兩 2 典 N ␮ 2共10兲 and ⑀ T共k兲⫺ ⑀ ⬁⫽y 具 兩 M ជ T共k ជ ,0兲 兩 2 典 2N ␮ 2,共11兲 where y⫽ ␳ ␮ 2 kBT ⑀ 0,共12兲 with ␳ the molecular number density and ␮ the molecular dipole moment. Another quantity of interest,16 ⑀ R(k), results by dividing both sides of Eqs. 共10兲and 共11兲 ⑀ R共k兲⫽ 具 兩 M ជ T共k ជ ,0兲 兩 2 典 2 具 兩 M ជ L共k ជ ,0兲 兩 2 典 ,共13兲 or equivalently ⑀ R共k兲⫽xT共k兲 xL共k兲,共14兲 where xT共k兲⫽ 具 兩 M ជ T共k ជ ,0兲 兩 2 典 2N ␮ 2共15兲 and xL共k兲⫽ 具 兩 M ជ L共k ជ ,0兲 兩 2 典 N ␮ 2.共16兲 In the low klimit the following condition must be fulfilled: lim k→0 ⑀ L共k兲⫽lim k→0 ⑀ T共k兲⫽ ⑀ .共17兲 Alternatively, lim k→0 ⑀ R共k兲⫽ ⑀ .共18兲 Equations 共9兲and 共17兲provide a route to determine ⑀ 共 ␻ 兲and ⑀ , respectively. However, due to the periodic boundary conditions used during the computer simulations, the kvectors which may be studied are restricted to k ជ ⫽2 ␲ L共l,m,n兲,共19兲 where l, m, n are integers and Lis the length of the cubic box. Thus, the smallest wave-vector which may be studied is kmin⫽2 ␲ /L. Whether the low klimit is reached can be tested through the consistency of the calculated ⑀ L(kmin , ␻ ) and ⑀ T(kmin , ␻ ). Another approach for the determination of ⑀ 共 ␻ 兲 is based on the calculation of the total dipole moment of the primitive cube of the simulation M ជ (t)⫽兺j⫽1 N ␮ ជ j(t). If the Ewald summation technique is used to handle with the longrange interactions and conducting walls boundary conditions are assumed, the relation between ⑀ 共 ␻ 兲and the total dipole moment correlation is given by17,18 ⑀ 共 ␻ 兲⫺ ⑀ ⬁⫽ 具 兩 M ជ 共0兲 兩 2 典 3VkBT ⑀ 0关1⫹i ␻ ⌽共 ␻ 兲兴,共20兲 where ⌽共 ␻ 兲is the Fourier transform of the M ជ autocorrelation function 2815J. Chem. Phys., Vol. 113, No. 7, 15 August 2000 Dielectric properties of liquid ethanol This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.83.95.18 On: Wed, 14 Jan 2015 14:53:57 ⌽共t兲⫽ 具 M ជ 共t兲•M ជ 共0兲 典 具 兩 M ជ 共0兲 兩 2 典 .共21兲 Finally, the static dielectric constant is given by ⑀ ⫽ ⑀ ⬁⫹y 3Gk,共22兲 where Gkis the finite system Kirkwood correlation factor Gk⫽ 具 兩 M ជ 共0兲 兩 2 典 N ␮ 2.共23兲 Gkaccounts for the correlations of the total dipole moment. It can be expressed as Gk⫽ 具 兩 M ជ 共0兲 兩 2 典 N ␮ 2⫽1⫹N⫺1 ␮ 2 具 ␮ ជ i• ␮ ជ j 典 ,共24兲 where the averages of the product ␮ ជ i• ␮ ជ jare performed for different molecules (i⫽j). Gkreflects the degree of correlation between the orientation of neighboring molecules. So, Gk⫽1 if there is no correlation, Gk⬎1 when molecular dipole moments tend to orient themselves parallel to each other, and Gk⬍1 when molecular dipole moments tend to orient themselves in opposite directions. Once ⑀ has been obtained from Eq. 共22兲the Kirkwood factor gkwhich describes orientational correlations in an infinite sample may be obtained by means of the relationship gk⫽共2 ⑀ ⫹ ⑀ ⬁兲 3 ⑀ Gk.共25兲 III. COMPUTER SIMULATION DETAILS Molecular dynamics simulations of liquid ethanol were performed at room temperature (T⫽298 K) and the experimental density at normal pressure (d⫽0.7873 g/cm3). The number of molecules in the central box was N⫽125, which corresponds to a cubic box of side length L⫽22.99Å. Periodic boundary conditions were used. We adopted the optimized potentials for liquid simulations 共OPLS兲potential and the molecular model proposed by Jorgensen in Ref. 8. Each ethanol molecule consists of four interaction sites corresponding to the two methyl groups and the oxygen and the hydrogen atoms of the hydroxyl group. Hydrogen atoms of the methyl groups are not explicitly considered. Bond lengths and bond angles are kept constant, and the only intramolecular motions considered in the simulations are the torsions around the central C–O bond. Interactions between sites are described by Lennard-Jones and Coulomb terms. The partial charges and molecular geometry yield a molecular dipole moment ␮ ⫽2.22D (1 D⫽3.335⫻10⫺30 ccm) that is significantly larger than the experimental gas-phase value19 ␮ gas⫽1.69D. In recent work,11,13 it has been shown that the OPLS potential reproduces fairly well the thermodynamic properties, the structure, and single particle dynamics of liquid ethanol at different thermodynamic states. The equations of motion were integrated using the algorithm proposed by Berendsen and co-workers in Ref. 20, which consists of a leap-frog Verlet algorithm with a weak coupling to a thermal bath. We used a time step of 2.5 fs. Constraints were handled by means of the SHAKE method.21 We performed two MD runs consisting of an initial equilibration period of about 100 ps and an equilibrium period of 1300 ps during which the properties were calculated. The results reported in the next sections are the average of those obtained from these two independent runs. The short-range forces were truncated at half the box length and the Ewald summation with conducting boundary conditions22 was used for Coulomb interactions. IV. STATIC DIELECTRIC PROPERTIES Let us first examine the static properties obtained from the MD simulations following the scheme described previously. To this purpose we consider the more relevant quantities, namely, the Kirkwood correlation factor, the static dielectric constant, and the longitudinal and transverse wavevector dependent static dielectric constants. A. k Ä0 properties As discussed in Sec. II, one of the approaches to the computation of the static dielectric constant by MD simulations is from Eqs. 共22兲,共23兲, and 共25兲共in our case y ⫽15.49兲. The values of Gk, ⑀ , and gkfound in this way are given in Table I and they are compared with the available experimental data. The Kirkwood gfactor Gkobtained from MD was greater than 1, as expected since there is a parallel alignment of the dipole moments of near-neighboring molecules.11 The convergence of the cumulative average of the static dielectric constant as a function of the simulation FIG. 1. Convergence of the static dielectric constant as a function of the simulation length for each run. The dashed line corresponds to the final value of ⑀ obtained by averaging the two runs. TABLE I. Static dielectric constant and Kirkwood gfactor. Method Gk ⑀⑀ ⬁ ⑀ ⫺ ⑀ ⬁gk This worka2.9⫾0.4 16⫾21 15⫾22.0⫾0.2 Expt.b24.32 2.69 21.63 Expt.c24.35 1.93 22.42 aStatistal errors have been computed following the blocking method of Ref. 23. bReference 6. cReference 7. 2816 J. Chem. Phys., Vol. 113, No. 7, 15 August 2000 Saiz, Gua `dia, and Padro ´ This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.83.95.18 On: Wed, 14 Jan 2015 14:53:57 length is plotted in Fig. 1. This picture shows that the MD simulations are long enough to reach a plateau value for ⑀ and that this plateau value is quite sensitive to the initial conditions. The value obtained by averaging the two runs, i.e., ⑀ ⫽16⫾2, is lower than the experimental data, but the agreement is similar to that achieved for liquid methanol.2–4 Since part of the discrepancies are due to the value of the high-frequency dielectric constant, we find it more significant to compare the difference ( ⑀ ⫺ ⑀ ⬁). In this case, the agreement between experimental and MD simulation results is better. It should be noticed that the experimental value of ⑀ ⬁corresponds to a polarizable system 共i.e., the real fluid兲, whereas the simulations have been carried out with a system of nonpolarizable molecules for which ⑀ ⬁⫽1. B. Wave-vector dependent dielectric constant The other approach to the computation of the static dielectric constant is from the low klimit of the longitudinal and transverse components of the permittivity tensor. In this study, kmin⫽0.273Å⫺1. We have calculated the wave-vector dependent quantities for the three smallest kvalues, namely, k⫽kmin ,&kmin , and )kmin , as well as for k⫽ 冑 6kmin and 3kmin . The values of xLand xTfound from the MD simulations by using Eqs. 共15兲and 共16兲, the resulting ⑀ L(k) and ⑀ T(k) components of the permittivity tensor computed from Eqs. 共10兲and 共11兲and the ⑀ R(k) values computed from Eq. 共14兲are collected in Table II. The statistical errors have been evaluated by using the blocking method as described in Ref. 23. The general trends found in the kdependence of the ⑀ T(k) and ⑀ L(k) functions agree with previous theoretical results for model polar liquids.24–26 In the low klimit, we can see that the relation ⑀ T(kmin)⯝ ⑀ ⯝ ⑀ R(kmin) is satisfied within the statistical errors. The good agreement between the results obtained following the two approaches to the calculation of ⑀ is an indication that the size of the simulated system is sufficiently large to be considered as representing a real isotropic fluid. As was already pointed out by Fonseca and Ladanyi,2 ⑀ T(k) provides a route much better than ⑀ L(k)to determine the dielectric permittivity. While ⑀ T(k) is a slowly varying function of knear k⫽0, the transition ⑀ L(k→0) ⯝ ⑀ occurs very sharply. Furthermore, ⑀ L(k) has a singularity at xL(k)⫽1/yand it is negative for yxL(k)⬍1关see Eq. 共10兲兴. As we can see in Table II, we have been able to obtain ⑀ L(k) for values of klower and higher than that corresponding to this singularity. To the best of our knowledge the only previous computer simulation results showing the divergence of ⑀ L(k) at small kare those reported for liquid water in a recent work of Bopp et al.27 V. DYNAMIC DIELECTRIC PROPERTIES In this section we are concerned with the calculation of the dipole density time correlation functions and their power spectra. Special attention has been dedicated to the investigation of the origin of the different contributions to the total dipole moment autocorrelation function. A thorough comparison with experimental dielectric relaxation data, where collective motions are of special interest, has been also performed. A. Dipole density time correlation functions The normalized time correlation functions of the longitudinal ⌽L(k,t) and transverse ⌽T(k,t) components of the dipole density defined by Eq. 共8兲have been computed during the MD simulations for the five kvalues listed above. In Figs. 2 and 3, the ⌽L(k,t) and ⌽T(k,t) functions are plotted and a detail of the short-time regime is provided in the inset. The ⌽(t) function resulting from the MD simulations, as well as its different contributions are represented in Fig. 4. The relaxation of the longitudinal dipole density plays a central role in solvation dynamics.28 As can be seen from Fig. 2, the main feature of ⌽L(k,t) is a fast initial decay with a strongly oscillatory character. The same behavior was previously observed for methanol4and water.29 The initial decay of ⌽L(k,t) has been attributed to inertial and librational dynamics.30 It is worthy to note that other polar nonhydrogen bonded liquids show a fast decay of ⌽L(k,t) to zero and then overdamped oscillations.16 TABLE II. Longitudinal and transverse components of the wave-vector dependent static dielectric constant. (k/kmin)2k/Å⫺1xL(k)xT(k) ⑀ L(k) ⑀ T(k) ⑀ R(k) 1 0.273 0.0581⫾0.0008 0.93⫾0.45 10.0⫾1.2 15.5⫾0.7 16.1⫾0.5 2 0.386 0.0605⫾0.0003 0.80⫾0.01 15.9⫾1.2 13.3⫾0.1 13.2⫾0.1 3 0.473 0.0656⫾0.0004 0.65⫾0.02 ⫺62.3⫾24.1 10.1⫾0.3 9.98⫾0.2 6 0.669 0.0960⫾0.0006 0.490⫾0.005 ⫺2.05⫾0.04 8.58⫾0.08 5.10⫾0.02 9 0.819 0.136⫾0.001 0.394⫾0.003 ⫺0.90⫾0.01 7.10⫾0.05 2.902⫾0.001 FIG. 2. Normalized time correlation function of the longitudinal component of the dipole moment density for (k/kmin)2⫽1共full line兲,(k/kmin)2⫽2共dot- ted line兲,(k/kmin)2⫽3共dashed line兲,(k/kmin)2⫽6共long dashed line兲, and (k/kmin)2⫽9共dotted–dashed line兲. The inset shows the short-time behavior of ⌽L(k,t). 2817J. Chem. Phys., Vol. 113, No. 7, 15 August 2000 Dielectric properties of liquid ethanol This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.83.95.18 On: Wed, 14 Jan 2015 14:53:57 Comparison of Figs. 2 and 3 reveals important differences in the relaxation of the transverse and longitudinal components of the dipole density. The decay of the transverse component is slower than that of its longitudinal counterpart. On the other hand, the behavior of the ⌽T(k,t) functions is similar to that of the total dipole moment autocorrelation function ⌽(t)共see Figs. 3 and 4兲. For the transverse and total functions displayed in Figs. 3 and 4, it is possible to distinguish three separate regimes. The long-time regime which takes place for t⬎1.5ps can be modeled as an exponential function and describes the relaxation of the principal dielectric band of ⌽T(k,t) and ⌽(t). This regime takes place when inertial and librational dynamics do not contribute appreciately to the collective reorientation. The short-time regime, i.e., for t⬍0.5ps, which is shown in detail in the inserts of Figs. 3 and 4, has a fast initial inertial decay and then an oscillatory behavior which is related to the librational dynamics of the molecule. The oscillatory behavior is characteristic of hydrogen bonded liquids and it is absent in polar liquids without hydrogen bonding interactions such as acetonitrile.16 The fast initial decay and the oscillatory behavior are both mainly due to the individual particle reorientation. This aspect will be discussed in more detail in Sec. VB. The intermediate regime can also be modeled as an exponential function but with a characteristic time smaller than the one corresponding to the long-time regime. The large and intermediate regimes for the different ␺ functions, where ␺ is any of the ⌽(t)or⌽T(k,t) correlation functions, have been fitted to a function of the type ␺ 共t兲⬃兺 j⫽1 2 Ajexp共⫺t/ ␶ j兲,t⬎0.5 ps. 共26兲 For the whole time interval considered, these functions can be expressed by ␺ 共t兲⬃兺 j⫽1 2 Ajexp共⫺t/ ␶ j兲⫹ 冉 1⫺兺 j⫽1 2 Aj 冊 ␾ librational-inertial共t兲. 共27兲 The fitted parameters are summarized in Table III. We observe that the amplitude A1of the main relaxation process decreases with increasing k, while that of the second process A2increases. Both relaxation times decrease as kincreases, but the dependence on kis much more important in the case of ␶ 1than for ␶ 2. On the other hand, the relaxation times corresponding to ⌽T(kmin ,t) are smaller than those of ⌽(t). The difference is specially significant in the case of ␶ 1. The so-called Debye relaxation time ␶ Dis also reported in Table III. It was obtained from the following expression:31 ␶ D⫽lim ␻ →0 ⑀ ⫺ ⑀ 共 ␻ 兲 i ␻ 关 ⑀ 共 ␻ 兲⫺ ⑀ ⬁兴.共28兲 For a Debye dielectric32 it should be ␶ D⯝ ␶ 1. According to the results given in Table III, liquid ethanol does not behave like a Debye fluid since ␶ Dis significantly lower than ␶ 1. B. Self- and distinct contributions to the total dipole moment autocorrelation function Decomposing 具 M ជ (t)•M ជ (0) 典 into its components arising from the correlations of the dipole moment of a molecule and the correlations between distinct 共different兲molecular dipole moments, one obtains the following expression for ⌽(t)关see Eq. 共21兲兴 ⌽共t兲⫽N ␮ 2 具 兩 M ជ 共0兲 兩 2 典 ⌽s共t兲 ⫹N共N⫺1兲 ␮ 2 具 ␮ ជ i共0兲•u ជ j共0兲 典 具 兩 M ជ 共0兲 兩 2 典 ⌽d共t兲,共29兲 where ⌽s(t) is the reorientational autocorrelation function of the individual dipoles and can be expressed by ⌽s共t兲⫽ 具 ␮ ជ i共t兲• ␮ ជ i共0兲 典 具 兩 ␮ ជ i共0兲 兩 2 典 ⫽ 具 u ជ i共t兲•u ជ i共0兲 典 ,共30兲 FIG. 3. Same as Fig. 2 for the transverse component of the dipole moment density ⌽T(k,t). TABLE III. Two-exponential fit parameters for the dipole density time correlation functions resulting from MD simulations 关components: self, total, and transverse 共T兲兴, Debye relaxation time, and single-particle correlation time. Component (k/kmin)2A1 ␶ 1/ps A2 ␶ 2/ps ␶ D/ps ␶ s/ps Self 0 0.731 26.0 0.105 1.50 20 Total 0 0.871 60.5 0.075 1.60 49 T1 0.861 46.5 0.080 1.50 T2 0.834 32.8 0.105 1.30 T3 0.785 20.9 0.125 1.25 T6 0.693 12.9 0.181 1.20 T9 0.575 9.70 0.250 1.18 2818 J. Chem. Phys., Vol. 113, No. 7, 15 August 2000 Saiz, Gua `dia, and Padro ´ This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.83.95.18 On: Wed, 14 Jan 2015 14:53:57 with u ជ i(t) a unit vector along the molecular dipole moment, ␮ ជ i(t)⫽ ␮ u ជ i(t). The ⌽d(t) function is given by ⌽d共t兲⫽ 具 u ជ i共t兲•u ជ j共0兲 典 具 u ជ i共0兲•u ជ j共0兲 典 ,i⫽j,共31兲 and, using Eq. 共23兲, 具 u ជ i共0兲•u ជ j共0兲 典 ⫽Gk⫺1 N⫺1,i⫽j.共32兲 Thus, Eq. 共29兲can be written as ⌽共t兲⫽1 Gk⌽s共t兲⫹ 冉 1⫺1 Gk 冊 ⌽d共t兲.共33兲 In our system, the prefactors of the self and distinct contributions are (1/Gk)⫽0.34⫾0.02 and (1⫺1/Gk)⫽0.65 ⫾0.02, respectively. Thus, the main contribution to ⌽(t) corresponds to the correlations among different molecular dipoles. Figure 4 shows the two contributions to ⌽(t). In the inset of the same figure we show the short-time behavior of the total and self-autocorrelation functions. The selfcomponent for short times consists of a fast initial inertial decay with oscillations due to librational dynamics. For longer times, ⌽s(t) shows a nearly exponential decay. The long-time decay of ⌽s(t) is faster than that of the total dipole moment autocorrelation function ⌽(t). The distinct component for short times does not present a fast inertial decaying behavior. Yet, ⌽d(t) decays more slowly than ⌽s(t). Therefore, the behavior of ⌽(t) at long times will be dominated by ⌽d(t). Thus, the long-time behavior of ⌽(t) should be associated with collective motions. As expected, the relaxation times ␶ 1and ␶ 2found for ⌽s(t) according to Eq. 共26兲are smaller than those corresponding to ⌽(t)共see Table III兲. Also given in Table III is the single-particle correlation time ␶ sdefined as ␶ s⫽ 冕 0 ⬁ ⌽s共t兲dt.共34兲 The relation between the total and single-particle correlation times is interesting from a theoretical point of view. According to our findings, liquid ethanol approximately obeys the Kivelson–Madden relationship14 ␶ D/ ␶ s⯝Gk. C. Dielectric relaxation Dielectric relaxation of hydrogen bonded liquids has been extensively studied by experimental techniques. Data for liquid ethanol at room temperature are available in the microwave zone 共up to 3 cm⫺1兲of the spectra6and up to frequencies of THz 共34 cm⫺1兲reaching the far infrared region.7The experimental dielectric relaxation provides information on the real and imaginary parts of the frequency dependent permittivity ⑀ ( ␻ )⫽ ⑀ ⬘( ␻ )⫺i ⑀ ⬙( ␻ ). In the case of monohydric alcohols, experimental data are usually analyzed by considering an empirical model which consists of three Debye-type processes.33 With this assumption, ⑀ 共 ␻ 兲can be separated into the real and imaginary parts ⑀ ⬘共 ␻ 兲⫽ ⑀ ⬁⫹共 ⑀ ⫺ ⑀ ⬁兲兺 j⫽1 ngj 关1⫹共 ␻ ␶ j兲2兴共35兲 and ⑀ ⬙共 ␻ 兲⫽共 ⑀ ⫺ ⑀ ⬁兲兺 j⫽1 ngj共 ␻ ␶ j兲 关1⫹共 ␻ ␶ j兲2兴,共36兲 respectively. In Eqs. 共35兲and 共36兲, ␶ jare the relaxation times and gj⫽( ⑀ j⫺ ⑀ ⬁j)/( ⑀ ⫺ ⑀ ⬁), ⑀ ⬁j⫽ ⑀ j⫹1, are the weights of processes jcontributing to the total dispersion, beginning at the static permittivity ⑀ ⫽ ⑀ 1and ending at the extrapolated high frequency permittivity ⑀ ⬁ ⫽lim ␻ →⬁ ⑀ ⬘( ␻ ). The two available sets of experimental dielectric relaxation parameters for liquid ethanol at room temperature are listed in Table IV. It is worthy to note that Kindt and Schmuttenmaer7obtained significantly faster relaxation times for the second and third Debye processes than those reported by Barthel and co-workers.6The ␶ jrelaxation times are comparable to the values reported in Table III for the total dipole moment autocorrelation function. In general, the relaxation times resulting from MD are smaller than the experimental estimations. The disagreement is especially significant in the case of the main dispersion ␶ 1value. A similar discrepancy was observed in the case of methanol.3,4 On the other hand, the ␶ 2value resulting from ⌽(t) is quite close to the second relaxation time given in Ref. 7. According to our findings, the relaxation time corresponding to the third Debye process cannot be determined from the decay of ⌽(t). This is probably due to the very different order-of-magnitude of ␶ 3and ␶ 1. It should thus be determined from the decay of a function directly related to the physical origen of the corresponding relaxation process 共see, for instance, Ref. 34兲. FIG. 4. Normalized total dipole moment autocorrelation function ⌽(t) 共solid line兲, self-contribution Gk ⫺1⌽s(t)共dotted–dashed line兲, and distinct contribution (1⫺Gk ⫺1)⌽d(t)共long dashed line兲. The inset shows the shorttime behavior of ⌽(t)共full line兲and ⌽s(t)共dotted–dashed line兲. TABLE IV. Dielectric relaxation parameters obtained from fitting three superimposed Debye equations to the experimental ⑀ ⬘( ␻ ) and ⑀ ⬙( ␻ ) data of liquid ethanol at room temperature. g1g2g3 ␶ 1/ps ␶ 2/ps ␶ 3/ps Ref. 6 0.917 0.031 0.052 163 8.97 1.81 Ref. 7 0.901 0.064 0.035 161 3.3 0.22 2819J. Chem. Phys., Vol. 113, No. 7, 15 August 2000 Dielectric properties of liquid ethanol This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.83.95.18 On: Wed, 14 Jan 2015 14:53:57 The MD results for ⑀ ⬘( ␻ ), ⑀ ⬙( ␻ ) and the corresponding Cole–Cole plot are compared with the experimental data in Figs. 5 and 6, respectively. The deviation of the Cole–Cole plot from a semicircle confirms that liquid ethanol is not a Debye fluid. In the high frequency region of the Cole–Cole plot, MD results exhibit the typical loops characteristic of fast librational motions. The dependence of the permittivity with frequency is qualitatively reproduced by the MD results, but we can observe some quantitative disagreements between the simulation results and experimental data. These disagreements are partly due to the discrepancies in the values of the static permittivity and the permittivity at optical frequencies. In Figs. 7 and 8, this is overcome by plotting the quantities ( ⑀ ⬘( ␻ )⫺ ⑀ ⬁)/( ⑀ ⫺ ⑀ ⬁) and ⑀ ⬙( ␻ )/( ⑀ ⫺ ⑀ ⬁). The discrepancies are now effectively reduced. The disagreements in the decay of the ( ⑀ ⬘( ␻ )⫺ ⑀ ⬁)/( ⑀ ⫺ ⑀ ⬁) function as well as in the position of the peak of ⑀ ⬙( ␻ )/( ⑀ ⫺ ⑀ ⬁) are directly related to the small value of the first relaxation time ␶ 1resulting from MD. The discrepancies in the value of the second relaxation time ␶ 2originate the differences observed in the intermediate frequency region of the Cole–Cole plot. VI. INFRARED ABSORPTION The experimental infrared absorption spectrum of a system is given in terms of the infrared absorption coefficient ␣ 共 ␻ 兲and the imaginary part ⑀ ⬙( ␻ ) of the frequency dependent dielectric constant. These quantities are related to the absorption line shape, I( ␻ ), by35 ␣ 共 ␻ 兲⫽4 ␲ 2 ␻ 3Vបcn共 ␻ 兲关1⫺exp共⫺ ␤ ប ␻ 兲兴I共 ␻ 兲,共37兲 ⑀ ⬙共 ␻ 兲⫽c ␣ 共 ␻ 兲n共 ␻ 兲 ␻ ⫽4 ␲ 2 3Vប关1⫺exp共⫺ ␤ ប ␻ 兲兴I共 ␻ 兲, 共38兲 where ␤ ⫽(kBT)⫺1,ប⫽h/2 ␲ ,his the Planck constant, cis the speed of light in vacuum, n( ␻ ) is the refractive index of the medium, and I( ␻ ) is the Fourier transform of the unnormalized total dipole moment autocorrelation function I共 ␻ 兲⫽1 2 ␲ 冕 ⫺⬁ ⬁dt 具 M ជ 共t兲•M ជ 共0兲 典 exp共⫺i ␻ t兲.共39兲 For classical systems, 具 M ជ (t)•M ជ (0) 典 is real and symmetric, and I( ␻ ) can be written as I共 ␻ 兲⫽1 ␲ 冕 0 ⬁dt 具 M ជ 共t兲•M ជ 共0兲 典 cos共 ␻ t兲.共40兲 Equations 共37兲,共38兲, and 共40兲can be used to relate infrared data with MD results. In the classical limit (ប→0), Eq. 共38兲 has the following form: ⑀ ⬙共 ␻ 兲⫽c ␣ 共 ␻ 兲n共 ␻ 兲 ␻ ⫽共 ⑀ ⫺ ⑀ ⬁兲 ␻ 冕 0 ⬁dt⌽共t兲cos共 ␻ t兲, 共41兲 FIG. 5. Real ⑀ ⬘( ␻ ) and imaginary ⑀ ⬙( ␻ ) parts of the frequency dependent permittivity. Molecular dynamics results 共solid line兲, experimental data from Ref. 6 共dotted line兲and from Ref. 7 共dashed line兲. Experimental lines are calculated from the parameters of Table IV. The arrows indicate the vertical scale corresponding to each curve. FIG. 6. Cole–Cole plot of the frequency dependent permittivity. Molecular dynamics results 共solid line兲, experimental data from Ref. 6 共dotted line兲and from Ref. 7 共dashed line兲. Experimental lines are calculated from the parameters of Table IV. FIG. 7. Same as Fig. 5 with the real and imaginary parts of ⑀ 共 ␻ 兲scaled as ( ⑀ ⬘( ␻ )⫺ ⑀ ⬁)/( ⑀ ⫺ ⑀ ⬁) and ⑀ ⬙( ␻ )/( ⑀ ⫺ ⑀ ⬁), respectively. FIG. 8. Same as Fig. 6 with the real and imaginary parts of ⑀ 共 ␻ 兲scaled as ( ⑀ ⬘( ␻ )⫺ ⑀ ⬁)/( ⑀ ⫺ ⑀ ⬁) and ⑀ ⬙( ␻ )/( ⑀ ⫺ ⑀ ⬁), respectively. 2820 J. Chem. Phys., Vol. 113, No. 7, 15 August 2000 Saiz, Gua `dia, and Padro ´ This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.83.95.18 On: Wed, 14 Jan 2015 14:53:57 where ⌽(t) is defined by Eq. 共21兲and we have used the definition of ⑀ given in Eq. 共22兲. Since we are concerned with frequencies higher than that of the principal dispersion band, we neglect the frequency dependence of the refractive index36 to write ␣ 共 ␻ 兲⬀ ␻ ⑀ ⬙共 ␻ 兲⬀ ␻ 2Re共⌽共 ␻ 兲兲.共42兲 An alternative procedure to calculate ␣ 共 ␻ 兲is through the k →0 limits of ⑀ L(k, ␻ ) and ⑀ T(k, ␻ ). As in the ␻ ⫽0 case 共see Sec. IVB兲, the most straightforward route to ⑀ ⬙( ␻ )is provided by the transversal component. In the k→kmin limit, Eqs. 共2兲and 共6兲lead to ⑀ ⬙共 ␻ 兲⯝Im共 ⑀ T共kmin , ␻ 兲兲⬀ ␻ Re共⌽T共kmin , ␻ 兲兲.共43兲 The absorption coefficient is thus proportional to ␻ 2Re(⌽T(kmin , ␻ )). Due to the ␻ 2factor, ␣ 共 ␻ 兲emphasizes the high frequency part of the spectra. Resonant motions lead to maxima in the absorption coefficient whereas exponential relaxations result in slowly increasing functions that reach a constant value for frequencies larger than the inverse of the characteristic time. The MD results for ␻ 2Re(⌽T(kmin , ␻ )) and ␻ 2Re(⌽( ␻ )) are depicted in Figs. 9 and 10, respectively. In Fig. 9, the ␻ 2Re(⌽L(kmin , ␻ )) function is also plotted for the sake of completeness. The ␻ 2Re(⌽T(kmin , ␻ )) spectrum exhibits a broad maximum at frequencies located between 500 and 750 cm⫺1with a peak at approximately 640 cm⫺1. Similar characteristics are shown by the ␻ 2Re(⌽L(kmin , ␻ )) function in the same region. However, the maximum is located at slightly higher frequencies 共⬃690 cm⫺1兲. On the other hand, the very good agreement between ␻ 2Re(⌽T(kmin , ␻ )) and ␻ 2Re(⌽( ␻ )) corroborates that the k→0 limit has effectively been reached in our simulations. To the best of our knowledge, there are not experimental infrared data to compare with our results. Experimental data corresponding to liquid methanol display a maximum at similar frequencies 共around 570 cm⫺1according to Ref. 37 and around 700 cm⫺1accord- ing to Ref. 38兲. To gain insight into the dynamic origin of the spectral features, we have analyzed the self- and distinct contributions to ␣ 共 ␻ 兲. By using Eq. 共33兲, we obtain the following relation ␻ 2⌽共 ␻ 兲⫽1 Gk ␻ 2⌽s共 ␻ 兲⫹ 冉 1⫺1 Gk 冊 ␻ 2⌽d共 ␻ 兲.共44兲 As can be seen from Fig. 10, the self-contribution Gk ⫺1 ␻ 2Re(⌽s( ␻ )) has a broad band in the infrared region located between 400 and 800 cm⫺1, with a maximum at ⬃600 cm⫺1. This frequency corresponds to an oscillation with a characteristic time of 0.05 ps. This high frequency maximum can be assigned to librational motions of the hydroxyl H atoms around the central methyl–oxygen bond and, more specifically, to hydrogen bonded molecules. This assumption is supported by a previous analysis of the spectral densities in the liquid.11 This analysis shows that the power spectrum of the velocity autocorrelation function of H-bonded hydrogen atoms has a band whose maximum is around 630 cm⫺1, while the H atoms that do not participate in a hydrogen bond give a band whose maximum is around 370 cm⫺1. Because of the very low percentage of nonbonded molecules,11 this band at intermediate frequencies can not be appreciated in Fig. 10. The Gk ⫺1 ␻ 2Re(⌽s( ␻ )) function also presents a maximum at low frequencies in the far infrared zone of the spectra. This peak is located at about 70 cm⫺1. It corresponds to an oscillation with a characteristic time ten times greater than the previous one 共⬃0.5 ps兲. In the case of water, the peak found in this frequency range has been ascribed to different physical origins, namely, hindered motions of the molecules in the cage formed by its neighbors,39 induced dipole contributions,40 and bending vibrations of the hydrogen bonds 共O¯O¯O units兲.41 In our ethanol simulations, dipole induced contributions have not been taken into account. On the other hand, a maximum at about 70 cm⫺1also appears in the power spectra of the velocity autocorrelation functions of the oxygen atoms of liquid ethanol, and this maximum is independent of the hydrogen bonding state of the molecule.11 We may then conclude that the low frequency peak shown by Gk ⫺1 ␻ 2Re(⌽s( ␻ )) can be ascribed to hindered motions of the molecules in the cage formed by its neighbors. The distinct component (1⫺Gk ⫺1) ␻ 2Re⌽d( ␻ )共see Fig. 10兲presents no spectral features until frequencies between 600 and 800 cm⫺1with a maximum at ⬃670 cm⫺1. Thus, the contribution of the distinct component produces a shift to higher frequencies of the librational band of the absorption FIG. 9. ␻ 2times the real part of the transverse 共solid line兲and longitudinal 共dotted line兲components of the dipole moment densitity spectrum for k ⫽kmin . Notice that the curve corresponding to ⌽L(kmin , ␻ ) has been multiplied by 0.1. FIG. 10. ␻ 2times the real part of the power spectrum of the total dipole moment autocorrelation function 共solid line兲, self-contribution Gk ⫺1 ␻ 2⌽s( ␻ )共dotted–dashed line兲, and distinct contribution (1 ⫺Gk ⫺1) ␻ 2⌽d( ␻ )共dashed line兲. 2821J. Chem. Phys., Vol. 113, No. 7, 15 August 2000 Dielectric properties of liquid ethanol This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 147.83.95.18 On: Wed, 14 Jan 2015 14:53:57