6192
J. Phys. Chem. B 1997, 101, 6192-6195
Molecular Dynamics Calculations for Ice Ih C. J. Burnham† and J.-C. Li* Department of Physics, UniVersity of Salford, Salford M5 4WT, U.K., and UMIST, PO Box 88, Manchester M60 1QD, U.K.
M. Leslie Daresbury Laboratory, Warrington WA4 4AD, U.K. ReceiVed: October 11, 1996; In Final Form: April 4, 1997X
Since the observation of two intermolecular translational optic peaks at 28 and 37 meV using inelastic neutron scattering techniques, there has been a renewed effort to reproduce phonon density of states with these peaks by computer simulation. In this paper, we demonstrate our attempt to calculate the measured inelastic neutron spectrum using molecular dynamic techniques for a number of potentials, such as TIP4P, MCY, and KKY potentials. The results show that none of them can reproduce the double peaks seen in the measured spectrum.
1. Introduction Li and Ross have proposed a “two strength hydrogen bond” model1 for ice in an attempt to explain the observed splitting at 28 and 37 meV in the phonon density of states (PDOS) spectra (Figure 3) of ice Ih and other phases in the translational region. Based on the hypothesis of two strengths of hydrogen bonding,1 lattice dynamics (LD) calculations1 have been carried out which appear to give good agreement with the experimentally measured results.2 This in itself is not proof that the two-strength model is the right one; it may be that other potentials in the literature can give equally good or better PDOS spectra. An important question at the present time to ask is whether any existing potential can also produce this splitting. In this paper molecular dynamics (MD) was used as the simulation technique in order to test this. The motive behind the two-strength model comes from the PDOS peak splitting at 28 and 37 meV. Although it is common to encounter what appears to be two such peaks in MD runs, this could result from poor statistics and insufficiently sized cells, which sample less phonons, producing a very spiky PDOS at frequencies corresponding to these peaks. Any true features of the PDOS will have to be invariant to the cell size, number of steps, and step size of the q-point in the first Brillouin zone (BZ). It was found that reasonable sampling of the BZ is in the region of 500 molecules,3 with typically 8000 steps of 10-15 s. Bearing these factors in mind, we have demonstrated that none of the more popular ice potentials are able to reproduce the peak splitting using MD. 2. Technical Details and Simulation Techniques Molecular dynamics as a simulation technique has definite drawbacks in comparison with lattice dynamics when calculating PDOS spectra. It is generally orders of magnitude slower than LD and needs a much larger unit cell in order to produce PDOS spectra of equal quality. Nevertheless, MD is a useful tool, which goes beyond the harmonic approximation of LD and can also be used to obtain structural information. The MD program used in the paper is called “MDCSPC4”.4 MDCSPC4 uses an * Corresponding author. † University of Salford. X Abstract published in AdVance ACS Abstracts, June 15, 1997.
S1089-5647(96)03259-2 CCC: $14.00
Ewald sum5 to carry out the force and energy evaluations for the charges. The motion of the molecules/atoms is found using a fifth-order Gear corrector algorithm to integrate the equation of motion. The proton disordered structures were created by using a “random walk” program which moves pair “ionic” defects around a geometrical representation of the simulation cell for a sufficient number of steps to produce random orientations of the water molecules consistent with the Bernal-Fowler ice rules. A PDOS can be calculated by the Fourier transform of the velocity autocorrelation function (VACF) using as input velocities from an MD run
g(ω) ) F {〈v(t)‚v(0)〉i,t} A nearly equivalent “fast Fourier transform” method was used, conferring an enormous advantage in speed.6 This method calculates the autocorrelation in the frequency domain using the relationship
vacf(t) ) 〈vi(t)‚vi(0)〉i,t ) 〈F -1 {|Vi( f)|2}〉i where Vi( f) ) F {Vi(t)}. By using the above relationship, a fast Fourier transform method can be used to find the VACF in the frequency domain. The PDOS spectra were created using a complete velocity file of the velocities of all the particles at every time step. A trajectory file was used to examine structural properties of the potentials tested. A good test of the stability of the crystal simulated is to examine the radial distribution function (RDF) to check the first, second, and third nearest neighbor distances. 3. Water-Water Potentials There are many water pair potentials in the literature.7 We have chosen a few typical ones that represent the rigid, flexible, and polarizable pair potential classes. They are described as follows. 3.1. Rigid Molecule Potentials. The TIP4P and MCY potentials7,8 were chosen as two commonly used examples of rigid molecule potentials. Both are four-site models, with two H sites holding a charge of +q/2 each, an “M” site (along the H-O-H bisector) with a charge of -q but no mass, and an O site with no charge. Both the O site and the H sites have mass. © 1997 American Chemical Society
Molecular Dynamics Calculations for Ice Ih All of the four sites are fixed on the molecule so no bending or stretching frequencies will be found. Six coordinates were used to specify each molecule for MDCSPC4, three spatial coordinates and three Euler angles. The H-O-H bond angle was set to 104.52° and the O-H length at 0.9572 Å. Both potentials were modeled with a 480-molecule simulation cell. These simulations, as is the case for all the simulations presented here, were carried out at a temperature of 10 K. The TIP4P potential uses a 6-12 interaction between the oxygen sites of molecules.
J. Phys. Chem. B, Vol. 101, No. 32, 1997 6193 3.2. The Kumagai-Kawamura-Yokokawa (KKY) Potential.9 This potential is a relatively recent addition to the literature, and some interesting results have been found using it.10,11 It differs from most of the existing potentials in that it contains both intermolecular and intramolecular motion; that is, it is not a rigid molecule potential. This gives extra degrees of freedom to the motion of the molecule; it can exhibit bending and stretching modes. We refer the reader to the original paper.9 Only the three-body term is given below.
UHOH(θHOH) ) -fk{cos[2(θHOH - θ0)] - 1}k1k2
A C V (rOO) ) 12 - 6 + electrostatic interactions rOO rOO (the parameters can be found in ref 7) It was found when running the MD simulation with a constant volume and temperature that the system became structurally unstable. Even though the first nearest neighbor distance was preserved at about 2.9 Å, the second and higher nearest neighbor peaks were washed out in the radial distribution function. This behavior would seem reasonable; a simple 6-12 potential is only capable of fixing the first nearest neighbor distance and has no preference for a tetrahedrally bonded structure. However, it was found that by keeping the volume of the cell fixed, the simulation was able to remain stable. This is not surprising; if the volume and rOO are fixed, the system will find it very hard to deviate from its initial configuration. Hence the input gives a tetrahedrally bonded structure and the molecules stay bonded in this way. Nearest neighbor molecules will move around in the energy minimum created by the pair potential and pairwise additive electrostatic forces. At low temperatures these nearest neighbors will sample only the parabolic part of the minimum, and anharmonic effects should play a very small role. Different potentials have different curvatures about their minima and thus different frequency positions, but it would be difficult to see how any modification to a simple pair potential could produce the qualitative difference we are interested in of two sharp peaks in the translational optic region. The MCY potential is of a slightly more complex form than the TIP4P:
VOO ) 1864.27 exp(-2.75311ROO) VOH ) 2.68445 exp(-1.43979ROH) 0.675342 exp(-1.14149ROH) VHH ) 0.662712 exp(-1.29998RHH) with
q ) 0.7175e on each H atom q ) -1.435e on H-O-H bisector 0.486 982 bohr from the O atom (Energies are given in Hartrees and distances are given in bohrs.) It was found by doing a simple energy calculation of a water dimer that the minimum in energy was around rOO ) 3.12 Å for the O-O separation. In fact, running the MD for the MCY potential using an initial O-O distance of 2.74 Å and a constant volume fails to give a stable structure. It was found that the simulation cell needs to be scaled up so that the O-O separation equaled 2.96 Å in order to stabilize the structure. This agrees fairly well with the result of 2.9663 Å for the O-O separation found by Morse and Rice for water.8
ki )
1 exp[gr(rOH(i) - rm)] + 1
where
θ0 ) 99.5°; gr ) 7.0 Å-1; rm ) 1.40 Å; fk ) 1.1 × 10-11 J It is helpful to go into some detail regarding the intramolecular motion if the potential is to be understood. The term is meant to give the bending modes from oscillation of the HOH angle. The -cos[2(θHOH - θ0)] factor has a minimum at θHOH ) θ0, so θ0 fixes the HOH mean angle. This angle is actually a few degrees higher than θ0 in practice; because of the HH repulsion and the hydrogen bonding, fk can be thought of as a force constant for the bending; the higher fk is, the higher the resultant frequencies. The k1k2 factor creates a step function, the purpose of which is to create two strengths of bonding at the covalent (r < 1 Å) and the hydrogen bond (2 Å > rOH > 1.5 Å) distances. It was found to be extremely important to the stability of the structure to include H-O-H and H-O-H three-body interactions in the cutoffs used. We were unable to produce a stable structure in the simulation for the ice Ih using the potential as quoted. It was suspected that the fk factor in the three-body part was a few orders of magnitude too high. To check this, a lattice dynamics program, PHONON12 was used. The major alteration to the program needed was the inclusion of the three-body term into the FORTRAN code. A simulation was carried out for a cell of 512 molecules. The PDOS was obtained from a histogram of the phonon frequencies found. We tried fitting the value of fk to the bending frequency peak found in the PDOS and arrived at a value of 1.2 × 10-19 J. The MD simulation using this potential used 512 molecules and was run with 16384 (×4 fs) steps. It was found that the simulation, unlike the other potentials in this paper, had to be run with a constant pressure option using Brown Clarke algorithm. 3.3. The Sprik-Klein (SK) Potential.13 The SK potential is one of a group of polarizable potentials; that is, the polarization of the individual molecules can change according to their relative environment. This is certainly worthy of our investigation, as the two-strength hydrogen bond model assigns the strengths according to the relative orientation of the water molecules, but one would expect a difference in bonding anyway due to the dipole-dipole interactions found in a polarizable model. To some extent the KKY potential of the last section behaves as a polarizable potential. The HOH bond angles and OH bond lengths can change, which will thus affect the dipole moment of the molecule. Sprik and Klein argue that there are two main contributions to the dipole moment of a water molecule, the gas phase dipole moment and the electronic polarization. The gas phase dipole moment is fixed by the OH distance and the HOH angle, giving a dipole moment of 1.85 D as measured experimentally. The remaining dipole moment
6194 J. Phys. Chem. B, Vol. 101, No. 32, 1997
Burnham et al.
TABLE 1: Frequency Positions for Different Potentialsa first trans. peak last trans. peak libration HOH-bending stretching a
KKY (Itoh et al.)
KKY (Burnham)
TIP4P
MCY
12.38 43.9 97-164.6 210.0 372.0-379.5
9.5 45.5 103.0-160.0 216.0 360.0-370.0
7.1 33.0 68.0-118.0
8.0 32.0 52.0-134
SK 6.95 37.0 63.0
expt 7.05 37.87 68.0-124.5 199.5
396.0-413.0
All energies are given in meV.
of around 1 D arises from the electronic polarization. Since the electronic distribution can alter on a time scale orders of magnitude smaller than the molecules and since this distribution is quantum mechanical in nature, the electrons cannot be simulated by standard MD techniques. Instead, the electronic polarization is allowed to change adiabatically; that is, the motion is ignored and a self-consistent charge distribution at each time step is found. Finding a perfectly self-consistent charge distribution is a major problem in itself. In some models (e.g. ref 14) this is brought about by an iterative process. In order to save on computational time and complexity, however, the SK model uses a fictitious mass equivalent to 1/2mh corresponding to the speed at which the dipole moment can change. At each time step the motion of the molecules and the change of the dipole moments are integrated. However, since the fictitious mass is the smallest mass in the simulation, the dipole moments change on a much quicker time scale than the molecules, thus approximating to the correct self-consistent polarization. The SK potential is based on the TIP4P rigid molecule potential dealt with in section 2. The polarization is represented by four sites containing charges displaced along tetrahedral directions from the M site (0.15 Å along the molecule bisector from the oxygen toward the hydrogen sites). The magnitude of the four charges on the sites is allowed to change with the condition that the molecule still has zero net charge, but the positions of the charges with respect to the molecule are fixed. The “force” that acts on the fictitious mass giving the charge “acceleration” is such that it acts to minimize the electrostatic and dipole energy of the molecules. A truly self-consistent dipole distribution would lie at the minimum of this energy every time step. In the limit of no intermolecular interaction they should have the gas phase dipole moment of 1.85 D. Since the OM separation is fixed to give this moment, there should be no additional electronic contribution. The only way this can occur is if the four tetrahedral sites all contain equal charge and intramolecular electrostatic interactions are not included. The four tetrahedrally positioned charges are modeled by Gaussians. The argument as given in Sprik and Klein’s paper is that dipole effects are only felt at long range, ‘long’ being in this case further than the OH distance. Closer to the molecule a test charge should not experience the four distinct charge sites but instead a smeared out distribution. 4. Results The PDOS spectra for the rigid molecule potentials were calculated from the center of mass VACF. For the TIP4P potential (Figure 1, Table 1) the first peak was found at 10.0 meV as compared to the experimental result of 7.1 meV. The last peak in the translation region is situated at around 33 meV (experiment gives 37 meV). There is very little intensity in the librational region because the calculated PDOS spectra for all the rigid molecule potentials were taken from the velocity data for the motion of the centers of mass. Figure 1 shows a librational band in the region 68-118 meV. This is in very close agreement to the experimental result of 68-124.5 meV.
Figure 1. Plot of g(ω) against energy for MD simulations of ice Ih using the TIP4P (upper curve) and MCY (lower curve) potentials.
The radial distribution function gives rOO as 2.75 Å, a result agreeing with the reported liquid rOO distance.6 The other rigid molecule potential reviewed here is the MCY potential (Figure 1, Table 1). The major difference is the first peak at 8 meV. The translational region ends at 32 meV, not very much different from the TIP4P result. The librational region is contained in the 52-134 meV range, which is higher than experimental results. The mean OO separation is at 2.95 Å. The PDOS obtained from the KKY potential (Figure 2, Table 1) is significantly different from the results from the rigid molecule potentials. Trivially, as mentioned earlier, the spectrum contains contributions arising from intramolecular bending and stretching peaks due to the fact that the oxygen and hydrogens are free to move with respect to each other within a molecule. However the translational frequencies are very different from those found in the rigid molecule simulations. The last peak is at 45.5 meV, much higher than produced by the TIP4P and the MCY potentials. The first peak however in the translational region is around 9.5 meV, comparable with the TIP4P potential. The librational region is between 103 and 169 meV, much higher than for the rigid molecule potentials. A possible reason for this is the inclusion of a hydrogen-bonded term. There is a peak due to bending at 216 meV. It must be borne in mind, however, that this is a result that derives from our fixing of the fk force constant (see section 5). The intramolecular stretching modes are found at 360 and 370 meV, compared with 396 and 413 meV for experiment.
Molecular Dynamics Calculations for Ice Ih
Figure 2. Plot of g(ω) against energy for MD simulations of ice Ih using the KKY (upper curve) and SK (lower curve) potentials.
J. Phys. Chem. B, Vol. 101, No. 32, 1997 6195 frequencies for the translational and librational regions and gives an rOO distance correct to 0.1 Å. Although the KKY potential contains extra degrees of freedom and can give internal bending and stretching modes, it does not appear to give especially accurate frequencies, though it seems the best potential in terms of producing the right structure for ice Ih. As mentioned earlier, the other potentials needed to be run using the constant volume option. The frequency positions found for the KKY potential differ somewhat from the results obtained by Itoh et al.10 This may be explained by our different values for fk, and until this question is resolved comparison between results is difficult. None of the potentials examined here have been found to reproduce the two sharp peaks found in experimental measurements taken by inelastic neutron scattering (see Figure 3). This does not mean that there aren’t potentials in the literature that can reproduce the measured PDOS. All we can say is that Li and Ross’s model at the moment produces a feature of the PDOS that we haven’t been able to produce from any other potentials we or anyone else to our knowledge have tested. Li and Ross’s model has a force constant dependent on the relative orientation of the water molecules. Dipole interactions also have an orientational interaction dependency, and it remains to be seen to what extent the dipoles affect the force constants. After completing the investigation of the SK potential, a next step is to see how Li and Ross’s ideas can be incorporated into a MD simulation. At the present, a simple harmonic spring approach has been used in LD simulations. This precludes the possibility of discovering anharmonic effects, but the harmonic approximation is usually a good one. It will be surprising if significant differences are found in MD simulations from the LD simulation. Also other phases of ice will need to be simulated, though they have not been experimentally investigated as fully as ice Ih. Acknowledgment. C.J.B. thanks his triumvirate of M.L. and J.C.L. (supervisors) and I. Morrison (advisor). This work has been funded by EPSRC. Also thank you to Sarah Burnham and Jessica Burnham for their patience. References and Notes
Figure 3. Plot of G(ω) against energy. G(ω) is the “generalized vibrational density of states” taken from inelastic neutron scattering data.1,2 G(ω) ) e-2WH ∑j,q〈|e Hj (q)|2〉, where e Hj (q) is the displacement eigenvector for the jth H atom at wavevector q.
The SK potential (Figure 2, Table 1) is still in the parametrization stage, but it is giving interesting results which we hope to publish in full shortly. It was found that the average dipole moments of the molecules is heavily dependent on the nearest neighbor molecular distances. By adjusting the parameters to obtain the correct 2.75 Å O-O distance, a total mean molecular dipole moment of 3.0 D was found, which can be compared to Sprik and Klein’s value of 2.88 D13 for liquid water. There are encouraging frequency positions in the PDOS. Also, there appear to be features at 28 and 37 meV though they do not manifest themselves as sharp peaks. This is very much a preliminary result, and no strong conclusions should be inferred at this stage. 5. Conclusions Overall perhaps the most satisfactory result from the nonpolarizable potentials is the TIP4P model. It produces good
(1) Li, J. C.; Ross, D. K. Nature 1993, 365, 327. Li, J. C. J. Chem. Phys. 1996, 105, 15. (2) Li, C. J.; Londono, J. D.; Ross, D. K.; Finney, J. L.; Tomkinson, J.; Sherman, W. F. J. Chem. Phys. 1991, 94, 6770. (3) Li, J. C. J. Phys. Chem., in press. (4) Smith, W. MDCSPC4; Daresbury Laboratory: Daresbury, Warrington, England WA4 4AD, 1991. Smith, W. Program MDCSPC4: Molecular Dynamics Computer Simulation of Molecular Ionic Crystals, Technical Memorandum DL/SCI/TM84T; Daresbury Laboratory: Warrington, 1991. (5) For example see: Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, 1990, p 156. (6) Smith, W. Information Quarterly for MD & MC Simulations, Number 7, Daresbury Laboratory: Warrington, 1982. (7) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D. J. Chem. Phys. 1983, 79 (2), 926. (8) Morse, M. D.; Rice, S. A. J. Chem. Phys. 1982, 76, 653. (9) Kumagai, N.; Kawamura, K.; Yokowawa, T. Mol. Simul. 1994, 12, 177. (10) Itoh, H.; Kawamura, K.; Hondoh, T.; Mae, S. Physica B 1996, 219,220, 469. (11) Itoh, H. Proceedings of PHONONS95; 1996. (12) Leslie, M. PHONON; Daresbury Laboratory, private communication. (13) Sprik, M.; Klein, M. L. J. Chem. Phys. 1988, 89, 7536. Sprik, M. J. Phys. Chem. 1991, 95, 2283. (14) Niesar, U.; Corongiu, G.; Clementi, E.; Kneller, G. R.; Bhattacharya, D. K. J. Phys. Chem. 1990, 94, 7949. (15) Coulson, C. A.; Eisenberg, D. Proc. R. Soc. A 1966, 291, 445.
