20518
J. Phys. Chem. C 2010, 114, 20518–20522
Hydrogen-Bond Defect in the Structure of Ice Ih† Jirˇ´ı Kolafa* and Milan Oncˇa´k* Department of Physical Chemistry, Institute of Chemical Technology, Prague, Technicka´ 5, 166 28 Praha 6, Czech Republic ReceiVed: May 31, 2010; ReVised Manuscript ReceiVed: August 19, 2010
We study the previously reported 5 + 7 defect in the structure of hydrogen bonds in water ice Ih. This defect appears in simulations with TIP4P water models, and it has not been clear whether it is an artifact of the model. We found that the TIP4P/2005 and TIP4P/ice models predict the probability of the 5 + 7 defect to be ∼10-5 (per water molecule), whereas the NE6 (six-site Nada-Eerden) model predicts it to be less than 10-10. The DFT/BHandHLYP and MP2 quantum chemical calculations suggest that the defect is real and its probability is around 4 × 10-8 to 2 × 10-9. Introduction
Methods
The common hexagonal ice is a complex system even at zero temperature, exhibiting a proton disorder leading to a residual entropy.1,2 To explain dielectric relaxation, occurring at a microsecond scale at temperatures slightly below the melting point,3 Bjerrum introduced two types of orientational defects: In the D-defect, there are two hydrogens between neighboring water molecules whereas in the L-defect, a hydrogen is missing. The Bjerrum defects are topological defects of the hydrogen bond network, that is, a single defect cannot be removed by a local bond rearrangement, and a single defect can be detected by summing the donor-acceptor directions over any surface enclosing it. These defects can be also viewed as pseudoparticles of effective charge of (0.38e.4 They also interact with vacancies and interstitials.5 In addition to modern spectroscopic6,7 and first principle8,9 studies of ice and its defects, ice structure and properties have been studied by classical molecular dynamics. Several simple models of water have been recently carefully parametrized to describe the structure of condensed phases of water10 and particularly ice polymorphs.11,12 Model TIP4P/200510 was fitted to a range of bulk fluid and crystal properties. Its ice Ih is slightly less stable (melting point 252.1 K).10 Ice melting in the TIP4P/ ice version with a higher quadrupole moment13 is well reproduced (272.2 K).11 Ice Ih in the six-site model NE6 is slightly more stable (280-285) K,14 which even facilitates homogeneous nucleation from supercooled water.15 Three-site models do not describe well water ices.16 The Bjerrum defects are well described by the classical models. However, an additional type of defect was found in ice Ih modeled by TIP4P water.17 This defect is not topological, that is, it is caused by a local rearrangement of hydrogen bonds. Since six-membered rings around a hydrogen bond are replaced by 5- and 7-rings, this defect was nicknamed the 5 + 7 defect. It was later rediscovered.18 Immediately, a question appears whether this defect is real or just an artifact of TIP4P models. To answer this question, we performed a study combining classical molecular mechanics (MM) and molecular dynamics (MD) with quantum calculations.
Models. For this study, we chose three classical nonpolarizable rigid water models. Two of the models, TIP4P/200510 and TIP4P/ice,11 have the partial negative charge moved from oxygen to a site called M on the molecule axis. The third model by Nada and Eerden,12 denoted as NE6, has two more interaction sites bearing a small negative partial charge and corresponding to lone electron pairs (denoted by L). Simulation Details. For MD calculations, we used the Verlet method with the SHAKE algorithm to keep the molecules rigid. Auxiliary site M in the TIP4P model was treated as dependent (its position was calculated from H and O, and forces were distributed back to H and O). Similarly sites O and M in the NE6 model were treated as dependent. To do this, we redistributed the masses of atoms in the NE6 model to external interaction sites (4 g/mol to both hydrogens and the rest to the L-sites) while O and M were kept massless. This trick allows for a longer time step (we used 1.25 or 2 fs) and does not affect any static quantity (including defect probability); however, kinetic properties are inevitably affected. The time step was 1 fs for the TIP4P models. We performed all MD calculations using a periodic sample of ice Ih containing 360 water molecules. This number allows for a simulation box that is approximately cubic. For a perfectly tetrahedral environment of each water molecule, this gives a box of dimensions (22.14, 22.59, 23.48) Å. Temperature was kept by a Berendsen thermostat with time constant of 2 ps (ideal gas equivalent). The diagonal components of the pressure tensor were kept separately at atmospheric pressure by a Berendsen barostat (with coupling time equivalent to 100 ps for ideal gas). The simulation box thus changes shape, but remains rectangular. The electrostatic forces were calculated by the Ewald summation (for 360 molecules the separation parameter was R ) 0.3477/Å and the maximum k-space vector was κ ) 0.3646/Å). It is a common practice to truncate the Lennard-Jones interactions and to add a correction calculated under the assumption that the atoms are distributed uniformly beyond the cutoff distance.19 This assumption can be questioned for a crystal. We use the following short-range approximation of the full Lennard-Jones potential ULJ,20
†
Part of the “Mark A. Ratner Festschrift”. * E-mails: (J.K.)
[email protected]; (M.O.)
[email protected].
10.1021/jp1049815 2010 American Chemical Society Published on Web 09/01/2010
Hydrogen-Bond Defect in the Structure of Ice Ih
{
uLJ(r)
u(r) ) A(r - C2 0 2
2 2
)
J. Phys. Chem. C, Vol. 114, No. 48, 2010 20519
for r < C1 for C1 < r < C2 for C2 < r
where C2 is the cutoff (we used C2 ) 10.5 Å) and C1 ) 8.13 Å was calculated from C2 so that the forces are continuous. Since the short-range potential switches smoothly from the full value to zero, the interference errors are small. We tested our truncation scheme on “system a” doubled in each direction using several longer cutoffs. The maximum error per one water molecule for C2 g 10 Å is less than 6 × 10-5 eV, and the defect energy error is less than 0.0003 eV. The MM ground state geometries were obtained by simulated annealing; that is, a MD simulation at zero thermostat temperature. The optimum thermostat time constant was 0.13 ps, and the barostat time constant was 0.06 ps (assuming the bulk modulus of ice 17 GPa; the actual bulk moduli are 10, 11, and 13 GPa for the TIP4P/2005, TIP4P/ice, and NE6 models, respectively). MM and MD calculations were performed by the MACSIMUS package.20 Proton Disorder. Proton disorder of the sample was generated by simulated annealing with the same penalty assigned to any Bjerrum defect. The lattice system was simulated by the standard Metropolis MC algorithm with decreasing temperature until all Bjerrum defects recombined. It was suggested to use zero dipole moment of the ice crystal.10 We therefore paid attention to the cell dipole moment. From 10 000 configurations obtained by simulated annealing (of the 360 molecule crystal), we determined the mean quadratic dipole moment (square root of its fluctuation) as 75.5 D. It does not differ too much from the expected value of 81.6 D (calculated from ice dielectric constant 95 by the fluctuation formula21). The percentage of cells with zero dipole moment obtained by the simulated annealing algorithm was 2.7%. We increased the number of samples with zero dipole moment to 7 in 26 and did not find any apparent effect on the results. Quantum Chemical Calculations. To approach the 5 + 7 ice defect by the quantum chemical calculations, we have employed the DFT/BHandHLYP/6-31+g* method for optimization and MP2/6-31+g* approaches for single-point recalculations. For easier applicability of these methods, we calculated smaller water clusters, which were cut out from the periodic ice. The first such cluster was composed of 26 water molecules, corresponding to the second neighbor shell around the central pair. The central pair is connected by hydrogen bonds to six molecules (see Figure 1). These six molecules are the same if the defect relaxes, although the topology of bonds is different. There are 18 second neighbors (see the graphical abstract). These second neighbors (the environment) are kept fixed during optimization, and only the central pair and its six neighbors are allowed to relax. These clusters are further denoted as 8 + 18. Similarly, we may fix the environment of the third neighbors (there are 32 of them) and optimize the inside of 26 water molecules or even the environment of the fourth neighbors (there are 58 of them) and optimize the inside of 56 molecules. These clusters are named 26 + 32 and 58 + 56, respectively. Both the defect and crystal DFT calculations must share the same fixed outer shell coordinates because the coordinates based on MM are not accurate enough for a quantum calculation (e.g., a small deviation of the H-O distance causes large energy difference). We used the arithmetic average of one calculation with the defect and one with the perfect crystal. Another
Figure 1. The 5 + 7 defect (defect “a” from Table 2) viewed from the hexagonal (z) axis. Hydrogen bonds 4 f 3 (donor f acceptor) and 7 f 5 of the original lattice (left) are broken, and new hydrogen bonds 7 f 4 and 5 f 3 of the defect are created (right); acceptor and donor of the central pair 5-4 are interchanged.
possibility is to perform two sets of calculations, one with the ideal crystal coordinates and the other with the optimized configuration with a defect that could bracket the true energy; however, we experienced convergence problems and had to abandon this approach. To ease the convergence of the DFT calculations, the initial configurations and environment for the DFT/BHandHLYP calculations were obtained from configurations of the periodic ice sample of 360 molecules preoptimized by the SCC-DFTB-D method22,23 with the recommended set of O and H dispersion constants. Unlike for the classical MM, only box volume (not shape) was optimized in the DFTB calculations (in addition to all coordinates), that is, the whole box was scaled by the same factor until the energy minimum was found. The DFT results for the MM-optimized environment for the 8 + 18 model systems were systematically lower by ∼0.05 eV with respect to the DFTB-optimized environment. We consider the DFTB-optimized environment as a better one because the O-H bonds and H-O-H angles are allowed to relax and therefore guarantee faster convergence as the cluster size increases. Without further discussion, we note that the defect energies calculated by the DFTB method in periodic boundary conditions are about 0.05 eV lower than ones for the TIP4P/2005 force field (see Table 2). All the calculations were repeated for four different defects and compared with classical calculation for the same cluster and in periodic boundary conditions. The DFT and MP2 calculations were performed in the Gaussian program package,24 and the DFTB+ program25 was employed for SCCDFTB-D calculations. Results and Discussion Classical Calculations. For each model (TIP4P/2005, TIP4P/ ice, and NE6), we performed 26 MD simulations with different initial configurations, lasting (after a short initial relaxation) 4 ns each. The temperature was set to 253 K for the least stable TIP4P/2005 ice; otherwise, it was 273 K. The trajectories were recorded at 1 ps intervals. The defects are apparent on the time dependence of the potential energy (see Figure 2). (The total energy, although less noisy, is delayed by the correlation time of the used Berendsen thermostat, and therefore, it is less suitable.) Lifetimes of the defects were “measured” simply by clicking a suitably rescaled plot at the energy in the middle between crystal and defect energy levels. At the same time, a configuration of a supposed defect was selected and optimized; in several cases of energy peaks shorter than a few picoseconds, the optimized configu-
20520
J. Phys. Chem. C, Vol. 114, No. 48, 2010
Figure 2. Dependence of the potential energy on time, 1 ns sample of the TIP4P/ice model.
ration corresponded to the ideal crystal, and therefore, this event was discarded. The defect represents a local energy minimum, and its energy is defined as the difference from the crystal configuration. The results on the defect statistics and energetics are collected in Table 1. The 5 + 7 defect is almost twice as probable (8.6 × 10-6 per one water molecule) in the general TIP4P/2005 model as in its ice-tuned cousin TIP4P/ice (4.9 × 10-6), although in the former case, the temperature was lowered to the model melting point. In addition, the averaged defect energy is smaller. No defect of this type was found for the NE6 model; its probability is thus less than ∼10-10. During simulations, we also detected simpler defects with one hydrogen bond broken and one of both affected molecules displaced into a cage of the ice Ih structure. In TIP4P models, this defect is much less probable than defect 5 + 7, and its energy is higher. One instance of this defect appeared also in the NE6 model. We will not study this defect in detail and will reserve the word “defect” for the 5 + 7 defect. In the most common 5 + 7 defect, a pair of neighboring water molecules (denoted 4 and 5 in Figure 1) partially reconnect their hydrogen bonds and are displaced toward the neighboring cavities in the ice structure. Each molecule of the pair breaks one hydrogen bond to one of its lattice neighbors (4-3 and 7-5) and then cross-connects (7-4 and 5-3). The mutual bond (5-4) is broken and reconnected during the process. In some cases (as in Figure 1), but not always, donor and acceptor are interchanged. In all cases, the hydrogen bonding in the original lattice around the central bond of the defect has a staggered conformation. The defect with eclipsed conformation is in principle possible, but it has higher energy and, thus, lower probability. We detected one in a testing run with a flexible variant of the TIP4P/ice model. We have never observed proton disorder relaxation, that is, in all cases, the lattice relaxed to the same configuration of protons (however, we observed rotation of a water molecule around its symmetry axis). The dielectric relaxation occurs in real ice at microsecond scale,3 which is far from available computer time. Indeed, the proton configurations of the six first
Kolafa and Oncˇa´k neighbors are unchanged so that the defect must return to the same proton disorder state, and the same arguments apply to any defect unless it spans a whole six-membered ring. The considered defect may be nevertheless the first step to such a relaxation because it breaks two six-membered rings and replaces them with two five-membered rings. Another mechanism includes pairs of Bjerrum defects. The found probabilities of defect occurrence (8.6 × 10-6 and 4.9 × 10-6 for TIP4P/2005 and TIP4P/ice, respectively; see Table 1) are in approximate agreement with the Boltzmann probabilities calculated from the defect energies. The simplest calculations from the averaged energies give 3.5 × 10-6 and 4.3 × 10-6 for the TIP4P/2005 and TIP4P/ice models, respectively (note that there are two bonds per one water molecule). However, the defects with lower energy are more probable and vice versa. As the first approximation, we may assume that the obtained set of defects is a representative sample of all possible 5 + 7 defects (more precisely, a defect is created with a probability given by the energy barrier, and its living time is given by the barrier of the backward reaction). In addition, we assume that the distribution of defect energies is Gaussian and employ the MD simulation results for averaged energies and standard deviations (see Table 1). By integration of the Gaussian distribution, we then get slightly increased probabilities 6.9 × 10-6 and 6.7 × 10-6, respectively. These probabilities match the directly determined probabilities within an accuracy of 20-30%, which is an acceptable agreement taking into account the small number of events detected. Since the probability is given by the Gibbs energy and we used the potential minima (zero-temperature internal energies), this means that the vibrational contributions to entropies of both the crystal and defect states are approximately the same (with uncertainty R ln 1.3 ≈ 2 J K-1 mol-1). From the NE6 defect energy, we can estimate that the defect probability of this model is less than ∼10-10 per molecule. Quantum Chemical Calculations. The main question, of course, is what is the real energy of the defect and, consequently, its concentration in real ice. Energy minima of four selected defects as calculated using the MM approach, and the DFT/BHandHLYP and MP2 methods are collected in Table 2. We will use the classical results to estimate the convergence of quantum calculations to an infinite system ab initio value. The defect energies calculated by MM for both the periodic system and clusters are lowest for the TIP4P/2005 model, slightly higher for the TIP4P/ice model, and highest for the NE model. The DFT/BHandHLYP and MP2 values for the same cluster lie between the TIP4P/ice and NE results (with one exception), with the MP2 energies shifted down by ∼0.07 eV with respect to the DFT/BHandHLYP ones. For the MM results, a change in the relative defect energies with the size of the model cluster is highest when passing from cluster 8 + 18 to 26 + 32. For cluster 58 + 56, the values are essentially converged to the results of periodic calculations. Trends in the DFT/BHandHLYP/6-31+g* energies are less obvious, partly because we experienced convergence problems for two largest clusters. However, energies seem to be basically converged (within 0.05 eV) already for 26 + 32, with 58 + 56 arrangement providing certain improvement. The averaged difference DFT minus TIP4P/ice is 0.157 eV for the largest cluster 58 + 56 and defects “c” and “d” and 0.215 eV for cluster 26 + 32. Transformed into probability, it predicts the defect probability 1.3 × 10-3 to 1.1 × 10-4 times smaller than that for the TIP4P/ice; namely, 6 × 10-9 to 5 ×
Hydrogen-Bond Defect in the Structure of Ice Ih
J. Phys. Chem. C, Vol. 114, No. 48, 2010 20521
TABLE 1: Statistics of Defects Detected in 26 Ih Ice Samples of 360 Water Molecules in 4 ns Runsa 5 + 7 (two reconnected bonds) -6
model
T/K
n
p/10
TIP4P/2005 TIP4P/ice NE6
253 273 273
81 88 0
8.6 4.9 0
one broken bond
τmax/ps
E/eV
n
p/10
320 178 0
0.289 ( 0.036 0.307 ( 0.036
6 8 1
7 9 1
-8
τmax/ps
E/eV
10 11 3
0.424 ( 0.029 0.428 ( 0.016 0.45
a Numbers n denote the numbers of defects; p, the defect probability per one water molecule; τmax, the maximum detected lifetime; E is the averaged potential energy (weighted by the lifetime) with the standard deviation of the sample. Only the 5 + 7 defect is studied extensively in this work.
TABLE 2: Energies (in eV) of Four Selected Defects Calculated by Various Methodsa system
method
TIP4P/2005
TIP4P/ice
NE6
a
360 per 8 + 18 26 + 32 58 + 56 360 per 8 + 18 26 + 32 58 + 56 360 per 8 + 18 26 + 32 58 + 56 360 per 8 + 18 26 + 32 58 + 56
0.348 0.359 0.352 0.355 0.247 0.163 0.213 0.248 0.273 0.193 0.240 0.242 0.322 0.271 0.289 0.303
0.377 0.389 0.381 0.384 0.268 0.177 0.232 0.269 0.297 0.210 0.260 0.264 0.349 0.293 0.313 0.328
0.572 0.583 0.589 0.577 0.540 0.432 0.497 0.529 0.563 0.510 0.533 0.528 0.560 0.529 0.526 0.555
b
c
d
DFT
MP2
0.522 0.608
0.468
0.385 0.481
0.303
0.383 0.456 0.402
0.311
0.428 0.490 0.504
0.359
a Calculation with 360 molecules in a periodic box is denoted as “360 per”. Symbol 8 + 18 means that 8 inside molecules were optimized and outer 18 neighbors constrained; other abbreviations are defined analogously. TIP4P/2005, TIP4P/ice, and NE6 denote classical water models; DFT stands for DFT/BHandHLYP/6-31+g* (including optimization); and MP2, for MP2/6-31+g* (single-point recalculation). Due to convergence problems, only “c” and “d” defect energies are given for the largest 58 + 56 cluster for the DFT approach.
10-10 (geometric average 1.7 × 10-9). Direct calculations based on DFT energies give approximately the same range, or slightly more if the distribution of energies is included. The systematically lower MP2 values give 20 times higher probabilities. Taking all these hints into account, we estimate the defect probability around 2 × 10-9 and 4 × 10-8 for the DFT and MP2 approaches, respectively, suggesting the degree of uncertainty of the ab initio calculations. Reaction Coordinate. Our observation of several defects during classical molecular dynamics trajectories suggests that there are at least two reaction paths connecting crystal and defect geometries. They differ by the rotation of one of the central water molecules (identity of both hydrogens). We analyzed in detail one of the paths of defect “a” from Table 2. The results are based on MM of the full periodic box with N ) 360 TIP4P/ 2005 molecules. The x, y, and z box dimensions were optimized separately. Figure 3 shows the energy along a suggested reaction pathway, and Figure 4, six selected distances of hydrogenoxygen pairs that change their bonding state: three bonds are broken and three new established. The x axis is the massweighted distance integrated along the reaction coordinate curve in the 3N-dimensional (N ) 360) Cartesian space (more precisely, since the energy minimization included separate relaxation of all three box sides, both configurations were rescaled to the averaged box before the mass-weighted Cartesian metric was calculated). For compatibility, the reaction coordinate
Figure 3. Energy along the reaction coordinate of “defect a” as modeled by TIP4P/2005 water and recalculated by DFT/BHandHLYP/ 6-31+g* (26 + 32 cluster). The single-point recalculations of DFT geometries at the MP2/6-31g* level are also supplied.
is given in atomic units, 1 au ) me1/2a0, where me is the electron mass and a0 is the Bohr radius. It is seen that the exchange of hydrogen bonds proceeds gradually in several steps. In the transition state, two hydrogens of the central atoms 4 and 5 are close together (1.95 Å). It was not easy to find the transition state. The usual initial approximation of the reaction pathway by a linear interpolation between the respective local minima lies too far from the transition state. We first constrained 5 f 3 and 7 f 4 bonds (see Figure 4) and made 2D scans starting from both local minima. However, we obtained two crossing, but not connected, sheets. We therefore also scanned the third bond, 4 f 3, with a resolution of 0.05 Å in the region of the intersection of the sheets and thus located a 3D saddle-point, the transition state. The reaction coordinate was then reconstructed and verified by simulated annealing with very high friction coefficient from around the transition state; another small maximum at reaction coordinate of 1700 au (secondary transition state) had to be overcome. The scan errors in positions are responsible for small kinks seen on the curves. The TIP4P/2005 activation energy to create “defect a” is 0.524 eV, to relax back 0.176 eV, and the energy difference is 0.348 eV. Since the typical defect energies are smaller, the respective activation energies likely are also smaller. (To get the observed rate of defect creation of 0.8/ns for TIP4P/2005, provided that the typical correlation time is 0.1 ps and there are 720 affected bonds in the system, we need the activation
20522
J. Phys. Chem. C, Vol. 114, No. 48, 2010
Kolafa and Oncˇa´k was supported by the quantum chemical calculations. The defect may play a role in the dielectric relaxation of ice. The defect probability as retrieved from the quantum chemical calculations, around 4 × 10-8 to 2 × 10-9, is overestimated by TIP4P-like water models and underestimated by the NE6 model. This observation is in accordance with knowledge that TIP4P models are “less tetrahedral” than real water (and give less stable ice Ih unless the hydrogen bond strength is increased, as in TIP4P/ice), but NE6 is “too tetrahedral”. The defect energetics may thus serve as a test of a force field quality. Acknowledgment. This work was supported by the Czech Science Foundation (project P208/10/1724) and computing facilities by the Czech Ministry of Education (Center for Biomolecules and Complex Molecular Systems, project LC512). References and Notes
Figure 4. Selected distances along the reaction coordinate.
energy as low as 0.35 eV; it may be slightly higher if we take into account entropic effects). The classical reaction coordinate served as a basis for DFT calculations on the 26 + 32 cluster. The initial coordinates of 26 inside molecules were taken from the classical (TIP4P/2005) calculation; the environment (outer shell of 32 molecules) was the average of defect and cluster (at the TIP4P/2005 level, i.e., without the DFTB step). The same constraints were applied (in addition to freezing the environment): 5 f 3 distance (H of molecule 5, O of molecule 3) in the path from crystal to the transition state, three distances (of 5 f 3, 7 f 4, and 4 f 3) at the classical transition state, and 7 f 4 in the path from the transition state to defect. No bonds were constrained at the ends of the path (crystal minimum and defect local minimum). The results were plotted at the reaction coordinate values corresponding to the underlying classical configurations. It is seen that the activation energies calculated by the DFT approach are higher than the corresponding TIP4P/2005 energies, with values of 0.916 and 0.333 eV for the crystal to defect and back, respectively (defect energy is 0.583 eV). The singlepoint recalculations of the DFT geometries at the MP2/6-31g* level (note that the diffuse functions are omitted from the basis set) give the respective activation energies as 0.855 and 0.359 eV (defect energy, 0.496 eV). This is in accord with the energetics discussion in the previous section. Conclusions We found that the 5 + 7 defect in the structure of hydrogen bonds in water ice Ih is most probably real, that is, its existence
(1) Pauling, L. J. Am. Chem. Soc. 1935, 57, 2680–2684. (2) Nagle, J. F. J. Math. Phys. 1966, 7, 1484–1496. (3) Johari, G. P.; Whaley, E. J. Chem. Phys. 1981, 75, 1333–1339. (4) de Koning, M.; Antonelli, A.; da Silva, A. J. R.; Fazzio, A. Phys. ReV. Lett. 2006, 96, 075501. (5) de Koning, M.; Antonelli, A. J. Phys. Chem. B 2007, 111, 12537– 12542. (6) Devlin, J. P.; Buch, V. J. Chem. Phys. 2007, 127, 091101. (7) Tse, J. S.; Shaw, D. M.; Klug, D. D.; Patchkovskii, S.; Vanko, G.; Monaco, G.; Krisch, M. Phys. ReV. Lett. 2008, 100, 095502. (8) de Koning, M.; Antonelli, A. J. Chem. Phys. 2008, 128, 164502. (9) de Koning, M.; Antonelli, A.; da Silva, A. J. R.; Fazzio, A. Phys. ReV. Lett. 2006, 97, 155501. (10) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (11) Abascal, J. L. F.; Sanz, E.; Fernandez, R. G.; Vega, C. J. Chem. Phys. 2005, 122, 234511. (12) Nada, H.; van der Eerden, J. P. J. M. J. Chem. Phys. 2003, 118, 7401–7413. (13) Abascal, J. L. F.; Vega, C. J. Phys. Chem. C 2007, 111, 15811– 15822. (14) Nada, H.; Furukawab, Y. J. Cryst. Growth 2005, 283, 242–256. (15) Pluharˇova´, E.; Vrbka, L.; Jungwirth, P. J. Phys. Chem. C 2010, 114, 7831–7838. (16) Vega, C.; McBride, C.; Sanz, E.; Abascal, J. L. F. Phys. Chem. Chem. Phys. 2005, 7, 1450–1456. (17) Grishina, N.; Buch, V. J. Chem. Phys. 2004, 120, 5217–5225. (18) Lindberg, G. E.; Wang, F. J. Phys. Chem. B 2008, 112, 6436– 6441. (19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1986. (20) http://www.vscht.cz/fch/software/macsimus; last accessed on August 31, 2010. (21) Perram, J. W.; Smith, E. R. J. Stat. Phys. 1987, 46, 179–190. (22) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. ReV. B 1998, 58, 7260–7268. (23) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. J. Chem. Phys. 2001, 114, 5149–5155. (24) Frisch, M. J. et al. Gaussian 03, Revision E.01”, Gaussian, Inc.: Wallingford, CT, 2004. (25) Aradi, B.; Hourahine, B.; Frauenheim, T. J. Phys. Chem. A 2007, 111, 5678–5684.
JP1049815
