Article pubs.acs.org/JPCA
Dielectric and Terahertz Spectroscopy of Polarizable and Nonpolarizable Water Models: A Comparative Study M. Sega* and C. Schröder Department of Computational Biological Chemistry, University of Vienna, Währinger Strasse 17, 1090 Wien, Austria S Supporting Information *
ABSTRACT: Using extensive classical molecular dynamics simulations, we compute the dielectric and far-infrared spectra of nine popular water models, including polarizable and nonpolarizable ones. We analyze the dielectric spectra using a two-relaxation model that allows one to extract the characteristic time of both the main dielectric relaxation and the fast relaxation. The use of a Cole−Cole functional form permits also quantitative assessment of the absence of deviations from the Debye form of the main dielectric peak. In the THz region of the spectrum, we compute the infrared absorbance caused by molecular libration, which appears to be qualitatively different for three main groups of molecular models. The complexity of the librational band is further investigated by decomposing the spectrum into the contributions of water fractions with a different number of hydrogen-bonded neighbors.
■
“broad-band dielectric spectroscopy”.10 The need for probing a wide frequency range, up to 15 decades (reaching the farinfrared part of the electromagnetic spectrum), is usually related both to the dielectric peaks being very broad with respect to other spectroscopic techniques and to the sensitivity of the resonant frequency to the temperature. For computer simulations, this means that particularly long sampling times are in principle needed to obtain a well-resolved spectrum, but often, dielectric spectra can be characterized by a relatively small number of parameters (relaxation times and amplitudes). For this reason, a considerable number of computed spectra reported in literature are the result of a fit of the collective dipole autocorrelation function to a predefined set of analytical correlation functions, and a growing technological effort in this direction is ongoing.11 However, if even small biases have to be avoided, in order to produce model-independent data, it is necessary to compute the dielectric spectrum directly from the sampled dipole autocorrelation functions. From the early 1970s to the present time, a large number of computer simulation investigations have been dedicated to the dielectric properties of water and in several cases to its dielectric spectrum, of which we can only cite a few representative examples,12−17 raising no claim to completeness. Some studies have also presented a comparative analysis of the dielectric spectrum of several nonpolarizable18 as well as polarizable models,19 but the dipole autocorrelation function relaxation has always been assumed to be described by an exponential decay.
INTRODUCTION The investigation of dielectric spectra is a fundamental step toward the interpretation of the collective molecular dynamics of liquids and liquid mixtures, and in recent years, a growing number of systems with increasing complexity have been investigated with the aid of molecular dynamics simulations in order to establish a connection between the microscopic dynamics and the resulting dielectric spectrum. Water mixtures represent an important fraction of the systems of interest, and the precise knowledge of the dielectric spectrum of bulk water (for the different computational models available) is an important prerequisite for understanding how the solvent− solute interactions influence spectral features of the solution. Water is a peculiar molecular liquid in many respects; the strong directional attraction provided by hydrogen bonds plays an important role in determining several thermodynamic, dynamic, and structural anomalies1 such as its temperature of maximum density2 or its characteristic shrinking upon melting up to more subtle effects such as the surface tension inflection.3 The spectral properties of water are also atypical. Water is also probably the only known liquid whose color originates only from molecular vibrations4,5 instead of a combination of vibrational, rotational, and electronic excitations (as is the case of most dyes6). Electronic transitions in water, on the other hand, are active only in the ultraviolet region of the spectrum (the lowest, nondissociative states in small cluster having an energy in the range of 6−9 eV7) and can be safely ignored for frequencies in the far-infrared region or lower. However, a coupling between the translational degrees of freedom of the water molecule and the electronic cloud is believed to give rise to the characteristic 200 cm−1 peak, which is in fact much suppressed in classical simulations.8,9 The collective dielectric response of a material is usually investigated using various techniques known collectively as © XXXX American Chemical Society
Special Issue: 25th Austin Symposium on Molecular Structure and Dynamics Received: July 24, 2014 Revised: August 13, 2014
A
dx.doi.org/10.1021/jp507419e | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
One should not forget that in the experimental measurements, water displays an almost perfect Debye-like main peak and that deviations from this behavior are a signature of additional relaxation processes. This happens, for example, when other substances are present in solution.20,21 The accuracy of experimental measurements can, for example, detect deviations on the order of 0.1% from the Debye behavior upon addition of small quantities (0.1 mol/L) of salt.21 The main aim of this work is to provide a reference for future simulation studies of the dielectric properties of aqueous solutions and mixtures, but we also investigate the influence of the number of hydrogen-bonded neighbors on some features of the infrared absorption spectrum. We report dielectric and THz spectra sampled over 100 ns trajectories for nine among the most popular water models, namely, SPC,22 SPC/E,23 TIP4P,24 TIP4P-2005,25 TIP5P,26 TIP5P-E,27 COS/G3,19 SWM4-DP,28 and SWM6.29
Figure 1. Schematic representation of a water molecule model: the oxygen (OW) and the two hydrogen (HW1, HW2) sites are present in all water models; the lone pair charges (LP1,LP2) are present only in the TIP5P, TIP5P-E, and SWM6 models; the virtual site (MW) is present in the TIP4P, TIP4P-2005, SWM4-DP, SWM6, and COS/G3 models. The Drude particle of the polarizable models is not shown. The values of the geometric and electric parameters are reported in Tables 1 and 2
■
COMPUTATIONAL DETAILS We have performed molecular dynamics simulations with a modified version of the gromacs 4.5.7 package30,31 that saves at every step the total electric current to disk. We have used an integration time step of 1 fs and simulated the isobaric− isothermal ensemble at temperature T = 298 K (β = 1/(kBT) = 2.47638 kJ/mol, where kB denotes Boltzmann’s constant) and at a pressure of 1 atm using the Nosé−Hoover32,33 and Parrinello−Rahman34 algorithms with time constants of 2 and 5 ps, respectively. The structure of water molecules has been kept rigid using the SETTLE algorithm.35 This removes internal vibrations, which are at wavenumbers larger than 1500 cm−1 and thus have little influence on the librations, which we investigate here. The position of the mobile Drude (or shell) particles generating the induced dipole moment is determined at every step using a steepest descent energy minimization with a tolerance of 0.1 kJ/mol/nm on the force acting on that particle. Simulation boxes containing approximately 4000 water molecules were prepared for each of the nine water models. For the COS/G3 model, a tolerance of 0.1 kJ/mol/nm was not achievable with the single-precision version of the gromacs code, most probably due to the large polarization charge of the model. For this model only, we used the double-precision version, which allowed us to reach the target accuracy. Note also that for this model, the estimated static dielectric permittivity, ϵ ≃ 107, is fairly large, especially if compared with the result in the original publication, where, however, the electrostatic interaction has been computed using the reaction field approach.19 The original potentials were parametrized using different combinations of electrostatic force treatment (for example, Ewald sum versus reaction field), the cutoff scheme, and longrange correction for Lennard-Jones interactions. We decided here to use a unified approach by truncating Lennard-Jones interactions at a distance of 1 nm and by calculating the electrostatic interaction using the smooth variant of the particle Ewald mesh method,36 with tinfoil boundary conditions, a mesh spacing not larger than 0.15 nm, a third-order polynomial interpolation, and the relative strength of the direct part of the potential at the cutoff of 10−5. This choice has been used for all water models, independent of the setup used in the original publications. In Figure 1, we report a schematic model that is general enough to encompass the characteristics of all nine models
under investigation. The oxygen (OW) and two hydrogens (HW1, HW2) are the only atoms used in the SPC and SPC/E models. These three atoms define a plane in which the additional charge (MW) for the TIP4P, SWM4-DP, and SWM6 models is located. The lone pair charges (LP1, LP2) are located in the plane orthogonal to that of the other charges and are present only in the TIP5P, TIP5P-E, and SWM6 models. MW, LP1, and LP2 are virtual sites having no real dynamics, and their position is determined at every step in order to fulfill the geometrical prescriptions of the model. The polarizable models have an additional Drude particle attached to the oxygen atoms (SWM4-DP,SWM6) or to the virtual site MW (COS/G3). The geometrical parameters of each model, as well as the value of the different partial charges, are listed in Tables 1 and 2. The dielectric spectrum can be obtained either from the collective dipole moment correlations ϵ(ω) = ϵ∞ +
4πβ {⟨P2⟩ + iω 3V
∫0
∞
e iωt ⟨P(0) ·P(t )⟩ dt } (1)
where the angular frequency is ω = 2πν or from the current autocorrelation function using the expression ϵ(ω) = ϵ∞ +
4πβ i 3ωV
∫0
∞
e iωt ⟨J(0) ·J(t )⟩ dt
(2)
where ϵ∞ is usually set to 1 for the nonpolarizable models and to the respective Clausius−Mossotti value ϵ∞ = 4πα/3vm for the polarizable one, where α is the molecular polarizability and vm is the molecular volume. Equations 1 and 2 are equivalent if no free charges are present in the system. The latter is more compact, easier to implement as no information on the molecular structure is needed (e.g., to unfold the atoms of molecules that are crossing the simulation box boundaries), and, in addition, more general than its dipolar counterpart as it can readily be used also in the presence of charged species. Often, however, particles with no real dynamics are employed to model the molecules, as is the case of virtual sites and Drude particles. In these cases, the current associated with the site is not directly available, and eq 1 is the only feasible option to compute the dielectric spectrum. B
dx.doi.org/10.1021/jp507419e | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
Table 1. Parameters for the Nonpolarizable Water Models qOa qHa qMa qSa qLPa dO−Hb dO−Mb dLP−Ob θH−O−Hc θLP−O−LPc σOb ϵOd a
SPC
SPC/E
TIP4P
TIP4P-2005
TIP5P
TIP5P-E
−0.82 0.41 − − − 0.1 − − 109.47 − 0.316557 0.650194
−0.8476 0.4238 − − − 0.1 − − 109.47 − 0.316557 0.650194
0.0 0.52 −1.04 − − 0.09572 0.015 − 104.52 − 0.315365 0.64852
0.0 0.5564 −1.1128 − − 0.09572 0.015 − 104.52 − 0.31589 0.7749
0.0 0.241 − − −0.241 0.09572 − 0.07 104.52 109.47 0.312 0.66944
0.0 0.241 − − −0.241 0.09572 − 0.07 104.52 109.47 0.30876 0.75838
Units of elementary charge. bUnits: nm. cUnits: deg. dUnits: kJ/mol.
produced a noticeable degradation of the signal because the value of the autocorrelation function was close to zero, with a corresponding decrease of the signal-to-noise ratio. For this reason, we used only the first 2.5 ns of the autocorrelation function to obtain the dielectric spectrum (see eq 1), again using a fast Fourier transform. Even applying this procedure, the high-frequency statistical noise at large time delays in ⟨P(0) P(t)⟩ induced a large, rapid fluctuation of the spectral signal after approximately 80 GHz. We decided to apply (only for the part of the spectrum with ν > 80 GHz) a simple low-pass filter by performing averages over progressively increasing block sizes, reducing in this way also the density of points to be displayed in the double-logarithmic plots.
Table 2. Parameters for the Polarizable Water Models qOa qHa qMa qSa qLPa dO−Hb dO−Mb dLP−Ob θH−O−Hc θLP−O−LPc σOb ϵ Od αe a
COS/G3
SWM4-DP
SWM6
0.0 0.450672 7.098656 −8.0 − 0.1 0.015 − 109.47 − 0.31746 0.94449 1.25
1.71636 0.55733 −1.11466 −1.71636 − 0.09572 0.024034 − 104.52 −
1.91589 0.53070 −1.13340 −1.62789 −0.10800 0.09572 0.0247 0.0315 104.52 101.098 0.319838 0.677808 1.04252
b
0.87998 c
■
MAIN DIELECTRIC AND FAST RELAXATION PEAKS The spectra resulting from the analysis are shown in Figure 2 (a) SPC, (b) SPC/E; Figure 3 (a) TIP4P, (b) TIP4P-2005; Figure 4 (a) TIP5P, (b) TIP5P-E; and Figure 5 (a) SWM4-DP, (b) SWM6, and (c) COS/G3. The frequency ranges from ν = 0.5 to 25 THz, and the spectra are decomposed into their real (ϵ′, squares) and imaginary (ϵ″, circles) parts. The qualitative, overall appearance of the spectra, shared by all models, encompasses a well-defined main dielectric peak in the vicinity of ν = 15 GHz (τ ≃ 10 ps), a second, weak relaxation at around ν = 3 THz (τ2 ≃ 0.1 ps, about 20 cm−1), and a third clear process in the far-infrared region approximately located at ν = 15 THz (10 fs, about 500 cm−1), associated with molecular
d
Units of elementary charge. Units: nm. Units: deg. Units: kJ/mol. Units: 10−3nm3
e
After 2ns of equilibration, during which the pressure and temperature relaxed to their equilibrium values, we sampled the collective dipole moment P(t) at each integration time step for 100 ns. We computed the autocorrelation function ⟨P(0)P(t)⟩ for the complete 100 ns trajectory using a fast Fourier transform approach37−39 (see the Supporting Information for the correlation functions for all models). The resulting autocorrelation functions were very well resolved up to about 2−3 ns. After this time, even tiny statistical fluctuations
Figure 2. Spectra of the three-site water models: SPC (a) and SPC/E (b). Real (squares) and imaginary (circles) parts are shown, together with the fitting to eq 3. The dotted line represents the imaginary part of the second relaxation contribution. Gray rhombuses and triangles are the experimental points taken from ref44 and represent the imaginary part of the spectrum and the residual imaginary part after subtracting the main peak contribution, respectively. C
dx.doi.org/10.1021/jp507419e | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
Figure 3. Spectra of the four-site water models: TIP4P (a) and TIP4P-2005 (b). The meaning of the symbols is the same as that in Figure 2.
Figure 4. Spectra of the five-site water models: TIP5P (a) and TIP5P-E (b). The meaning of the symbols is the same as that in Figure 2.
Figure 5. Spectra of the four-site water models: TIP4P (a) and TIP4P-2005 (b). The meaning of the symbols is the same as that in Figure 2.
librations.40 The translational peak observed experimentally at about 5 THz (200 cm−1) is not reproduced by the simulation, as was first noticed by Neumann,14 and the possible origin of its nature is discussed later on; it is worth noticing here that even polarizable models fail in reproducing this feature. The
intermediate, small relaxation process has been the subject of several studies and is related to the collisional relaxation of the non-hydrogen-bonded structure41 or to the relaxation time of water dipoles between two successive large-angle jumps.42 Early experimental estimates of τ2 ≃ 1 ps43 were later found to be D
dx.doi.org/10.1021/jp507419e | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
Table 3. Best-Fit Results Δϵ SPC SPC/E TIP4P TIP4P-2005 TIP5P TIP5P-E COS/G3 SWM4-DP SWM6 a
63.3 68.5 49.3 56.8 90.1 99.2 103.6 76.0 68.0
± ± ± ± ± ± ± ± ±
τ 0.4 0.7 0.4 0.5 1.2 1.4 0.8 0.5 1.1
6.6 11.0 6.6 11.7 15.4 17.4 11.7 10.9 12.9
± ± ± ± ± ± ± ± ±
γ 0.1 0.1 0.0 0.1 0.2 0.3 0.1 0.1 0.3
0.981 1.021 0.991 0.992 1.00 0.995 0.980 0.985 0.991
± ± ± ± ± ± ± ± ±
Δϵ2 0.005 0.008 0.006 0.008 0.010 0.011 0.007 0.006 0.014
0.7 1.5 0.7 1.0 1.2 1.3 1.2 1.2 1.1
± ± ± ± ± ± ± ± ±
τ2
0.4 0.4 0.3 0.3 0.6 0.7 0.5 0.3 0.6
0.05 0.34 0.11 0.26 0.34 0.39 0.13 0.15 0.16
± ± ± ± ± ± ± ± ±
ϵinf 0.05 0.16 0.08 0.14 0.33 0.37 0.10 0.08 0.17
1.67 1.98 1.36 1.78 2.29 2.96 2.65 1.70 2.83
± ± ± ± ± ± ± ± ±
ϵ∞a 0.47 0.23 0.25 0.22 0.41 0.43 0.47 0.29 0.58
1 1 1 1 1 1 1.6 1.4 1.5
Determined from the ensemble average.
biased by a insufficient coverage of the frequency domain, the more recent estimate at 298.15 K being τ = 8.32 ps and τ2 = 0.42 ps.44 To provide a quantitative analysis of the main dielectric peak, we performed a simultaneous fit (in the range of 10−4 < ν < 3 THz) of the real and imaginary parts of the measured spectra to the empirical functional form ϵ(ν) = ϵ inf +
Δϵ2 Δϵ + 1 + (i2πντ )γ 1 + (i2πντ2)
Table 4. Static Dielectric Permittivity ϵa SPC SPC/E TIP4P TIP4P-2005 TIP5P TIP5P-E COS/G3 SWM4-DP SWM6
(3)
Equation 3 describes two processes; the first is a Cole−Cole process with relaxation time τ, characterized by a symmetrical deviation from the Debye form (with γ = 1 corresponding to the Debye case), and the second, with relaxation time τ2, is a plain Debye process. Notice that ϵinf is a phenomenological offset, which differs from ϵ∞. The difference ϵinf − ϵ∞ takes into account the amplitude of all of the processes at frequencies higher than those considered in the fit, that is, at ν > 3 THz. We have chosen the Cole−Cole relaxation in order to have a clear, quantitative estimate of how well the models are reproducing the Debye form of the main dielectric peak, as the parameter γ measures the deviation from it. In principle, one could use a Cole−Cole relaxation to model also the faster process, but the introduction of an additional degrees of freedom turned out to make the fit unstable. However, the addition of the second process was fundamental to achieve a good fit of the spectrum and could not be neglected. Without entering the debate on the interpretation of the Cole−Cole functional form as an approximation for a superposition of a continuum distribution of a Debye relaxation process or as an expression of a fundamentally different type of relaxation, here we limit ourself to notice that it is a very convenient empirical formula to summarize the (symmetric) deviations from the Debye process, and as such, we use it to characterize more in detail the features of the main peak at around ν = 15 GHz. As free parameters for the Levenberg− Marquardt45−47 best-fit procedure, we used the exponent γ, the intensities of the relaxation processes Δϵ and Δϵ2, the characteristic times τ and τ2, and the offset ϵinf. The results of the fitting procedure along with the respective error estimate are presented in Table 3, while the resulting static dielectric permittivities are presented in Table 4.
■
a
65.7 72.0 51.4 59.6 93.6 103.5 107.5 78.9 71.9
± ± ± ± ± ± ± ± ±
1.3 1.3 1.0 1.0 2.2 2.5 1.8 1.0 2.3
ϵ = Δϵ + Δϵ2 + ϵinf.
the simulation data in terms of the infrared absorbance for this specific band. At such short wavelengths, quantum effects cannot be neglected.8,48 However, it is possible to use one of the several available classical limits of the exact formulas.49,50 We report here two approximations, namely, the so-called standard approximation, in which the absorption spectrum takes the form ⎡ 4πω ⎤ 1 − exp( −β ℏω) ASA (ω) = ⎢ ⎥ ⎣ 3V ℏcn(ω) ⎦ 1 + exp( −β ℏω) ∞
×
∫−∞ eiωt ⟨P(0)P(t )⟩ dt
(4)
and the harmonic approximation, in which the absorbance reads AHA (ω) =
β ℏω ⎡ 4πω ⎤ 1 − exp( −β ℏω) ⎢ ⎥ 2 ⎣ 3V ℏcn(ω) ⎦ 1 − exp( −β ℏω) ∞
×
∫−∞ eiωt ⟨P(0)P(t )⟩ dt
(5)
Here, ℏ is the reduced Planck constant, c is the speed of light, and n(ω) = [(ϵ′ + (ϵ′2 + ϵ″2)1/2)/2]1/2 is the frequencydependent refractive index.51 The harmonic approximation, which is exact in the harmonic regime,52 is the only approximation that satisfies the fluctuation−dissipation theorem and was found to perform generally better than the other approximations.50 Here, we report the computed spectra in both the standard and the harmonic approximation to provide a means of comparison. The computed absorbance spectra multiplied by n(ω), as a function of the wavenumber λ−1 = c/ν, are presented in Figure 6 (standard approximation) and Figure 7 (harmonic approximation) for all nine water models, normalized to the peak height and shifted along the absorbance axis for the sake of
INFRARED LIBRATIONAL BAND
At the far-infrared end of the dielectric spectrum, one can clearly see the typical libration peak, at around ν = 15 THz. This part of the spectrum is often probed with infrared spectroscopic techniques, and here, we provide an analysis of E
dx.doi.org/10.1021/jp507419e | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
region due to the close similarity of the models (which differ only by the values of the respective Lennard-Jones parameters), it is interesting to note that also the three polarizable models have remarkably similar spectra, despite their quite different geometrical and electrical characteristics. All spectra are rather complex, but those of the TIP5P and TIP5P-E models show a particularly clear split into two peaks. Using two Gaussian functions to model the TIP5P spectrum model leads to an almost perfect fit from 300 to 1200 cm−1 (dashed line) when the peaks are centered at 525 ± 3 and 864 ± 2 cm−1 in the standard approximation. In the harmonic approximation, on the contrary, the spectrum of the TIP5P model can be nearly perfectly fitted in the whole range by a superposition of three Gaussian functions centered at 59 ± 1, 653 ± 2, and 893 ± 2 cm−1.
Figure 6. Infrared absorbance spectra in the standard approximation: the results of the nine models are shifted (baselines shown as dotted lines) along the vertical axis, for the sake of clarity. From the bottom to the top: SPC; SPC/E; TIP4P; TIP4P-2005; TIP5P; TIP5P-E; COS/ G3; SWM4-DP; SWM6. Dashed line on the TIP5P curve: fit to the sum of two Gaussian functions.
■
CONTRIBUTION OF HYDROGEN-BONDED SUBNETWORKS TO THE LIBRATIONAL BAND In order to shed some light on the complex shape of the librational band, we now focus on the SPC/E model, which is the one that more closely resembles the experimental one, and analyze the role of the hydrogen bond network. Here, we consider two molecules to be hydrogen bonded if the distance between the two oxygens and that between one of the possible oxygen−hydrogen pairs is smaller than the distance of the first minima in the respective radial distribution function,55 namely, 0.335 nm for the oxygen−oxygen distance and 0.245 nm for the oxygen−hydrogen distance. We analyzed the contribution to the total spectrum coming from atoms that are hydrogen bonded to a number nHB (with nHB = 0−6) of neighbors. Because the correlation functions appearing in all previous expressions are collective quantities, the change in the number of hydrogen-bonded partners of one molecule does not pose any problem to the analysis of the correlation function itself (as in the case of single-particle properties). Following the linear response theory, in order to calculate the partial susceptibility of the fraction of molecules having nHB hydrogen-bonded neighbors, one needs to correlate the respective partial dipole moment PnHB with the quantity that couples to the external field, which is, in this case, always the total dipole moment P. This approach of separating different contributions to the total dipole moment (or current) can be particularly helpful in understanding the origin of different spectral features.56 Here, we compute the partial correlation functions ⟨PnHB(0)P(t)⟩ that are associated with the subnetworks of molecules having nHB hydrogen-bonded neighbors, and we use them to compute the respective contribution to the infrared spectrum in the harmonic approximation using the analogous of eq 5. In Figure 8, we report the contributions of the fractions of water molecules with different nHB to the infrared absorption spectrum in the harmonic approximation. The curves corresponding to nHB = 0, 1, and 6 have a very small amplitude due to the small relative population (see the inset of Figure 8) and are omitted. The total contribution is mainly determined by the nHB = 3 and 4 fractions, which have their maximum intensity at different wavenumbers. Even though water shows most of the time a coordination number of four, compatible with the tetrahedral structure of the network, the less connected water molecules still give a sizable contribution to the librational spectrum. However, to understand the effect of the number of hydrogen-bonded neighbors on the spectral
Figure 7. Infrared absorbance spectra in the harmonic approximation: the results of the nine models are shifted (baselines shown as dotted lines) along the vertical axis, for the sake of clarity. From the bottom to the top: SPC; SPC/E; TIP4P; TIP4P-2005; TIP5P; TIP5P-E; COS/ G3; SWM4-DP ;SWM6. Dashed line on the TIP5P curve: fit to the sum of three Gaussian functions.
clarity. While all of the models roughly reproduced the main qualitative characteristics of the dielectric peak, in the 600 cm−1 band, they present specific quantitative and qualitative differences that in some cases also strongly depend on the underlying approximation. In particular, the SPC, SPC/E, TIP4P, and TIP4P-2005 show the appearance of one main peak in the 600−700 cm−1 region with a pronounced tail at lower wavenumbers, which well reproduces qualitatively the experimental findings.53,54 An exception to this is represented by the translational band at 200 cm−1, which has been suggested to originate from the coupling of the electronic cloud and intermolecular motions and is therefore absent in classical simulations without any corrections.8 This coupling has however little or no influence8 on the librational band around 600 cm−1, which can thus be well reproduced by classical simulations. On the contrary to the SPC and TIP4P model families, the TIP5P and polarizable ones show more complex spectra, which have little experimental counterpart. They are characterized by a plateau at wavenumbers lower than 200 (TIP5P, TIP5P-E) and 400 cm −1 (polarizable models) in the standard approximation. In the harmonic approximation, this plateau is strongly damped but still survives in the form of a broad distribution. While it is not unexpected that the TIP5P and TIP5P-E models show a virtually identical spectrum in this F
dx.doi.org/10.1021/jp507419e | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
■
CONCLUSIONS
■
ASSOCIATED CONTENT
Article
We have performed a comparative analysis of the dielectric spectra of nine different water models with different complexity, ranging from rigid, nonpolarizable, three-site models (SPC, SPC/E) up to the more complex six-site polarizable model SWM6, aiming at providing a reference and guidance for future studies involving the simulated spectroscopy of water mixtures. The overall spectral features of the different models agree qualitatively in the presence of two clear relaxation processes, namely, the main dielectric relaxation at around ν = 15 GHz and a fast process in the region around ν = 0.8−2.4 THz. While the SPC, SPC/E, SWM4-DP, and SWM6 performed best in reproducing the static dielectric permittivity of water (ϵ = 78.457), the characteristic time of the main dielectric relaxation peak was best approximated by SPC, TIP4P, and SWM4-DP models. The fast relaxation time τ2 is not determined by the fitting procedure, with a precision high enough to allow a meaningful comparison, but was always found to be of the correct order of magnitude, with the exception of SPC water. All water models provided a main dielectric relaxation that is compatible within two standard deviations with a Debye process. Note that some of the measured values differ from those reported in the literature, even taking properly into account statistical inaccuracy due to different simulation conditions, for example, electrostatic force calculation, Lennard-Jones cutoff, and long-range corrections. The far-infrared absorption spectrum, instead, seems to be characterized by a superposition of different processes for all models. For the three- and four-site models, there is one clear peak with a long tail at lower frequencies that is qualitatively compatible with that in the experimental spectrum, while for the more complex models, a different structure is revealed. The TIP5P family shows two well-separated peaks, while the polarizable models, although much different from each other, show an unexpectedly similar far-infrared spectrum, especially due to the presence of the same plateau in the lower wavenumber region of the spectrum. One of the possible origins of the complex spectral features in the libration band is the presence of fractions of water molecules having a different number of hydrogen-bonded partners; we have investigated this issue and found that different fractions have a spread in their peak position that can reach 200 cm−1, with the largest contribution coming from the fractions with three and four hydrogen-bonded neighbors. The spectrum of each fraction, however, can in turn be described by a superposition of two Gaussian functions, indicating that other underlying processes are contributing to create the complex librational band. Investigations in this directions are ongoing.
Figure 8. Contribution (in the harmonic approximation) to the librational band of the infrared spectrum of molecules with a different number nHB of hydrogen-bonded neighbors. Dashed line: the total contribution. Inset: probability P(nHB) for a molecule to have nHB hydrogen-bonded neighbors.
properties of water, it is more convenient to divide the spectra by the respective population of the fraction, P(nHB), as shown in Figure 9. Here, it is seen that the curves follow a clear trend
Figure 9. Contribution (in the harmonic approximation) to the librational band of the different fractions of water molecules, divided by P(nHB) and normalized to the height of the peak for nHB = 5.
that sees the peak height and position increasing with nHB (apart from the height of the nHB = 6 fraction) and the width correspondingly narrowing. Most noticeably, the peaks shift by about 200 cm−1, going from nonbonded water molecules to the more constrained molecules with five or six neighbors. This behavior well fits the picture of an increasingly restricted rotational motion of water molecules in the presence of a high coordination number nHB. The present decomposition of the infrared spectrum, with the large spread in peak position observed for the different fraction, partially explains the complexity of the librational band. However, it has to be noted that the different contributions that we have calculated here cannot be described by a single process. Each of the contributions can be approximated rather well with two Gaussian functions, therefore showing that there are other levels of complexity in the librational spectrum that still need to be identified and are possibly related to the presence of two different infrared-active librations (in the molecular plane and out of the molecular plane, excluding the rotations around the dipole vector) or to the asymmetry between donors and acceptors.
S Supporting Information *
The time correlation function of the global dipole moment for all water models. This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest. G
dx.doi.org/10.1021/jp507419e | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
■
Article
(22) Berendsen, H. J.; Postma, J.; Van Gunsteren, W.; Hermans, J. Intermolecular Forces; Springer: Dordrecht, The Netherlands, 1981; pp 331−342. (23) Berendsen, H.; Grigera, J.; Straatsma, T. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (24) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926. (25) Abascal, J. L.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (26) Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910. (27) Rick, S. W. A Reoptimization of the Five-Site Water Potential (TIP5P) for Use with Ewald Sums. J. Chem. Phys. 2004, 120, 6085. (28) Lamoureux, G.; MacKerell, A. D., Jr.; Roux, B. A Simple Polarizable Model of Water Based on Classical Drude Oscillators. J. Chem. Phys. 2003, 119, 5185. (29) Yu, W.; Lopes, P. E.; Roux, B.; MacKerell, A. D., Jr. Six-Site Polarizable Model of Water Based on the Classical Drude Oscillator. J. Chem. Phys. 2013, 138, 034508. (30) Hess, B.; Kutzner, C.; van Der Spoel, D.; Lindahl, E. Gromacs 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (31) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; et al. Gromacs 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845− 854. (32) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (33) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695. (34) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182. (35) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the Shake and Rattle Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952−962. (36) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (37) Futurelle, R. P.; McGinty, D. J. Calculation of Spectra and Correlation Functions from Molecular Dynamics Data Using the Fast Fourier Transform. Chem. Phys. Lett. 1971, 12, 285−287. (38) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, 1st ed.; Oxford Science Publications; Clarendon Press: Oxford, U.K., 1987. (39) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992. (40) Draegert, D. A.; Stone, N.; Curnutte, B.; Williams, D. FarInfrared Spectrum of Liquid Water. J. Opt. Soc. Am. 1966, 56, 64−69. (41) Yada, H.; Nagai, M.; Tanaka, K. Origin of the Fast Relaxation Component of Water and Heavy Water Revealed By Terahertz TimeDomain Attenuated Total Reflection Spectroscopy. Chem. Phys. Lett. 2008, 464, 166−170. (42) Zasetsky, A. Y. Dielectric Relaxation in Liquid Water: Two Fractions Or Two Dynamics? Phys. Rev. Lett. 2011, 107, 117601. (43) Buchner, R.; Barthel, J.; Stauber, J. The Dielectric Relaxation of Water between 0 and 35 °C. Chem. Phys. Lett. 1999, 306, 57−63. (44) Fukasawa, T.; Sato, T.; Watanabe, J.; Hama, Y.; Kunz, W.; Buchner, R. Relation Between Dielectric and Low-Frequency Raman Spectra of Hydrogen-Bond Liquids. Phys. Rev. Lett. 2005, 95, 197802. (45) Levenberg, K. A Method for the Solution of Certain Problems in Least Squares. Q. Appl. Math. 1944, 2, 164−168. (46) Marquardt, D. W. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. J. Soc. Ind. Appl. Math. 1963, 11, 431−441.
ACKNOWLEDGMENTS The authors acknowledge support from the European Community’s Seventh Framework Programme (FP7-PEOPLE-2012-IEF) funded under Grant No. 331932 SIDIS and generous allocation of computer time on the Vienna Scientific Cluster (VSC-2).
■
REFERENCES
(1) Errington, J. R.; Debenedetti, P. G. Relationship Between Structural Order and the Anomalies of Liquid Water. Nature 2001, 409, 318−321. (2) Poole, P. H.; Sciortino, F.; Essmann, U.; Stanley, H. E. Phase Behaviour of Metastable Water. Nature 1992, 360, 324−328. (3) Sega, M.; Horvai, G.; Jedlovszky, P. Microscopic Origin of the Surface Tension Anomaly of Water. Langmuir 2014, 30, 2269. (4) Braun, C. L.; Smirnov, S. N. Why Is Water Blue? J. Chem. Educ. 1993, 70, 612. (5) Pope, R. M.; Fry, E. S. Absorption Spectrum (380−700 nm) of Pure Water. II. Integrating Cavity Measurements. Appl. Opt. 1997, 36, 8710−8723. (6) Zollinger, H. Color Chemistry: Syntheses, Properties, and Applications of Organic Dyes and Pigments, 3rd ed.; Wiley-VCH: Weinheim, Germany, 2001. (7) Chipman, D. M. Excited Electronic States of Small Water Clusters. J. Chem. Phys. 2005, 122, 044111. (8) Torii, H. Intermolecular Electron Density Modulations in Water and Their Effects on the Far-Infrared Spectral Profiles at 6 THz. J. Phys. Chem. B 2011, 115, 6636−6643. (9) Heyden, M.; Sun, J.; Forbert, H.; Mathias, G.; Havenith, M.; Marx, D. Understanding the Origins of Dipolar Couplings and Correlated Motion in the Vibrational Spectrum of Water. J. Phys. Chem. Lett. 2012, 3, 2135−2140. (10) Kremer, F.; Schönhals, A. Broadband Dielectric Spectroscopy; Springer: Berlin, Germany, 2003. (11) Schrö der, C.; Steinhauser, O. Using Fit Functions in Computational Dielectric Spectroscopy. J. Chem. Phys. 2010, 132, 244109−244109. (12) Barker, J.; Watts, R. Monte Carlo Studies of the Dielectric Properties of Water-Like Models. Mol. Phys. 1973, 26, 789−792. (13) Neumann, M. The Dielectric Constant of Water. Computer Simulations with the MCY Potential. J. Chem. Phys. 1985, 82, 5663. (14) Neumann, M. Dielectric Relaxation in Water. Computer Simulations with the TIP4P Potential. J. Chem. Phys. 1986, 85, 1567. (15) Alper, H. E.; Levy, R. M. Computer Simulations of the Dielectric Properties of Water: Studies of the Simple Point Charge and Transferrable Intermolecular Potential Models. J. Chem. Phys. 1989, 91, 1242. (16) Bopp, P. A.; Kornyshev, A. A.; Sutmann, G. Frequency and Wave-Vector Dependent Dielectric Function of Water: Collective Modes and Relaxation Spectra. J. Chem. Phys. 1998, 109, 1939−1958. (17) Rick, S. W.; Stuart, S. J.; Berne, B. Dynamical Fluctuating Charge Force Fields: Application to Liquid Water. J. Chem. Phys. 1994, 101, 6141. (18) van der Spoel, D.; van Maaren, P. J.; Berendsenb, H. J. A Systematic Study of Water Models for Molecular Simulation: Derivation of Water Models Optimized for Use with a Reaction Field. J. Chem. Phys. 1998, 108, 10220−10230. (19) Yu, H.; van Gunsteren, W. F. Charge-on-Spring Polarizable Water Models Revisited: From Water Clusters to Liquid Water to Ice. J. Chem. Phys. 2004, 121, 9549. (20) Levy, E.; Puzenko, A.; Kaatze, U.; Ishai, P. B.; Feldman, Y. Dielectric Spectra Broadening as the Signature of Dipole−Matrix Interaction. I. Water in Nonionic Solutions. J. Chem. Phys. 2012, 136, 114502. (21) Levy, E.; Puzenko, A.; Kaatze, U.; Ishai, P. B.; Feldman, Y. Dielectric Spectra Broadening as the Signature of Dipole−Matrix Interaction. II. Water in Ionic Solutions. J. Chem. Phys. 2012, 136, 114503. H
dx.doi.org/10.1021/jp507419e | J. Phys. Chem. A XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry A
Article
(47) Press, W.; Flannery, B. P.; Teukolsky, S.; Vetterling, W. Numerical Recipes; Cambridge Univ Press: Cambridge, U.K., 1992; Vol. 1, p 989. (48) Chen, W.; Sharma, M.; Resta, R.; Galli, G.; Car, R. Role of Dipolar Correlations in the Infrared Spectra of Water and Ice. Phys. Rev. B 2008, 77, 245114. (49) Borysow, J.; Moraldi, M.; Frommhold, L. The Collision Induced Spectroscopies: Concerning the Desymmetrization of Classical Line Shape. Mol. Phys. 1985, 56, 913−922. (50) Ramírez, R.; López-Ciudad, T.; Kumar, P. P.; Marx, D. Quantum Corrections to Classical Time-Correlation Functions: Hydrogen Bonding and Anharmonic Floppy Modes. J. Chem. Phys. 2004, 121, 3973. (51) Wooten, F. Optical Properties of Solids; Academic Press: New York, 1972. (52) Bader, J. S.; Berne, B. Quantum and Classical Relaxation Rates from Classical Simulations. J. Chem. Phys. 1994, 100, 8359−8366. (53) Guillot, B. A Molecular Dynamics Study of the Far Infrared Spectrum of Liquid Water. J. Chem. Phys. 1991, 95, 1543−1551. (54) Brubach, J.-B.; Mermet, A.; Filabozzi, A.; Gerschel, A.; Roy, P. Signatures of the Hydrogen Bonding in the Infrared Bands of Water. J. Chem. Phys. 2005, 122, 184509. (55) Jedlovszky, P.; Brodholt, J.; Bruni, F.; Ricci, M.; Soper, A.; Vallauri, R. Analysis of the Hydrogen-Bonded Structure of Water from Ambient to Supercritical Conditions. J. Chem. Phys. 1998, 108, 8528− 8540. (56) Sega, M.; Kantorovich, S.; Holm, C.; Arnold, A. Communication: Kinetic and Pairing Contributions in the Dielectric Spectra of Electrolyte Solutions. J. Chem. Phys. 2014, 140, 211101. (57) Fernández, D. P.; Mulev, Y.; Goodwin, A.; Sengers, J. L. A Database for the Static Dielectric Constant of Water and Steam. J. Phys. Chem. Ref. Data 1995, 24, 33−70.
I
dx.doi.org/10.1021/jp507419e | J. Phys. Chem. A XXXX, XXX, XXX−XXX
