Interpretation of inelastic neutron scattering spectra for water ice by lattice and molecular dynamic simulations

P.B. Balbuena and J.M. Seminario (Editors) Molecular Dynamics. From Classical to Quantum Methods Theoretical and Computational Chemistry, Vol. 7 9 Elsevier Science B.V. All rights reserved

471

Chapter 12

Interpretation of inelastic neutron scattering spectra for water ice by lattice and molecular dynamic simulations Jichen Li a and John Tomkinson b aDepartment of Physics, University of Manchester Institute of Science and Technology (UMIST), PO Box 88, Manchester, M60 1QD, UK bISIS Facility, Rutherford Appleton Laboratory, Chilton, Didcot, Oxon, OX11 0QX, UK

In this article we describe a wide range of simulations of the inelastic incoherent neutron scattering spectra of ice Ih. These simulations use a variety of different water-water potentials from simple pair-wise (rigid and non-rigid molecule), to more sophisticated polarisable potentials. We demonstrate that in order to reproduce the measured neutron spectrum, greater anisotropy (or orientational variation) is required than these potentials presently provide.

1. INTRODUCTION Hydrogen bonding is one of the most important and intriguing molecular interactions. It is the physical basis for proton conductivity, it imposes the biologically important tertiary-structure necessary for life and is responsible for the unique properties of water. Scientists, for many years and across several different disciplines, have endeavoured to understand the complex nature of water and other hydrogen bonded systems. Despite considerable scientific effort there is still no coherent explanation for most of the properties of water, often referred to as its 'anomalies' [1-3]. Water has a large bond energy and an asymmetric hydrogen bond geometry. Moreover, the electrons in its oxygen (2sZp3) orbitals can easily rehybridise in response to the relative orientations of adjacent molecules. These properties give water, and ice, a number of abnormal properties, which cannot be explained by the 'ordinary rules' of physics and

472

chemistry. As a consequence, a large number of models have been proposed in attempts to interpret some of these properties of water (such as water's high heat capacity, high melting and boiling temperatures, and it's density and entropy fluctuations [2,3,8]). Meanwhile, a significant number of water potentials have been proposed in order to reproduce these properties by computer simulations. Some of these are based on ab-initio quantum mechanical calculations for the water dimer (e.g. MCY [4] and ST2 [5]), others are intuitive (e.g. TIP4P [6] and KKY [7]). Some succeed in reproducing the structure of water while others work better in reproducing its thermodynamic properties. However, no potential can yet provide a coherent explanation of or a complete model for the anomalies of water. Gradually it has been realised that pair-wise potentials are insufficient for water, partly because many-body interactions play an important role in organising the electronic obitals around the hydrogen bonds. Unfortunately these are inadequately taken into account in the classic pair-wise type of potentials. Hence a new class of water potentials has emerged, the polarisable potentials [8], which are much more successful in many ways than the pair-wise additive ones (for details see section 3). The quantitative studies of the properties of water and ice require detailed consideration of the forces acting on the atoms and the molecules. Experimental information about the strength of the hydrogen bond interaction can be obtained directly by measuring vibrational spectra. A particular vibrational mode (or phonon) is determined by the interatomic force constants, which in turn are the double differentials of the potential function. Therefore, measuring dynamic properties constitutes one of the most powerful ways of investigating interatomic potentials in a given material. Such investigations are traditionally carried out by means of optical spectroscopy, such as IR absorption and Raman scattering. These are very powerful techniques, which have been highly refined, and their usage has resulted in extensive and valuable data for water [ 10,11 ] and ice [ 12-15]. In water and ice, however, the normal selection rules governing the interaction of radiation with matter are broken due to the local structural disorder (or proton disorder in crystalline structures of ice), and analysis of the spectral intensities is difficult in general. On the other hand, although IR and Raman spectra are very sensitive to the intramolecular modes involving the O-H stretching and bending, they are less sensitive to the intermolecular modes involving the vibrations of water molecules against each other. Therefore, under normal circumstances, optical spectroscopy provides only limited data in the translational region which is vitally important for obtaining direct information about the hydrogen bond interaction. For instance, the acoustic and some of the optic frequencies have not been observed by the standard IR and Raman techniques; the spectra in this region show a predominate peak at 27 meV (or 220 cm 1) with an additional shoulder in the right hand side of the main peak at

473

37 meV (or 300 cm 1) for Raman spectrum as shown in Fig. 1. The difference between optic and inelastic neutron scattering (INS) spectra shown in the figure is still not fully understood. ,,

,

,

,,

,

. . . . . . . . . . . . . .

RAMAN

-

-

*

,

'

.--

I

i

>03 Z Z

b

--

..,.

0

,

I

1000

,

I

2000

i

,

I

3OO0

.

I

401

ENERGY TRANSFER (cm-1)i Fig. 1. Comparison of the spectra of ice Ih measured by IR, Raman [9] and INS [14,15] techniques show that the IbiS gives more detailed information on the translational (<400 cm "l) and librational modes (400-1200 cml). The curve (a) is the spectrum measured on TFXA and (b) on HET at ISIS (UK). In contrast, neutron spectroscopy is a more powerful probe, its results are directly proportional to the phonon density of states (DOS) (see Fig. 2) which can be vigorously calculated by lattice dynamics (LD) and molecular dynamics (MD). Applying these simulation techniques provide an excellent opportunity for constructing and testing potential functions. Because optical selection rules are not involved, INS measures all modes (IR/Raman measure the modes at the Brillouin Zone (BZ) q = 0, see Fig. 2) and is particularly suitable for studying disordered systems (or liquids). It hence provides direct information on the hydrogen bond interactions in water and ice.

474

1.1. Inelastic neutron scattering techniques The development of high fluxes of low-energy neutrons from pulsed accelerator sources, such as IPNS at Argonne National Laboratory (USA), KEK (Japan) and, especially, ISIS at Rutherford Appleton Laboratory (UK) provides a most suitable probe for studies of the vibrational dynamics of hydrogenous materials. The advantages of neutron scattering for the study of molecular dynamics stem from their remarkable properties.

Neutron Dispersion c u r v e s

Optic

=>

>

E >(.9 fig LLI

Z

LM

r

0

INTENSITY

r

0

qmax

Fig. 2. A schematic illustration of the difference of scattering intensities between the IR/Raman and neutron scattering techniques to the relationship of dispersion curves. For instance, Infrared spectroscopy measures frequencies at the BZ centre, q = 0, the peaks shown are relatively sharp, the width of the peaks is determined by the resolution of the instrument used. In an INS experiment, a broadened spectrum for each dispersion curve was observed; the spectrum has higher intensity at the fiat part of the dispersion curve at the BZ boundary. Hence the mode assignment is not appropriate for the INS spectrum. Thermal neutrons have energies comparable to the excitation energies of molecular solids and because of their mass they carry momentum. This momentum or wave-vector, k, is conventionally represented through the characteristic deBroglie wavelength, ~. and hence, k = 2n/X. Typical neutron wavelengths match the interatomic distances in solids, ca. 2 A and, unlike the

475 photon techniques, a single neutron scattering event is simultaneously sensitive to the structure and dynamics of the measured system. The molecular phenomena of relevance to the present study have energies of the order of 10~ to 103 meV (1 eV = 8050 cm 1) and wave-vectors about 10.3 - 10 &-l. Neutrons are, therefore, the ideal probe for revealing the full scope of the dynamics of ice. In an INS experiment it is the variation of scattering intensity with neutron energy transfer, h co, and momentum transfer, Q, that is observed. The energy and momentum transfer are related:

h.o= ~ - f : =(h 2 /2m).(k,2 -k: 2)

(I)

Q=(~ -k:)

(2)

where Ei~ is the initial (final) neutron energy and m is the neutron mass. The range of h co and Q in which measurements are performed characterise the INS scattering experiment. This can be compared with the other well-established experimental techniques of photon absorption (IR) and scattering (Raman). The ranges of energy transfer are identical but the wavelengths of light are much longer, by a factor of ca. 1000. Photon techniques probe much longer distances in solids and so are governed by bulk properties such as the electro-optic parameters of the sample. Neutrons are free from these restrictions and also the well-known optical selection rules because they penetrate the electron clouds around atoms and scatter from the nuclei. The total scattering observed from a system of atoms is given by the scattering length, b, and is composed of two components: first, incoherent scattering which is scattered almost isotropically from the sample, and; second, coherent scattering which is characterised by dramatic variations in intensity with scattering angle. This variation is nothing more than Bragg scattering and it occurs wherever a neutron wave-front is simultaneously scattered from several nuclei. The resulting, scattered, wave-fronts can undergo constructive, or destructive, interference and yield intense features in specific scattering directions. This fact is commonly exploited in elastic coherent neutron scattering measurements and it provides information on the equilibrium structure, this an extensive subject area and has been reviewed elsewhere [ 16]. The incoherent scattering is given by the difference between the squares of the average scattering length and the average of the squares of the scattering length.

476 Incoherence arises because samples are, usually, not isotopically pure or, alternatively, the atomic spins are not aligned. The latter is a very important consideration for all hydrogenous systems but especially for water or ice and it results in the dominant cross section for hydrogen (H) being incoherent. The incoherent signal contains no information on the position or motion of the scattering atom relative other atoms but rather refers to a single hydrogen nucleus. The magnitudes of the scattering cross sections of particular nuclei are a fact of experimental observation and there is presently no theoretical estimate of these values. For H atoms in normal water (H20) ~ n e n (= 80 barn, 1 barn = l xl 0 2s m 2) is much larger than trcohH (-- 2 barn), hence the coherent contribution to the scattering intensity is insignificant and can therefore be ignored. Neutrons interact with nuclei through the strong nuclear force, however, this force is so short ranged and neutrons interact (scatter) so infrequently that the results of the process are accurately described within a weak perturbation calculation, namely the Born model. Here we make no attempt to derive the relevant scattering functions since there is an extensive body of literature that addresses this subject directly [16]. Instead we shall exploit the situation and simply report the relevant results in a manner that allows a straightforward physical interpretation. The experimental observable is the rate at which neutrons are scattered into a detector, with a given solid angle, as a function of the energy exchanged with the sample. This is the double differential scattering cross-section:

q~z IV'(r)l g, q~, )126( E, - Ey + hco)

(4)

Where the neutron's wave-functions are ~ and the sample's wave-functions are qJ. The four terms in the expression are: first, the ratio of the incident and final neutron moments; second, the fundamental constants; third, the relationship between initial and final states, and finally; fourth, the condition of total energy conservation. The final term ensures that the difference between the incident and final neutron energies, El, equals a quantised energy state of the system, h0~, or is zero for elastic scattering. Only one functional form for V(r) successfully reproduces a spherical (S-wave) final neutron wave-function. This is the Fermi pseudo-potential, arising from a series of atoms, a, at positions R~ 2~h 2

V'(r)- ~ ~ - ' ~ b a 6 ( r - R a ) m

a

(51

477

S-wave scattering is the only practical outcome since P-wave final neutron states are not accessible to thermal neutrons, because these wave functions have negligible amplitude at the small radial values that are typical of atomic nuclei. It is convenient to rewrite the equation as a dynamical structure factor (or Scattering Law), which emphasises the dynamics of the sample.

(6) Here we are specifically ignoring all but the incoherent contributions and oo

S(Q, co) - ~ _1 ~

idt. exp(-icot)~(exp(-iQ. R~ (0)). exp(iQ-R a (t)) a

(7)

This refers the position, Ra(t), of the target H atom a at time, t, if it was originally at Ra(0) at t = 0. We choose the conventional decomposition of the position vector in terms of displacements away from its equilibrium position:

(8)

R ( t ) = R Equ~.b~u~ + U rr~., + U Rot + U Wb

The total dynamical structure factor is thus separated into three contributions (9)

STota I -- STranslati m i~) SRotatio n I~) Svibration

In the case of molecular solids the translational component, for instance, refers to displacements of the scattering atom by virtue of whole body translational vibrations of the molecule and provides oO

S(Q,CO)T= Idt.exp(-ico.t) exp(-iQ2.(UT(O)2)) exp(iQ2.(UT(O).Ur(t))) DO0

9

time independent

9

(10)

time dependent

where we have used the identity (exp A.exp B ) - exp(A2 ) 9exp(A.B). In performing the Fourier transform from the time to the frequency domain, it can be appreciated that the time independent terms will include contributions from all frequencies whilst time dependent terms will be specific to particular frequencies. The expression for the exchange of n quanta with an isotropic harmonic system is:

478 S(Q, nw) = exp(n| where |

(11)

e x p ( - ~ ) . I, (A)

h~.,

O=Q2

2keT

h coth(|

A=Q2~h

2,uco

2,uco

cosech(|

The (9 term is a Boltzmann temperature, T, factor; the 9 term is a Debye-Waller (DW) factor, well known from diffraction work; and finally, the I, are Bessel functions of the first kind. This single, simple, harmonic system is quite unrepresentative of any realistic molecular solid, however, it is instructive to express the equation with appropriate experimental parameters. The translational optic modes of water (see below) appear at about 35 meV (260 cm l) and typical experimental temperatures are about 20 K (= 2 meV). 19 = ~

hco

2kBT

35 ~ -- ~ 9

4

.'. coth(|

= 1.0

and

A<
and

cosech(|

<< 1.0 (12)

=Q2

h 2/2o9

Since the argument of the Bessel function is small, we expand it according to 1

13,

Substituting, expanding and simplifying yields S(Q, nco) oc (Q2Uz)" exp(-QZU z ) n[

(14)

This remarkably straightforward expression is the physical basis for the interpretation of INS spectra by lattice dynamic and molecular dynamic approaches. The Mean Square Displacement (MSD) of the scattering atom, U 2, is a function of the energy, co, and mass, ~, of the oscillator: [U2[-

h 2ktco

or

U2(A2) -

2.0717 ~(amu) 9co(meV)

(15)

479

We see that the observed INS intensity of a vibrational transition is directly related to the MSD of the scattering atom in the mode of interest and the momentum transferred to the system, Q2, during the scattering event. The MSD is a function only of the sample, i.e. of the forces acting on the atom and Q2 is a function only of the spectrometer, i.e. how the experimem was performed. Returning to our consideration of real molecular solids. The dynamical structure factor must encapsulate the richness found in the INS spectrum obtained from any realistic sample. This, fortunately, is a complexity only of degree. It arises from three effects; first, the large number of vibrational states available to a system; second, the few, if any, truly isotropic atomic displacements; and third, phonon dispersion. Anisotropy is to be anticipated in any system where strong, highly directional, bonds dominate the local force field. The vectors in the response function can be rewritten using dyadic notation [18], where qo is the angle between the momentum transfer vector and the atomic displacement vector. Below, the displacemem, B, is described as a weighted unit tensor, e 2. Where the reciprocal lattice is referred to the Cartesian system (i,j,k), and 9~2 is the real space (x,y,z) unit tensor. (Q.U)2 _ (QU. cos(~p))2 = ( Q Q r . UU r) (16) u u r - u 2 .e2(i,j,k)

- B

; QQr =Q2.~2(x,y,z)

Only the components of displacement along the scattering vector are effective in producing an INS response. If the vectors are parallel the maximum response is obtained but this falls to zero, as the vectors become orthogonal. Each of the modes, i, displaces the scattering atom in a characteristic direction as described by the tensor Bi. The total single quantum response is the sum over all atoms, a, in all modes.

S(Q, co, ) = E ~ QQ" B ~. exp(-QQ" A'* ) a

(17)

i

A o = ~ B ai i

Phonon dispersion expresses the phase relationship between different molecules in the lattice as they participate in the same mode. The repeat distance between molecules which vibrate totally in-phase is the phonon (deBroglie) wavelength, it is written in terms of the crystalline unit-cell dimensions and defines the characteristic momentum of the phonon, q. The

480

frequency of a dispersed mode is a function of its momentum and so is the displacement tensor, which could be rewritten accordingly as Bi(q). Some modes are more dispersed than others; translational modes disperse strongly, librational modes significantly less-so and internal molecular vibrations hardly at all. The full range of q values is available within the first BZ and it is one of the great strengths of neutron scattering that modes at q > 0 are accessible. (Since the photon has no mass it transfers no momentum and only interacts with BZ centre modes, q ~ 0). Unfortunately phase coherence between molecules cannot be directly explored through the INS spectroscopy of incoherent scatterers, or powdered samples. However, provided that the neutron has sufficient energy to excite the transition any momentum transfer value will elicit some spectral response. As we saw above the strength of this response depends on the number of particles in the oscillating system and the oscillator mass. As the number of particles increases the oscillator mass also increases. The scattering intensity would fall if not precisely compensated by the rising crossection. In dispersed modes, the q-range over which the mode is dispersed gradually and evenly becomes filled with bands until, in the limit of very large numbers, a continuous, smooth, distribution is obtained (see Fig. 2). This can be followed theoretically using a 'beaded-chain' model and demonstrated with the INS spectra of the low molecular weight alkanes. However, dynamic strcture factor, S(Q,c0), is determined by the number of final states as a function of the energy, not Q. This is usually written as the phonon density of states per unit energy, g(c0), often referred as DOS. f~x g(co), dco - 1

(18)

It must not be imagined that g(c0), because of its width, is devoid of sharp features. Indeed they are rather common and occur wherever the slope of the dispersion curve falls to zero, i.e., dav'dq = 0. Where this occurs, large numbers of oscillator states are cramped into narrow regions of energy, creating (van Hove) singularities [17] and the INS response increases sharply. This is even more dramatically demonstrated in the case of acoustic phonon branches, these phonons have a maximum propagation frequency and beyond this frequency the DOS falls abruptly to zero producing a very well defined edge, or cut-off. These edges are readily identifiable on high-resolution spectrometers and can be reported with good accuracy. For these features to be observable in INS spectroscopy the correct alignment, discussed above, must be obtained. Powders are ideal for such work since they present every orientation of B to the

481 scattering vector Q. Powdered samples yield almost as much information as single crystal samples and they are much easier to prepare. There are several, more or less sophisticated, methods of incorporating powder-averaging effects into the response function, see [1 8].

S(Q,(oi) = (QQ" Bi (q). exp(-QQ" A))Powder

(19)

At the crudest level is the isotropic approximation, which is almost always representative and is completely adequate in some cases.

Q2"TrBi ( Q2"TrA ) S(Q, co~) = -----f--- g(co,), exp 3

(20)

Three important consequences follow from the powder average; first, all vibrations have some component parallel to Q and are therefore observable; second, the intensity of each mode is weakened, and; third, combinations are observable. Simple combinations are two quantum events, which occur when Q has components along two atomic displacement directions of the same scattering atom. Combinations and overtones are members of the 'multiphonon' family of scattering events that involve the exchange of more than a single quantum of energy with the sample. In the case of undispersed molecular vibrations (i.e. very narrow DOS) these features are well described by dynamical structure factors similar to those for the single phonon events [19]. However, the broad single phonon DOS caused by dispersion leads to even broader spectral responses from the multiphonon events. Even the simplest two quantum overtone, whose shape is the convolution of the original g(03) with itself, has lost a great deal of spectral structure. A common treatment of multiphonon scattering is to regard its integrated contribution as a flat 'background', which is then subtracted. Such an approach is rarely justified when data over large ranges of energy transfer are available. An alternative method involves a selfconsistent calculation of the single and multiple terms (for details, see the next section). Some Instrumental Considerations: We have previously touched on some of the practical considerations, like sample temperature, that limit the generality of the theoretical forms outlined above. It is also appropriate to focus, briefly, on some aspects of the ISIS instrumentation used in collecting the neutron spectra discussed here [17]. The instruments fall into two general classes

482 conventionally described as direct-geometry (e.g. HET and MARI) and indirectgeometry (e.g. TFXA). Direct-geometry spectrometers have access to broad ranges of (Q,c0) and can measure the values of S(Q,c0) at constant energy of incoming neutrons. Well-established extrapolation procedures provide analysed results in terms of a spectral response proportional to the DOS. The indirectgeometry spectrometers used in the ice work form a special set of instruments, which explore a unique trajectory in (Q,c0) and, for a given co, there is no significant range of Q available. However, on these spectrometers every final state of a given oscillator-mass generates the same spectral response, irrespective of co. Therefore, each system of oscillators gives a signal proportional to its DOS, weighted by the DW factor. However, the DW factor remains almost constant over the restricted co range of the DOS, at least for the samples discussed here. Therefore, the observed results, although strictly S(Q,c0), are locally dominated by g(0~) and closely conform to its overall shape. 1.2. Lattice dynamic calculations The vibrational spectrum measured by neutron scattering can be simulated by either LD or MD methods. In the LD approach the calculation begins with the adiabatic approximation which enables us to treat the solutions of the electronic problem as interaction potentials of the nuclear problem. Although modem LD programs can adopt effective potentials, or use several force constants for longrange interactions (such as the second and third neighbour interactions), they are usually restricted to nearest neighbour interactions only. This is not usually a limitation, however, since the nearest neighbour interactions play a dominant role. Moreover this simplification in the calculation results in good quality estimates of the DOS without resort to powerful computers. By assuming harmonic forces and periodic boundary conditions, we can obtain a normal mode distribution function of the nuclear displacements at absolute zero temperature (under normal circumstances). The problem is then reduced to a classic system of coupled oscillators. The displacements of the coupled nuclei are the resultants of a series of monochromatic waves (the normal modes). The number of normal vibrational modes is determined by the number of degrees of freedom of the system (i.e.-~3N, where N is the number of nuclei). Under these conditions the one-phonon dispersion relation can be evaluated and the DOS is obtained. Hence, the measured scattering intensities of equations (10) and (11) can be reconstructed. The method frequently used to calculate phonon dispersion curves from a potemial model can be found in the literature [20] and only brief details will be given here. A potential model describes the potential energy of a crystal. It is assumed that the internal coordinates (basis atom positions) of the crystal are at

483 their fully relaxed positions. The potential energy of the system is expanded, as a Taylor series in atomic displacements. The total energy, V(R), of the system (here R = (rl, r2, r3,.. r,) and denotes the collective nuclear positions) is:

1 dV(R) 1 ~a' d 2 V(R) dradra, V(R + dR) - V(R) + ~ Z dra dra + ~ Z dradra, a a

(21)

where the indices run over all atoms a, a ' and a given position vector can be decomposed in the conventional Cartesian frame, r, =(Xa, y~, Zo). Within the harmonic approximation the expansion is truncated at the quadratic term. For a fully relaxed system the net force on each atom is zero, (i.e. V'(R) = 0) and the coefficient of the third term is the force constant: d2V(R)

F.., (R)= - ~

(22)

aroa,-o,

The equations of motion are constructed from this model and plane-wave solutions are sought. The frequencies of vibrational modes with a given q can be determined by the diagonalisation of F(q)X(q):

(23)

[F(q)X(q)- f2(q)], e(q) = 0

Here the force constant matrix, F(q), and the geometry matrix, X(q) (the matrix of atomic positions and masses) combine to give the dynamic matrix of the crystal. The eigenvalue matrix, f2(q), contains the oscillator frequencies, c0i, and the matrix of eigenvectors, e(q), is obtained directly. The order of the matrices is determined by the geometry, 3n x 3n where n is the number of atoms in the unit cell. Unit cells that contain larger numbers of atoms will require larger matrices to be diagonalised at each q. The calculated S(q,c0) is then given by S(q,co) ~:

1

co,(q))

(24)

a,!

Z le?(q)l - 1

(25)

Q

This is evaluated over all relevant q values see below. In the case of powdered materials S(Q,c0) is obtained directly from equation (19) by repeated evaluation

484 over all directions in space. through

This result is then compared with experiment

S(Q, co) - ISmeasured(Q, co)-background-multiphonon]

(26)

The background scattering (e.g. from the sample container and cryostat) can be measured separately but the calculation of multi-phonon contributions from equation (11) might not be straightforward. Alternatively a number of selfconsistent iteration techniques have also been used, for details see ref. [21]. In most of our work we have compared the observed (S(Q,c0) :=>) g(o) obtained from equation (20) with that calculated from 1 g(co) = ~ - ~ ~

6(09- co~)

(27)

Here we should note that the dispersion curve calculation has provided all the information required to obtain the response from a single crystal sample aligned along a specific direction in Q. Indeed, if such an experiment were realistically feasible it would be the preferred technique. This is because the dispersion curves would be measured directly and the detailed information about the force field could be extracted. However, this is often not practical, at least for the exotic phases of ice and powdered samples were used. For ice Ih, single crystals are widely available (but a large crystal of ice Ic has not been obtained), after many attempts [49,55], reliable dispersion curves have yet to be obtained. This is due to the proton disordering in the structure of ice Ih and hence all the optic modes are localised. As an example of the information lost by exploring the DOS of powdered samples we compare the calculated dispersion curves of ice Ih and Ic. They have quite different dispersion curves in the translational region due to their different symmetries, but these are inaccessible to neutron dispersion measurements due to the proton disorder in the structures. Only the DOS can be measured. As a result, the detail of the information in the dispersion curves is lost, or at least degraded, by comparing only the DOS. From both experiments [22,55] and calculations, we have found that the two ices share an identical spectrum as shown in Fig. 3. This is because they share the same local structure in their lattices (the tetrahedral symmetry) and the same local force field. If the integration over the first BZ is incomplete (i.e. if too few q-points were used), there would be a considerable difference between calculation and observation spectra.

485

1/cm

350 "~

meU

Ice Ih

300 "1

4o 35

250 -~

30

200

~

100 - ~

~

25

w(Q)

_

50 -

-5 ! i I i i i i i i

O-

0.5 A

UNIT 1/cm

0.0 I

i

1.0 K

O (l/A)

I

I

i

!

I

i

1.0 M

i

-0 0.0

I

300

UNIT m e U

lee Ic

25O

30

200

25

w(Q)

20

150

15 100 -

10

50

5

0

to 0.5 A

0.0 I

0.0 Q (l/A)

0.5 M

K

0.0

g(w)

9

50

100

150

200

!

250

,

!

i

300

350

300

340

ENERGY TRANSFER E 2 - E1 (1/cm)

+
50

100

150

200

250

ENERGY TRANSFER E 2 - E1 (1/cm)

Fig. 3. The upper two diagrams are the calculated dispersion curves for ice Ih and ice Ic based on a simple LD model containing O atoms only. An O-O-O bending and an O-O stretching force constants, G = 0.33 eV/Rad 2 and K = 1.1 eV/A 2 were used. The diagram second below shows three curves of g(m) for the three particular reciprocal directions in ice Ic. The lowest diagram is the completely BZ integrated g(m) for both ice Ic and Ih. It differ considerably from the curves above, indicating the incomplete BZ integration can be misleading.

486

2. M O L E C U L A R DYNAMIC SIMULATION OF NEUTRON SPECTRA Molecular dynamic (MD) calculations are based on the numerical integration of Newton's equations of motion for a many-body system. Using very small but finite time steps, the real positions and the associated velocities of all of the atoms in the system are followed. New atomic positions and velocities are determined from the numerical solution of the equations of motion with suitable potentials acting for fixed durations. For a lattice system a modest number of atoms covering several unit cells can be used to create a simulation cell. Moreover, if periodic boundary conditions are imposed, properties that depend upon distances much larger than the simulation cell can be calculated. The mean kinetic energy of all of the atoms in the system represents the system's temperature. Scaling the velocity co-ordinates and equilibrating can be used to control this temperature. With the increasing power of computers, this method of classical mechanical simulation has become very popular. Unlike LD, it simulates the dynamical processes at given temperature and pressure. It is therefore a widely used tool for the investigation of a range of microscopic properties in liquids and solids, such as structure, vibrational dynamics, diffusion and phase transitions. Its versatility is demonstrated by the wide variety of work related to water and ice. Of particular interest here is its use in the simulation of the INS spectra of lattices. MD calculations readily relate the velocity, v(t), of a given atom, at time to = O, to that of the same atom at some latter time, t. The development of this correlation as time passes provides access to the single-particle velocity autocorrelation function, VACF(t).

VA CF(t) - (Va (t)" V~ (0))

(28)

In a harmonic crystal the DOS is the real part of the Fourier transform of VACF(t), the mass weighted power spectrum Z(c0), which is almost identical to the DOS calculated by LD and measured by INS, the differences will be discussed latter.

_

mo)

/vo o

"

(29)

This makes MD a very powerful tool in the simulation of INS spectra. In the computer simulation process the VACF(t) must be calculated from the starting

487 value V(to), averaged over many initial times to from zero to the maximum. The standard convolution method of calculating Z(c0) is very slow indeed. An alternative method of calculating the power spectrum is available in the fast Fourier transform technique [23], because:

f < Va (t). Va (0) >. exp(icot), dt

~ V a (t) . exp(-icot) . dt[ 2

-

(30)

We therefore have:

Z'(ol)-~f~ a

m

Va(t )

a 0

exp(-icot) dt "

"

12

(31)

Using this expression enables us to make use of the fast Fourier transform algorithms, which provide an enormous gain in speed over the equivalent autocorrelation method of equation (29). The LD and MD simulations for ice Ih using the TIP4P potential are shown in Fig. 16 and 17. At low temperatures, LD and MD results share almost the same features in the translational and librational regions, indicating that both approaches are valid for DOS calculations. However, considerable differences remain, for instance the MD simulation obtains forces from the first derivative of the potential function (i.e. F~ = dV(R)/dr~) and LD uses its second derivative (i.e. restore force) which can be obtained from equation (22). Comprehending their differences in the simulation processes and the methods of obtaining a DOS, their distinct advantages might be combined. This will become essential if we are to simulate a range of effects and make further progress in our understanding of the vibrational dynamics of ice and water, as we discuss below.

2.1. Equipartition theorem Because of the classic approach involved in a MD simulation, energy is equipartitioned among all the vibrational modes (equipartition theorem [24]). In a system with N elements, the total kinetic energy, Ek, is the sum:

Ek _ _~13 Nk 8 T = ~ lma 2 v j (O = _~k s T f ~ ~ 6

- coj )do)

(32)

For a harmonic system, we also have (from equipartition theorem): 1

k , T - -~m,,(U2,~(O).co 2)

(33)

488

From these equations, we can calculate the averaged vibrational amplitude Ua, i(O) for each mode, i, or, if we calculate the averaged velocity, the instantaneous temperature of the system can be obtained. This is important if we wish to maintain the system at a constant temperature. A suitable thermostat (for details see ref. [25,26]) can be used for this purpose. Because of equipartition all frequencies should be equally weighted (i.e. W(o)) = 1 for all o)). This clearly contradicts our understanding of dynamics from quantum mechanics (QM) theory. Since QM is essential to the treatment of certain phenomena, quantum corrections must be introduced, regardless of whether structural or dynamic simulations are required. This discrepancy is illustrated if we calculate the total kinetic energy of the system quantum-mechanically: oo

Ek -

hco. W(co). 6(co - co,.)dco

(34)

i

and the factor W(0)) is given by a modified Boltzmann distribution [27]"

1/

1

/

W(co) = -~ + exp(hco/k.T)- 1

(35)

Thus the population of modes with energy above kaT would decay very rapidly. Most of the librational and all the intramolecular modes would be suppressed leaving only their zero-point motion, as expressed by the factor of 89 in equation (35). The RT motions of water are, therefore, effectively limited to molecular diffusion and those vibrations in the low energy part of the translational region. Hence the thermal displacement should be U a (t)

- ~

U~,, (0). sin(co~t + qoi). W(co~)

(36)

i

Which clearly differs from the classical MD result" U~ (t) - ~

U~,~(0)-sin(co~t + ~o~)

(37)

i

This difference would seriously hamper any attempt to simulate some of the dynamical aspects of water, such as diffusion and phase transitions. We should note that, if the amplitudes of zero-point motion are sufficiently small, the use of rigid-molecule models in the MD simulations of water is appropriate. This

489

approximation eliminates the internal degrees of vibrational freedom, where quantum effects are very strong. 10 8

~)~

50K

6 20

4

~

2

100K 1

0

0

20 4

200K

2

0

0.0

10.0

20.0

30.0

40.0

ENERGY TRA NS FE R (meV)

50.0

0

50.0

70.0

90.0

1 1 0 . 0 130.0 150.0

ENERGY TRANSFER (meV)

Fig. 4. The MD simulations for a 512-water cell using TIP4P [6] at T = 50 K, 100 K and 200 K in the translational (lift) and librational (right) regions. The spectrum for the highest temperature is smoother than others. But the effect is less dramatic than that observed experimentally.

2.2. Anharmonicity and temperature effects The MD calculated DOS is, to first-order, independent of temperature and only small anharmonic effects appear at non-zero temperatures. These anharmonic effects arise from the fact that most potential functions are not parabolic. As the displacements increase in magnitude the molecules explore non-parabolic regions of the potential and the overtone frequencies with perfect integers of the base frequency COo, such as 2(o0, 3C0o..., due to the Fourier expansion of the non-parabolic potential function. Experimentally, however, this is not observed and the overtones rarely fall at exact integer values. Overtones shiit to lower (or higher) frequencies dependent on the curvature of the potential function. By limiting the exploration of the well to the area around the local minimum the effects of anharmonicity can be

490 reduced. Indeed the smaller the range of this exploration the more harmonic will be the response. Equation (33) shows that the higher the vibration frequency COo,the smaller the vibrational amplitude U would be. This means that the anharmonicity effect is less significant for the intramolecular modes in the MD simulation, and we must pay particular attention to the low frequency modes in the simulations for water. Anharmonic effects increase with temperature since the amplitudes, U~, increase as U~ = (kBT)l/2/co~. The atoms displace further from their equilibrium positions, producing large anharmonic effects though the mechanism mentioned above. This is very different from the other temperature effects observed from real neutron scattering experiments, such as the effects from DW factor and phonon-population factor, n(c0) for energies less than kBT which also results of multi-phonon scattering. Hence the INS measurements at high temperatures (T > 200 K) become so much contaminated by these factors and hence their contributions are very difficult to be separated from the one-phonon term. On the other hand, at higher temperatures, the MD simulations give an overall increase of intensity for all the frequencies with small additional anharmonicity. Finally, the combination features in the INS spectrum have not been observed from the standard MD simulations. 2.3. Size Effects and the Calculated DOS The calculated velocity auto-correlation function only gives vibrational modes at q = 0, because of restrictions on the periodic boundary conditions. Other modes in the first BZ can, therefore only be introduced by the "zone-folding" process of super-cells (see Fig. 5). The larger the super-cell, the more q-points that are included. The super-cell size affects the quality of a DOS obtained from any integration across the BZ and such considerations are especially important for strongly dispersed modes (see section 5.2). There are very few MD simulations with more than 300 molecules and precise estimates of the q-points represented in the cell are needed. For instance, in order to include the boundary points, the size of the super-lattice would need to be at least 2x2x2 (= 8) primary cells (i.e. folding once in each reciprocal direction). The result of this calculation can not be considered as an "accurate" representation of the integrated DOS, since it contains only q-points at BZ centre and boundaries. A typical 256 molecule cell is equivalent to a super-cell of 4x4x4 (= 64) unit cells. This cell gives an additional q-point in the middle for each dispersion curve, or 3 wave-vectors one for each reciprocal direction, as shown in Fig. 5. We believe that the minimum requirement for MD calculations of the DOS is a super-lattice cell of 5x5x5 (= 125) unit cells or 512 water molecules for ice Ih. Ideally, 8x8x8 super-lattice cell (= 512 unit cells or over 2000 water molecules for

491

ice Ih) is more appropriate if the computing time is not restrictive. This super-cell provides 5 q-points on each dispersion curve (or a total of 5x5x5 q-points in the first BZ) for comparison with the measured DOS. A typical LD calculation integrates over 50x50x50 q-points in the first BZ. This represents a thousand fold improvement in the mode integration, which is equivalent to a super-lattice with a million water molecules.

l

I

....

o

,

~

q=,=/l

?-?q L_A_A 1.0

0.5

-0.5

-1.0

Fig. 5. Schematic illustration of the zone folding effect (above figure). Even when the super cell increases to 4 unit cells in each direction, which is a total of 4x4x4 primary cells (equivalent to a 256-molecule super-cell, see the size C), there are still only 3 wave-vectors allowed in each direction. As shown at the bottom of the diagram, the wave-vector at q = 0 (i.e. d = oc) is not included.

492

In order to demonstrate the size effect, Burnham [26] has made a series of MD calculations with different lattice-cells, having 64, 128, 256 and 512 water molecules at 100 K. The potential function used was again TIP4P. As one can see from Fig. 6, the size effect is quite dramatic. The 64 and 128 water cells give a DOS with highly structured noise. In more complex systems, some of these features could be mistaken for real peaks. Indeed, in the case of ice the noise at 28 meV was otten believed to be one of the two peaks observed in the INS spectrum. This incomplete sampling of the BZ is also demonstrated LD simulation in Fig. 3.

30

512 20~

10

20i 03 u.I F-I-03 I.i. 0 >I--03 Z u.i CI

256

4

10

\___

2 0

0

I

20

128

4

0 20

64

10

01._& 0.0

10.0

20.0

30.0

~ 40.0

ENERGY TRANSFER (meV)

50.0

0 I,,~-" 50.0 70.0

90.0

110.0

130.0

150.0

ENERGY TRANSFER (meV)

Fig. 6. MD simulations for ice Ih with different sizes of super-lattice cells, 64, 128, 256 and 512 water molecules using TIP4P potentials. The calculations show that intensities for 64 and 128 molecules are very "noisy". The 512 molecule cell shows a good agreement with LD simulation result see Fig. 16, indicating the BZ integration is about acceptable with at least 512 molecules.

493

2.4. The energy resolution and intensity statistics In order to compare with the experimentally measured spectrum, one would ideally like to have the spectra simulated with equal or better energy resolution and with a good statistical reliability for the predicted intensity. Therefore, a suitable estimate of these quantities is important, in order to minimise the output to a manageable level (i.e. to output the trajectories and velocities as less as possible). According to the sampling theorem [23], the smallest time step for the Z(03) calculations is a factor of n' times larger than MD simulation step which is determined by the maximum frequency, C0Ma~, of the system to be simulated. Hence the appropriate time step n 'At is given by the Nyquist sample rate, 203Max, as:

n' At -

(38)

03Max

In water and ice C0M~corresponds to intramolecular vibrations at about 500 meV (or 4000 cml). Hence, we estimate that n 'At is about 4 femto-seconds (fs) or n' = 10 for a MD simulation step of 0.4 fs for simulations of non-rigid-water models and n' = 20 for the rigid-water models. These values of n would allow enough MD steps to resolve the highest frequencies required for the simulations. The maximum time period in a MD simulation is T (=M'At, where M' is the total number of steps in the simulation). This is determined by the required energy resolution of the resulting spectrum and is given by the Fourier transform of auto-correlation functions [23]"

A E - h(Aco) -

h

2zT'

9

7/"

T - M ' At - ~ Aco

(39)

A typical INS spectrum is obtained with instrumental resolutions, AE/E, varying from 1 - 3 %. In the case of TFXA, the resolution in the low energy region is 1 % . If we wish to resolve the modes at ~ 30 meV, this requires a energy resolution, AE, of-~0.3 meV; which corresponds to M ' = 10,000 Md:) steps using 0.4 fs per step.

3. W A T E R - W A T E R POTENTIALS Over 100 years ago, Rontgen and co-workers [28] were already aware of the range of the anomalous properties of water. They postulated that the liquid was

494 an aggregate of two different types of constituent. The first type was ice like and the second type resulted in a decrease in the volume of the solution. This idea was later developed imo the physical models, which were popular in the 70's (e.g. flickering cluster model [29], and ice-like continuous model [30]). They provided a graphical explanation of the abnormal properties of water. With the development of modem computers these physical models were replaced by model computations based on suitable water potential functions. There are many advantages to these modem methods; for instance in MD simulations, a range of time dependent dynamic properties can be calculated on time scales from fractions of femto-seconds to thousands of pico-seconds (ps). Hence, the last twenty years has seen a sharp increase in the simulation of the properties of water and ice using a wide variety of water potentials [8,31,32]. These simulations demand better and more accurate water potentials to simulate complex phenomena, such as the vibrational dynamics, phase transitions and transport properties. The potential functions used in these calculations have gradually evolved, developing from very simple LennardJones type with 3-point charges (e.g. BF [34]), 4-point charges (e.g. ST2 [5]), poladsable potentials (e.g. SK [35] DC [36] and NCC [37]) to the very complicated anisotropic multiple polarisable potential (ASP [38]). The process was also associated with a gradual increase in the anisotropy of these potentials.

q

+q

+q

+q

-

q

-2q 3/3

3/4

+q 4/5

Fig. 7. Schematic diagrams of the point-charge arrangements in the classic pairwise potential functions. On the left-hand-side and in the centre is shown the 3/3point (the open cycle represents O and the large solid cycle is H), 3/4-point charge distributions, which are two categories of the 3-point charge models. On the right-hand-side is the 4/5-point (or 4-point) charge model. The early work considered water molecules as rigid entities. Both the attractive and repulsive parts of a core potential are needed and these were constructed in two principal ways. In the first approach the components are obtained by ab-initio quantum mechanical calculations for the ground state energy of the water dimer (e.g. MCY at Hartree-Fock level [4]). The analytical form of these potentials was fitted to the calculated energy surfaces (details of these potentials are given in ref. [8]). This resulted in unusually long O-O

495

distances and relative soft curvatures at the potential minima and made these models less successful in simulating the bulk properties of water and ice. In the altemative approach, a few physical parameters of bulk water, such as the measured O-O distance, the binding energy and the dipole, were fitted (e.g. TIP4P [6]). These potentials are best considered as effective potentials and are much more suited to the simulation of bulk properties. They have seen wide use in the simulation of the structure and dynamics of water and ice. Other potentials are all more or less similar but show variation in the values of the parameters they use these values are detailed in Table 1. The arrangements of the point-charges used in some of these potentials are show in Fig. 7. ,

,

,

~

,

,

Z m

0

I0

R, IX]

t --

,

-12d -~o'

,

.

d

.

.

.

sb

.

.

i;~o

.

!

fad

240

ANGLE

Fig. 8. On the left are shown plots of total energy vs O-O separations (upper) and the relative water dipole-dipole angle for the MCY potential. The relative dipoledipole orientational configurations for ice Ih is shown in the diagram. The type-B configuration at 180 ~ corresponds to the lowest energy in the V(r) plot, while the type-A is at 60 ~ type-C is at 120 ~ and type-D is at 0 ~

496

r

O9

c~

o

0

0 0

c~

..o

o

c~

..o

o r

0

0

~

~

9

c~

ox

'

,

~

~

'

~

--~ o n

s

~

r

'

,

mlm

! l

~

,~t- c q

cq

ox

~ ~,,,i

o

0

~c~

N

.~

I~

o

o

c~

o

~

~

~

'

o

~

~

Ox

,~t- c q

oo

Ox oo

,

9

oo

o

~o

,

c~

'

ox

'

Ox oo

c~

o

g,~

0

~

o

~ . -

m

9~

497 The curvatures across the energy minima of a number of the classic pair-wise potentials have been calculated (see Table 1). The force constants, of equation (22), have been extracted for the four proton configurations, namely; A, B, C and D for ice Ih shown in Fig. 8. These values can be compared with the force constants used in the classic LD calculations of section 5.1. As one can see from Table 1, the ratios amongst the four force constants for the configurations for the potentials listed are all less than 1.3 which is considerable small than the value of 1.9 required for reproducing the observed INS spectrum for ice Ih as (more details discussed in section 6.1).

3.1. Water clusters and polarisable potentials In recent years, there has been considerable research activity into the investigation of isolated water clusters, which are believed to be the basic building blocks of bulk water. An understanding of the structure and dynamics of clusters is, therefore, of considerable importance in the process of constructing more accurate water potentials. Moreover, the recent development of ultra-high resolution IR (vibration-rotation tunnelling) spectroscopy has provided high quality spectroscopic data on small clusters [39,40]. High-level quantum chemistry calculations and diffused Monte-Carlo methods [41,42] have also been used to interpret this data. The advantage of studies of the different sizes of water clusters is to identify the many-body interaction contribution to the potential function. The total energy of a water cluster depends on several contributing terms; e.g. a trimer consists of 2-body terms and a 3-body term and the tetramer has an additional 4 body-term and so on. From studies of the water dimer, trimer, tetramer to hexamer, the various many-body contributions can be separated and accurate estimates made of the 1-, 2-, 3- and 4-body contributions for these clusters and larger systems. These calculations frequently show that the 3-body term is very significant and contributes -~20% towards the total energy of the cluster. The 4-body term contributes ~5% and all higher orders, taken together, is less than 5%. The many-body terms constitute ~30% of the lattice energy. These additional terms have not been properly considered by the classic pair-wise potentials, which by their nature account only for the 2-body terms. In this type of potential, a fixed dipole moment was used (e.g. ~2.2 Debye (D) for TIP4P). Hence they are unsuitable for the simulation of water vapour or mixed vapour and liquid. By introducing polarisability into the potential the local electric field given from surrounding water molecules generates an additional dipole moment. Because the Bemal-Folwer ice rules [34] constrain the allowable orientations of the nearest neighbour water molecules, viewed from a central (i.e. a target) molecule, the electric field generated by each molecule cannot cancel as illustrated in the lower diagram of Fig. 9. This produces a strong effective field

498 which polarises the changes on the target molecule and gives rise an additional dipole for the target molecule. For molecules much further away from the target molecule, their contributions to the electric field are less, since the orientations of these molecules are more random and produce better cancellation. Under these conditions the polarisation effects are also short range, which is consistent with the experimental evidence discussed the later sections. This polarisation effect would probably produce a large orientational variation of the potential.

C

"""OH a

"""O H b

"~H c

P

P d Fig. 9. Schematic diagrams of the three types of polarisable potentials. The left-hand diagram shows a point polarisability model (e.g. SK [35] and DC potentials [36]). The centre diagram shows the polarisation on the two O-H bonds (e.g. NCC potential [37]). The fight-hand diagram shows the all-atomic (or three-) polarisation models (e.g. Bemardo et al [44] and Burnham [26]). The lower diagram schematically illustrates the relative orientations of molecular dipole moments of the four nearest neighbour molecules would in possible to cancel out due to the ice rule and give rise a strong local field. The success of this class of water potentials is that the additional dipole moment generated from the local electric field and polarisation successfully accounts for the difference between the gas phase value (-~1.8 D) and that of the liquid or solid phase (-~2.3-2.8 D). Moreover, an additional energy produced from the polarisation agrees well with ab-initio results from a number of polarisable potentials (e.g. SK. DC and etc). We believe that polarisable potentials correctly account for many-body effects in total energy and dipole

499 less effective in calculating the energy from polaraisation (only 14%) than recent additions, the new ones, such as SK, DC and ASP, are capable of providing polarisation energies between 25-30% for water, which is very close to the ab-initio results. However, the NCC potential [36] produces a polarisation energy in excess of 50%, which is well above realistic values.

3.2. Validating water potentials As we have described above, currently a large number of water potentials are available. Choosing the appropriate experimental data for the validating tests is therefore important. One of the conventional methods was to reproduce the partial radial correlation functions GHH(r), Goo(r) and Goa(r) for water obtained by neutron diffraction (either by using the isotopic substitution method or by combining X-ray or electron diffraction data). In general, the MD simulations of water structures using these potentials give good agreement with the ones obtained experimentally. It is often seen that the simple rigid point charge potentials, such as SPC~ [65], give almost identical results to the very complicated polarisable ones as illustrated by Dang and Chang [36]. In fact, the uncertainties (or errors) introduced in the partial correlation functions by the data reduction from the measured diffraction data were much greater than the The test of a w a t e r potential

[

,

I MicroscopicpropertiesI I

Strcturefactors

I

I Neutron diffraction

SHH(Q)

Soo(Q)

SHo(Q)

I

Fourier transform~ GHH(r) / Goo(r) | GoH(r) /

I

I

1~ Bulkproperties I

I I Vibrationaldynamics1

[ Diffusioncoefficient

I

--IHeatcapacity -1 4"(:3density maximum

Inelastic

neutron scattering g(w)

land etc.

500

errors (or differences) between the fitted data by the MD simulations and the measured data. The systematic errors in the data treatment rise from two main sources: inelasticity and non-equivalence of H and D in the isotopic substitution (for details see ref. [45]). In addition, the measured data were made in the reciprocal space (i.e. Q space) which requires Fourier transformation of the measured S0{Q) to the real space variable G,~(r) in order to compare with MD simulation results. Very detailed studies of these classic pair-wise potentials have been made by Finney et al [32] and Morse and Rice [33] some time ago for water and ice respectively. Finney et al showed that although these potentials are in general reproducing qualitatively the radial distribution functions of water, the classic pair-wise potentials remain rather primitive. From their studies, they concluded that the 4-point charge models have stronger angular constraints and produce better water structure than the 3-point charge models, because the 3-piont charge models give rather simple liquid-like structures. Morse and Rice's simulations of ices indicated that most of the classic potentials are capable of stabilising ice structures, such as ice II, XI and VIII [33]. However, their predictions of other physical properties are poor. Hence the structural simulations are insufficient for the validation of a given potential function. This is because the simulations use a potential involving a delicate balance of many competing effects. The short-range interaction, from the core of the potential, competes with the long-range charge interaction, and those terms dependent on orientation compete with distance dependent terms. Small adjustments in the arrangement of the charges result in insignificant changes in the structures finally produced as frequently showed in MD simulations [36]. These collective effects are impossible to isolate from one another and there is, as yet, no consensus as to which particular potential is most acceptable for the various simulation tasks. Hence, the exclusively structural investigation of water potentials is unable to provide clear guidance and is insufficient to validate a given potential function. Further progress on the determination of the water potential needs more precise experimental data to validate the potential functions. Since high resolution neutron spectroscopy became available in recent years, it has provided additional experimental data for the benchmark testing. Other the other hand, because both MD and LD provide the power spectrum (a sum of the normal modes) using the first and second derivative of the potential, the simulation results can be directly compared with the experimental spectrum without involving the Fourier transformation and other uncertainties (or errors) associated with diffraction data. Hence, simulating the INS spectra for a large variety of crystalline and amorphous phases is of considerable advantage in the process of further examining the potential function developed.

501

4. NEUTRON VIBRATIONAL SPETRA OF THE EXOTIC ICES Since dedicated neutron sources for scientific research became available in the 1960's, neutron scattering techniques have been widely used for the investigation of the structure and dynamics of water [46,47] and ice [48,49]. However, the early attempts at the measurement of the vibrational dynamics were compromised by instrumental limitations of neutron flux and poor energy resolution (see Fig. 10). They lacked sufficient detail in the translational region (< 40 meV), which is a crucially important area in providing information on hydrogen bonding in ice and in testing the accuracy of the existing water potentials by LD and MD simulations. The current development of intense, pulsed neutron sources such as ISIS has provided impetus to this work. Specifically the range of high resolution inelastic scattering instruments, such as TFXA, HET, MARI and PRISMA at ISIS (UK) [50], has made possible an accurate study of the dynamics of ice. These instruments have far superior resolution to any other available spectrometers in the world. The high neutron brightness of the source and highresolution of the instruments reduces backgrounds and improves the signal-tonoise ratio, to negligible proportions for scattering samples such as ice. The low, flat, background obtained from the sample container is measured separately and subtracted. As a result, we have been able to obtain the spectra of ices with unprecedented accuracy. We have shown, in later sections, how precise INS measurements of the DOS provide the most stringent means of testing the model potential functions that lie at the heart of any LD or MD simulation. In the last a few years, we have systematically studied the vibrational dynamics of a large verity of phases of ice using above instruments at ISIS. These spectra were obtained at very low temperatures (< 15 K) on the recoverable high-pressure phases of ice and a few forms of amorphous forms of ice, in order to reduce the Debye-Waller factor and avoid multiphonon excitations. Hence the one-phonon spectra, g(co), can be extracted from the experimental data for the theoretical simulations. Ideally, the measurement of the g(co) for normal ice (Ih) at different pressures would provide the information about the hydrogen bond interaction V(r) as a function of r. The difficulty with such measurements is that the structure of ice Ih readily transforms at only modest pressures, less than 3 kbar, and below this pressure there is little change in the hydrogen bond lengths. Hence, pressure measurements have to be performed on other phases of ice, such as ice II, III, V, VI and VIII in order to cover an extended pressure range. On the other hand, in these ice structures, the O-O distance varies even at the ambient pressure from 2.76 A to 2.965 A and the 'tetrahedral' O-O-O angle also changes, from 83.8 ~ to

502

126.2 ~ depending on the phases. These are significant distortions away from the values of ice Ih, with its O-O distance of 2.75 A and true tetrahedral angle, 109.47 ~ and provide broad scope for the measurement of the potential surface. 7i

....

I . . . .

I

....

I ' ' ' ' 1

....

I

D

FZ LLI

Z

NEUTRON ~ . , ,

0

I,

10

....

I .....

20

I , . , , I

30

....

40

I

50

ENERGY(meV)

Fig. 10. Comparison of typical IR, Raman and neutron spectra for ice Ih shows that there are significant differences in emphasis in the different techniques. The IR and Raman data show the main peak at 27 meV. The early neutron spectrum [48] shows a lack of detail due to poor resolution and intensity. The new neutron spectrum measured on TFXA clearly shows two peaks at 28 and 37 meV. The higher energy peak is twice as intense as the lower energy peak. The structures of the range of exotic crystalline phases of ice have been, for the most part, well known for many decades [9] and provide a suitable framework for the theoretical modelling. Moreover, by suitably choosing the

503

appropriate phase of ice, data for different values of the O-O separation distances from different structures are available.

L

'

I

i

I

I

I

J

~

i

I

i

i

i

I

i

I

a

I

=

i

i

I

I

i

I

r

II q

r o i

p

I

"i

r

9

=l=..

r-

~

t

I

"

i

i

20

i

I

i

i

i

I

,

,

,

I

i

40 60 80 100 120 energy fronsfer (meV)

i

i

140

Fig. 11. IINS measurements of HAD (top) and vapour deposited LDA (middle) ice (bottom) using TFXA on ISIS. The spectrum for ice Ih is also plotted for comparison. All measurements were made at temperatures below 15 K. In the pressure range below 25 kbar there are at least 10 different crystalline phases and a few amorphous forms of ice. Most of the phases can be measured at ambient pressure using the technique of 'recovery'. In the recovery process metastable phases of ice at high-pressure are quenched at liquid nitrogen temperatures and the structures are retained when the pressure is released. The structures of these phases would remain unchanged at temperatures below 120 K, but relax slightly to give small changes in their densities. The formation of the amorphous ices is complicated due to the uncertainty of their structures. In general, the community believes that there are two main classes of amorphous ice, the high density (HDA) and low density amorphous (LDA). The HDA can be produced by pressurising ice Ih at about 10 kbar at temperatures below 120K. This form of ice can be recovered using the technique described above and has density of 1.17 g/cm 3 in the recovered state. The classification of the LDA is

504

slightly complicated, because its structure can vary slight dependent on the preparation processes (or techniques, for details see ref. [51 ]). One of method is recovery from higher density crystalline phases of ice such as ice VI or ice VIII, or HAD [51 ]. The other way is to deposit water vapour on a cooled surface. We have found that these two forms of LDA have quite different g(o9), hence their structures can not be the same. In addition, a few crystalline phases, such as ice III, IV and VII, are not accessible through the recovery technique and measurements have to be made at the necessary pressures. So far only ice III (having an almost identical structure as ice IX) and VII has been measured under direct pressure [52]. Although this requires the presence of bulky metal pressure cells in the neutron beam, a correct choice of elements (e.g. A1 and a special ZrTi alloy) with high quality background measurements would minimise the problems associated with data reduction processes. The spectra obtained for ice Ih, LDA and HDA, using the TFXA spectrometer at---10K [53] is shown in Fig. 11. Ice Ih is the most common and readily obtainable phase of ice which has now been well studied [14,15,48,49]. Its spectrum has a very simple structure, the translational modes below 40 meV are well separated from the librational modes (or hindered rotations) in the energy region between 65-125 meV (very few system shows similar behaviour and this is due to the large mass difference between O and H). The observed acoustic phonon peak is at 7 meV. The two sharp peaks at 28 and 37 meV are the opticphonon bands and have an unusual triangular-shape. In contrast, only a single feature appears in the IR spectrum, at 27 meV, and the Raman spectrum has an additional shoulder at 36 meV (see Fig. 10). The observation of the two distinct triangular peaks for the molecular optic modes by the high resolution INS measurements represents a considerable challenge to the MD community. Even today, it remains a comroversial subject. The difficulty of the task is seen from the fact that a single optic feature dominates the spectra of other tetrahedral systems, such as Si and Ge. If then one assumes that all the hydrogen bonds in water or ice are equivalem (i.e. that the hydrogen bonding is isotropic and shows no changes with angular variation, see Fig. 9), the simulations produce a single optic peak between 30 and 40 meV for ice as shown in Fig. 3 and later sections. This observation indicates strongly that it must be the anisotropic components of the water potential that causes the optic peak splitting. Throughout this article, we have underlined that from the position of the two bands, the orientational variation is considerably larger than one would normally anticipate. Understanding the features of the INS spectrum holds the key to understanding the dynamical properties of water and ice and the mechanism of the water anomalies.

505 The INS spectra of ice Ic and the recovered LDA ice (obtained by annealing the HDA form at 120 K) have very similar features to those of ice Ih (see Fig. 11 and 10), which indicates that the force field and the local structures of these ices are almost identical. This is despite the significant differences in their longrange structures and symmetries: ice Ih has hexagonal symmetry, ice Ic has cubic symmetry and LDA has no long-range order. However, the spectrum of the HDA ice differs considerably from those of LDA, ice Ih or Ic. This reflects the fact that the HDA local structure has been crushed by the high-pressure and its density has risen to 1.17 g/cm 3. In the INS spectrum of vapour deposited LDA ice, it is the higher energy optic band, at 37 meV, that dominates [51 ], which is considerably different from the spectrum of recovered LDA. This indicates that porosity in the vapour deposited LDA has produced an increased surface area, where water molecules are not completely hydrogen bonded. These molecules may be able to relax to the lower energy configurations in Fig. 9. A similar phenomenon has also been observed in the INS spectra of water on the surface of porous solids, such as Vycor and silica gel [54]. Using a single crystal of ice Ih, INS spectra were measured with Q along the c-axis of the crystal and in the basal plane. The data show that, although there were small differences in the acoustic region (< 7 meV), the intensities of the optic peaks at 28 and 37 meV were independent of the crystal orientation [14,55]. Furthermore, no differences were observed for the librational or intramolecular band strengths, covering the region from 60 to 500 meV. This reinforces our view that the local field determines the spectral features and that this field is the same in ice Ih, Ic and recovered LDA. The INS spectra of the recovered high-pressure phases of ice, measured at ambient pressure, are quite different from ice Ih in the important translational and librational regions [55], see Fig. 12. This is because the local structures have been strongly distorted. The hydrogen bonding in these systems is different and changes the local force field. Little theoretical work has been done on these apart from a few studies of the simpler proton'ordered structures, such as ice II [57] and VIII [56,58], indicating that the distorted hydrogen bonds are considerably weaker than the normal ones. As a consequence, there would be a range of force constants among the hydrogen-bonded water molecules, the two optic peaks vary in position and spread considerably, depending on the local environment of the water molecules in the ice structures. However, the highenergy cut-off for the translational band remains. Except in the spectrum of ice VIII [56], where there is only one optic peak in the translational mode region, at 28 meV (see Fig. 12). This may be due to the fact that ice VIII has a proton-ordered structure and the local dipole configurations may correspond to the weaker interactions in ice Ih [53,55]. The

506

high-pressure measurement for ice VII (a structure is almost identical to ice VIII, but its protons are disordered, the degree of the proton-disordering depends on the temperatures), shows that the high energy peak appears when the protons in the structure is disordered [53].

'

'

'

I

"

'

'

I

'

'

i

'

'

'

I

'

'

'

I

'

'

'

I

'

'

'

I

i

10

,

~/I

ICE VIII

1i

B

~ 6

I E

4

Z

J

J 0

20

dO 60 BO I00 ENERGY TRANSFER (meV)

120

140

Fig. 12. INS measurements of all possible recovered crystalline phases of the exotic ice (H20) were made TFXA on ISIS at temperatures below 15 K. The intramolecular vibrational frequencies occur at higher energy transfer above 200 meV and the HET spectrometer at ISIS was used for this work in order to reduce molecular recoil, this has been described in detail elsewhere [55]. As can be anticipated from the covalent nature of the forces responsible

507 for these intramolecular modes, changes in the external structures (or pressures below 24 kbar) should produce little impact. Such weak effects may be accessible to more sensitive probes, such as optical techniques. The two main features are shown in Fig. 13; these are the intramolecular bending at ca. 204 meV and the symmetry and asymmetry stretching modes at ca. 410 meV. Only small changes were observed between the different phases of ice (for details see ref. [55]). The broad features at ca. 280 meV are combination bands between the bending modes, at 204 meV, and the strong librational bands about 70 meV. I

24

>"

4J .,--I {D

20

'

'

'

'

I

'

'

'

'

I

'

i

,

,

I

'

'

'

'

I

,

,

,

I

,

,

,

,

I

Ice Vl

-

r"

c

F-4

16Ice V

1:3)

E -r-4

c_

12-

Ice II

-IJ

U

m

8Ice IX

121 r 4J :3

m z

4

-

I

,

,

,

,

I

,

,

,

,

I

,

i

100

200

300

400

500

Energy Transfer (meV) Fig. 13. Neutron spectra for a number of recovered exotic phases of ice and ice Ih (H20) were measured using HET spectrometer on ISIS with incident energy of Ei = 600 meV at temperature T - 10 K. The data show very small differences among the different phases, indicating there is little effect to the intramolecular frequencies from the external structures. The high-energy transfer spectrum of ice VIII was again unique. After subtracting the background and multi-phonon contributions, we found a single, strong, feature at 426 meV for H20 (slightly higher than other phases of ice). This band was also significantly narrower (---20 meV) than those observed in the spectra of the other ice phases (typically,-~40 meV). This sharpness are difficult to understand, it may result from the ordered nature of ice VIII and suggests that

508

the two O-H stretching modes (i.e. symmetric and antisymmetric modes) in this phase may be separated by as little as l0 meV. Again, there was little theoretical work in this area. In comparison with the gas phase of water, where the stretching modes are at 465.7 and 452.8 meV (345.8 and 330.5 meV for D20), showing high frequencies and small splitting. The presence of neighbouring molecules in the condensed phases increases the length of the O-H covalent bond in different ice structures and couples the intermolecular and intramolecular phonons. i

I

i

!

i

I

m

n

m

m

m

m

n

m

m

m

H20

,l

J J m

ice-lh

(

tff

i.-vll,

i

J

i

J

j

I

J

I

J

J

I

I

i

i

I

I

I

1 O0 200 300 400 energy transfer (meV)

500

Fig. 14. Neutron spectra from ice VIII measured using HET (Ei = 600 meV, at T ~ 10 K) on ISIS. The spectrum for ice Ih is also plotted for comparison.

5. SIMULATIONS OF NEUTRON SPECTRA OF ICE lh In recent years, the spectra of ices have been frequently simulated by means of LD and MD techniques. The early LD models were based on an assumed set of force constants for the lattice vibrations. By adjusting the force constants to fit

509 the observed spectrum, usually obtained from IR absorption and Raman scattering, the relevant inter- and intra-molecular forces were obtained. The early LD calculations for ice Ih were made by Faure and Chosson [60] and in their model, the orientational disorder of water molecules was ignored, treated the H20 units as point masses. The model provides crude dispersion relations and the DOS in the translational region. Shawyer and Dean [61 ] corrected this by using a super-lattice cell consisting of over more than 500 atoms (or ~166 molecules), a size which remains difficult to do even today because the large memory requirements of the dynamic matrix (see section 1.2). In their model, only the nearest neighbour interactions were taken into account by using a range of force constants (including intra- and inter-molecular forces) which are listed in Table 2. The aim of their calculations was to fit well for the IR absorption data [ 12]. Since there is only one optic peak, at--28 meV (or 220 cm'~), present in the IR spectrum, the task of matching the experimental data was not very difficult. The high energy shoulder at 36 meV (290 cm "1) in the spectrum was not reproduced, because it was assigned to either a combination or overtone. Later, Wong and Whalley [62] improved Shawyer and Dean's calculation by introducing point dipole-dipole interactions. Because there are four different local proton arrangements in the ice Ih structure (see Fig. 9), four different O-O stretching force constants were used and listed in Table 2. However, the additional contribution from the dipole moments are relatively small, compared with the strong O-O force constant, and the resulting spectra were very similar to the Shawyer and Dean's calculation in the translational region. Since INS spectra became available in late 60's, scientists have realised that the peak at 37 meV is also part of the one-phonon DOS. There was, therefore, a requirement to improve the existing LD model to reproduce both features in the INS spectrum. This was no simple task, as we demonstrate below. The normal picture of hydrogen bonding (Lennard-Jones term + charge-charge interaction) can not reproduce the two peaks in the DOS. One approach, by Renker [49], was to assume that the interactions along different crystal directions are different for a proton ordered structure of ice Ih. The O-O stretching force constant along the caxis was approximately 1.7 times larger than the stretching force constant in the basal plane (see Table 2). This approach was able to reproduce the two peaks observed in the INS spectrum of ice Ih that he measured at the Institut Max yon Laue-Langevin, France. However, our detailed study of this model found that the high energy peak at 37 meV in the calculated spectrum exists along the c-axis only, while the spectrum integrated over all the vibrational amplitudes in the basal plane only has one peak at 28 meV [55]. This is incompatible with neutron spectra measured along these two directions in a single crystal of ice Ih [ 14], both peaks are present in both orientations.

510 The lack of long-range charge-charge interactions in these models is often criticised since the results are far from realistic. In recent years, a number of LD simulations of ice Ih spectrum (measured by either IR, Raman or INS) were made using potential functions without a priori restriction on the number of force constants. Importantly, the long-range interactions between water molecules were included (e.g. by the implementation of Ewald-summation). Table 2. Comparison of force constants used in LD models: (unit: eV/Rad 2 for g/G and eV/A 2 for k/K)

Lattice dynamical models for ice lh Force Constant Type KO.H

Sh/D [61] 33.58 -1.17 0.74

Prask et al [48] 33.9' . . 1.64

Renker [49]

W/W [62]

Bosi et al

[**]

Li/Ross

[68]

38.7 36.1 . . 1.8 -1.2 1.83/1.56 1.5/2.1 1.8 1.1/2.1 1.78/1.49 gH-O-H 4.57 3.55 4.1 3.2 GH-O---H 0.148 ....0.50 0.31 0.78 GH-o---o 0.148 0.26 0.3 0.31/.4 0.62 0.61 GH---O---H -0.34 0.65 0.45 k and g are the internal stretching and bending force constants, K and G are the external stretching and bending force constants. The values for Ko--.n are the hydrogen bonding force constant. ** P. Bosi, R. Tubino and G. Zerbi, J. Chem. Phys. 59 (1973) 4578.

KH.H Ko... H

.,

_

The early LD simulations by Nielson and Rice used several classic pair-wise potentials, such as MCY and SL2 [63], to study a number of high-pressure ice structures: including proton-ordered ice II, VIII, XI and the proton disordered ice Ih and LDA [64]. Only the zone centre vibrations were calculated and the DOS were very noisy. Later, a similar LD calculation made by Marchi et al [66] improved the quality of the DOS, of ice Ih, by using a larger super-lattice cell with 128 molecules and SPC potential [65], hence more vibrational modes were included. The large cell size also allowed adequate disordering of the proton configurations. (MD calculations were also made by this group using the same structures and the SPC potential function.) Because the calculation was again limited to the BZ centre, the results were not a true DOS and could not be

511 compared with the available INS data. However, it was probably adequate for its intended use in the subsequent calculation of lR and Raman intensities. Later, Criado et al [67] calculated the DOS using the same SPC potential but with a much better BZ integration, 40x40x40 q-points in the first BZ. The results produced two peaks in the optic-phonon region at 28 and 33 meV, which are in reasonable agreement with those measured 28 and 37 meV. However, a protonordered structure of ice Ih (i.e. only 4 molecules in the unit cell with space group: Cmc2) was used in this calculation. In this ordered arrangement the local relative dipole orientations are type-B (mirror symmetry) along the c-axis and type-D (mirror symmetry) in the basal plane which have very different curvatures, i.e. force constants (see Fig. 9) associated with the two configurations. We believe it is this difference in force constants that gives rise to two optic modes if a proton ordered structure is used. Hence the result is very similar to the early work by Renker [49]. Indeed, when a proton-disordered structure is introduced the splitting becomes invisible, as we show later in section 6.

i

i'

I,,

I''''

I''''

1''''

i ....

i.1

i

-

g

(b)

I--

)

m z Q

~I

....

0

1 ....

10

I ....

20

l,,j

30

, I , , , , I

40

50

ENERGY (meV}

Fig. 15. Plot of both IINS measured (a) and LD calculated (b) spectra for ice Ih based on the Li and Ross (LR) model [68]. In an attempt to explain the splitting of the optic modes, the 'two strength hydrogen bond model' for ice Ih was proposed by Li and Ross [68]. In this model, a pair of hydrogen bonded water molecules have two different force constants according to their relative orientations. The values are 1.1 or 2.1

512 eV/~ 2. The ratio of the strengths of the forces, 1.9, is slightly greater than that of Renker's early models [49]. The significant difference is that the two force constants are chosen randomly according to the local proton configurations. In the models of both Renker and Criado et al, the stronger value was applied along the c-axis and the other in the basal plane. The LD calculations were carried out with a large proton-disordered structure and a high quality of BZ integration, 50x50x50 q-points. This approach seams to produce a result in good agreement with the experimentally measured spectra [55,68] (see Fig. 15). It indicates that local orientational variations of the true water potential could be considerably greater than that estimated from the classic water potentials. Moreover, if the existance of large orientational variations in the hydrogen bonding could be proved, a whole range of anomalous properties of water might be explained (see section 6.3). Hence it is of interest to explore other classic potentials in the literature, to determine whether any existing potentials can also reproduce the molecular optic band splitting in the INS spectrum for ice Ih by means of LD and MD simulation methods. 5.1. Lattice dynamic simulations by use of the classic potentials With the advance of computing techniques classic LD programs have become more and more sophisticated. The PHONON program, provided from Daresbury Laboratory [69], is one such excellent example. PHONON uses the quasi-harmonic approximation and has a wide range of two body potentials embodied in the code. In addition, angular three-body bending potentials, fourbody torsion potentials are also included. The program has been widely used for simulations of a variety of properties, such as dispersion curves, defects and surface phonons of crystalline and amorphous materials. Using the classic water potentials given in Table 3, an extensive survey was conducted by Dong [70]. The aim was to understand the fundamental reasons behind of the splitting of the molecular optic modes and to see which features could be produced by the existing classic potentials. His study, therefore, provides a much-needed understanding of the vibrational dynamics arising from these classic potentials. The simulation cell used for all the potentials consisted of 32 water molecules with hexagonal symmetry (ice Ih). The protons in the structure were disordered by use of a random walk program. The number of bonds with proton configurations, A/D (weak) and B/C (strong) were 28 and 40 respectively, compared with the ideal ratio of 1 92. The structure is not quite large enough to fully represent the observed proton disordering, because of limitations on the size of the super-lattice used. (Larger cells were also used and we are awaiting the results of these calculations.) However, the modest size sufficiently represents the mixture of different configurations, as we demonstrated in the

513

section 6.2 using a series of lattice sizes from 4 - 32 molecules, the result for the 32 molecules reproduces the measured spectrum very well. In this series of calculations, the lattice constants were initially set for c = 7.32 A and a = 4.50 A and the program relaxes both the molecular structure and lattice constants to minimise the total energy of the crystal. A particular local energy minimum was achieved for each of the potentials listed in Table 3. The additional intramolecular force constants, k (= 36 eV//~2) and g (= 2.6 eV/A 2) were also introduced for all the potentials used to stabilise the internal structure of the molecule, they are similar to those of other LD models in Table 2 and give reasonable intramolecular frequencies. As a consequence of the nonrigidity of the water molecule in this model, the dipole and quadruple moments are slightly larger than those given in other work. Here TIP4P provides 2.35 D as compared to 2.177 D [6]. This 5% increase provides a value closer to the experimental value of--2.5 D and similar increases were also found for quadruple values. They are the direct result of a slight increase in the O-H covalent bond length, since the water molecule is not rigid and the internal force is not infinite. After relaxation, the dynamical matrix is resolved and integrals of the phonon modes across the first BZ were made. The plots of the DOS for each of the potentials are shown in Fig. 16 and the curve at the bottom of the figure is the experimental data measured on TFXA [14] for comparison. Although the calculations included the intramolecular frequencies only the translational and librational regions are plotted here. The quality of the integration across the BZ is shown from the initial curve of the acoustic band from 0 to 5 meV, this is in good agreement with the Debye behaviour, i.e. g(m) -~m 2. Using a larger number of sampling points in the BZ can improve the roughness of the curves. Only 7x7x7 points were used in these calculations in order to make a like-to-like comparison with the MD results discussed below. In the translational region, the main features of the spectra for TIP2, TIP4P, Rowlinson [5] and BF potentials are similar to spectra calculated using the simple force constant model shown in Fig. 2. A single peak at less than 40 meV is predicted to dominate the molecular optic modes (the exact energy positions for the main features are listed in Table 3). For the MCY potential the molecular optic peak is shifted to considerably lower energy, 28.6 meV, indicating the hydrogen bonding is much soft than the reality. The left cut-off for the librational band is also much lower than other potentials used. Both of which were confirmed by MD simulations, see below. For TIP3P and SPC potentials, the optic peak is much broader than that generated by other potentials, which may indicates that the hydrogen bond stretching forces are spread by the charge-charge interactions (or the structures were not fully relaxed [70]).

514

TIPS2

m

-!

~

Oz '

LINSON

MCY

BF

.

,. ~

50 100 ENERGY TRANSFER (meV)

Fig. 16. Calculated spectra for ice Ih using the classic potentials listed in Table 3. The bottom curve is the measured spectrum for comparison.

~

o

0

~D

0 ~D

0

0

r

~

I ,_...~ ~

~,~

9

c'q

~

~r i0

9

',0 ~

'~ ~

9

9

.r

9

9

ce~

t~

O~

t.~

0

.~..

0~

oo ce~

r

"7

~

~

d"

t~

~'

~.0

t "~

" ~

t..:'~

~

~

~

t~

~t-- .~" ~

t"q it ~

O0

9

~

9

t~

~i o ~

t-- ~

~

'~t" O0

e~

o, o

~. ~

~r

~

r

~

o

~

515

516 The width of librational band and the left and right cut-off values reflect the orientational restrictions of the potentials. These restrictions usually involve the charge-charge interaction of the pair-wise potential (for details see Table 3). From the simulation results, we can conclude that the classic pair-wise potentials are unlikely to reproduce the double peak structure of the observed translational spectrum. This implies that the differences between the forces of the different proton configurations are too small and the potentials are, therefore, much too isotropic, as we illustrated in Table 1. The maximum difference among the four configurations is only -~15%, for most of the potentials investigated. A few potentials show larger differences, such as-~25% for ST2, LS and RS. As demonstrated below, section 6.2, optic mode splitting requires differences of at least 50%. In comparison, Renker's model has a difference between the two force constants of 70%. In the next section, we show a series of MD simulations performed using some of the classic pair-wise potentials and a few polarisable potentials. Because the polarisable potentials can take into account the many-body contributions, this represents a significant advance in the simulation of ice Ih spectra. An extensive trial of the shell-model was made using the PHONON program, to emulate polarisation effects. Unfortunately, interactions between the shells and the charges were too strong. Stable structures could only be generated with polarisation values [70]. It is for the same reason that smeared, or distributed, charges were used in other MD simulations [26,35]. This device avoids the so called "catastrophe" effect [26], which arises from closely interacting point-charges. 5.2. Molecular dynamic simulations The MD simulations of the vibrational dynamics have considerable advantages above LD. This is because a variety of effects, such as long-range charge interactions, anharmonicity and many-body terms are introduced naturally. The early MD simulations often used very small lattice cells, of less than 300 molecules, and simple potentials, such as SPC [66]. Here the motivation was comparison with the optical spectra for ice Ih. The advantage of using MD for the calculations of an optical spectrum is that the difficult selection rule analysis can be avoided for disordered systems and the intensity can be obtained numerically. Using flexible polarizable potentials, MD simulations provide the dipole moment- and polarisability- derivatives upon which the optical intensities depend. However, lack of information in the measured optic spectra in the translational mode region less than 40 meV (see Fig. 10), makes the MD calculations less useful. Hence, to compare INS spectrum for ice (i.e. the DOS) is the only viable alternative.

517

Recently, there were many attempts of MD simulations for the vibrational dynamics of ice. In these calculations more realistic, either non-rigid or polarizable, potentials were used. One such calculation was made by Itoh et al [72] using the KKY potential [9] which has three separate pair-wise terms: Voo(r), VoH(r), VHH(r) and an extra three-body term for H-O-H and H-O---H bending. These calculations produced the all the fundamental modes up to 450 meV (or 3622 cml). The resulting spectra show very similar features to results from the MCY and TIP4P potentials in the translational and librational regions (see Fig. 16 and 17). Using a polarised potential developed from the MCY potential (namely NCC potential [37]), Sciortino and Corongiu [73] have calculated the DOS for ice Ih using a cell o f - 4 0 0 molecules. The DOS reproduces the double peak feature in the translational region at 28 and 34 meV, but with very poor statistics and a small energy separation in comparison with the experimentally measured one. This indicates an incomplete summation over the BZ. Indeed, as we discussed above, 300 molecules would provide less than 20 q-points in the first BZ. Further work on a larger super-cell and improved energy resolution is needed to compare a better DOS with the INS spectra. The true features of the DOS will be invariant to the super-lattice cell size, the number of steps, step size and the BZ integration as we discussed in section 2. Reasonable sampling of the BZ requires a minimum of 512 molecules, with 20,000-40,000 steps of 0.1-1 fs, depending on whether a rigid or non-rigid potential is used. Bearing these factors in mind, Burnham et al [74] have studied a number of water potentials from the simple TIP4P, MCY to a sophisticated allatom polarisable potential [75] in an attempt to reproduce the INS spectrum. The MD program used for these simulations was MDCSPC4 developed at Daresbury Laboratory [71]. It uses an Ewald-sum to carry out the force and energy evaluations over the charges. The motion of the molecules and atoms is found using a 5 th order Gear corrector algorithm to integrate the equation of motion. Certain modifications were made in order to accommodate the range of the potentials used here. Rigid molecule potentials: The TIP4P and MCY potentials were chosen as two commonly used examples of rigid molecule potentials. (Other classic potentials listed in the Table 1 varied slightly under parameterisation). Both potentials are four-site models: they have 2H sites each holding a charge of +q/2, and an 'm' site (along the H-O-H bisector with O-m distance of 0.15 A for TIP4P for instance) with a charge o f - q but no mass and an O-site with no charge. Because all four sites are fixed on a rigid molecule no intramolecular bending or stretching frequencies are obtained. The H-O-H bond angle was set to 104.52 ~ and the O-H length at 0.9572 A. Both potentials were modelled with

518

a 520-molecule super-lattice cell. These simulations, as is the case for all the simulations presented here, were carried out at a temperature of 100 K.

-,,

30

9

,

1

8 6 4 2

g

0 2O

,. ~ ^~

---TIP4P

4

2

0 0.0

10.0

20.0

30.0

40.0

ENERGY TRANSFER (meV)

5O.0

0 50.0

70.0

90.0

110.0

laO.O

'

150.0

ENERGY TRANSFER (meV)

Fig. 17. A plot of the spectra calculated using MCY, KKY, TIP4P and SK potentials in the translational (on the left-hand-side) and librational (on the right) regions. The features shown in the figure are very similar to the results obtained from LD (see Fig. 16). Using the MCY potential at constant pressure and temperature the system became structurally unstable as described in ref. [74], even though the first nearest neighbour distance was preserved at about 2.9 A. A considerable distribution was found for the local tetrahedral symmetry. This behaviour is reasonable since a simple 6-12 potential has no preference for a tetrahedrally bonded structure. However, with a fixed cell volume the simulation became stable. Nearest neighbour molecules move within the energy minimum created by the pair-potential and the pair-wise additive electrostatic forces. At low temperatures, these molecules only sample the parabolic part of the potential

519 well and anharmonic effects should play only a minor role as discussed in section 2.2. The main features of the MD spectrum obtained from TIP4P are very similar to the LD result. The maximum intensity of the acoustic peak appears at 7.1 meV (57 cm ~) which agrees well with the experimental value of 7.05 meV [14]. The main feature in the optic mode region is still one band predicted at 33 meV (265 cml). The MCY potential is slightly more complex than the TIP4P and was fitted to ab-initio total energy calculations for the water dimer. The minimum position in energy occurred for an O-O separation about 3.12 A but the potential was rather shallow. This allows attractive forces to compress the lattice cell to give a O-O separation of~2.96 A which is significantly larger than the experimental value of 2.76 A, but it is comparable with other studies [33]. This shallow potential provides only a weak hydrogen bond force constant, -~0.28 eV/A 2, for the dimer (see Table 1). Although this value increases to -~1 eV/A 2 in simulated ice structures, it is still very weak compared with other potentials, having the value of--2 eV/A 2. As a consequence the translational band is shifted to low energy, ca. 29 meV. In addition, the librational band is very broad, 54-124 meV, c.f. observed width of 67-125 meV. For the TIP4P, the calculated optic peak is at 33 meV which is very closer to the measured high energy optic peak at 37 meV and the width of librational band (from 67-12 meV) agrees well with experimental data. However, again only one optic peak is predicted. Non-rigid water potential" Non-rigid water potentials were also used to simulate the ice spectrum, e.g. KKY potential. Although there are a number of non-rigid water potentials available in the literature, such as RSL [76], its complex forms for Voo(r), Vole(r) and Vm4(r) make it less adaptable. The KKY potential has been well studied recently for the water structure [7] and ice dynamics [72] and has the advantage of a relatively simple form for all three pair-wise terms:

gji(r): L (b i -1-bj )exp~ j b+a,+bj- r~j) +foD~j exp["2B0(r~J- r/J*)- 2exp~-B0(r/J - r/j*)) (38)

_[_CiCJ .~_e2zizj 6 An additional three-body term is also given by

V,,o,XO,,o,,) =

[;.(o,_,o,,- Oo)]-

(39)

520

1

k~-

(40)

exp [g~(ro~,)-r,)]+ l The r~ is an interatomic distance, the parameters, z, a, b, c, are related to the atomic species and D, B, r*f~ @zoz, rm and g~ are related to O-H pairs. The values of these parameters are given elsewhere [7]. The distances between the H atoms of adjacent water molecules are very different and could lead to considerable orientational variation of the force constants, much as one would expect to find in a polarisable model. However, the original parameters [7] failed to produce a stable structure for ice Ih [74]. A low value for 3~ was used in order to stabilise the structure and frequencies distribution function similar to earlier work were obtained [72].

!

r

!

|

7 6

Q r.~

5 4

CD

2

!

' 0

50

100

|

15o

260

260 36o'3~o'4oo

Energy transfer (meV) Fig. 18. Plot of the MD simulation result using the KKY potential. The spectrum provides the whole range of inter- and intra- molecular vibrational frequencies up to 400 meV. This MD simulation was carried out for a cell of 512 molecules with ~15000 (x 0.4 fs per step) steps, and, because of this larger cell, the quality of the calculated DOS is much improved over earlier calculations [72]. The most significant result is that the calculation contains both intermolecular and

521

intramolecular motions. The predicted bending and stretching modes, at 216 and 372-380 meV (see Fig. 18) are only in modest agreement with the experimental values, of 200 and 390-430 meV [15]. In the translational region, the acoustic peak at 12 meV is at much higher energy than that observed, 7.1 meV. The molecular optic mode is present as a single peak at 42 meV. This is again consistent with other pair-wise potentials but not with INS data. This potential overestimates the molecular bending and gives very high values for the librational modes, 97-164 meV, and c.f. INS 67-124 meV. Polarisable potentials: So far we have discussed MD simulations using different pair-wise additive potentials where the bulk properties were introduced through "effective" functions. A major drawback of these potentials is their lack of flexibility when treating water or ice of different densities. When parameters such as charges and positions are fixed the dynamic properties like dipole moment and many-body interactions are largely ignored. This is also tree for non-rigid pair-wise potentials, such as RSL and KKY. These are very significant terms, as we indicated in section 3. Here we shall illustrate one such MD simulation based on the SK potential. Other polarisable potentials were also used, such as those of DC [37] and Burnham et al [75]. Broadly, however, the results are more or less the same as for the SK potential shown here. The SK potential is a rigid-polarizable potential, which was developed by Sprik and Klein based on TIP4P potential [35]. Four fixed charge sites, containing Gaussian distributed electronic clouds, are arranged tetrahedrally around the m-site in order to give the angular variation of the dipole momentum. The magnitude of the charge can be varied, within the constraint of zero net molecular charge. The force acting on the m-site is to minimises the electrostatic and dipole energy of the molecules. In the case of no intermolecular interactions, the gas-phase dipole moment should be obtained and since the O-m separation is fixed there is no additional electronic contribution. The spectrum calculated using SK potentials is shown in Fig. 17. The features closely resemble the results from TIP4P in the energy transfer less than 20 meV. The principal difference is that the polarisation has increased the energy cut-off of the translational band from 36 meV for TIP4P to 47 meV for SK and broadened the peak considerably. This broadening phenomenon was our primary interest and was also observed from other polarisable potentials, such as DC [37] and Burnham et al [75]. This may imply that polarisation affects the strengths of hydrogen bonds differently for different proton configurations, hence the orientational variations of the potentials are greater than the pair-wise ones as we would expect.

522

5.3. Summary Although the MD simulations discussed above used a modest number of different potentials they cover an ample selection. The main aim of the simulations was to reproduce the INS spectrum for ice Ih, especially the split of the optic modes in the translational region. The potentials produced, broadly, the same spectral features, with some variation in band position and band width. Those relevant LD simulations, which used the same potentials as shown in section 5.1, confirm the reliability of the MD results. So far we have found no potential, tested by either LD or MD, which was capable of reproducing the measured spectrum for ice. However, there is clearly a tendency for polarisable potentials to broaden the optic features beyond that obtained from simple potentials of the TIP4P type which usually has a single, narrow optic peak. So far, NCC is the only potential function able to reproduce the optic mode splitting of 6 meV (the peak positions are at 28 and 34 meV) which is still much smaller than that measured 9 meV, despite the over estimation of the dipole moment (giving a value of 3.3 D) and the polarisation energy.

6. THE TWO STRENGTHS OF HYDROGEN BOND MODEL

In order to reproduce both features in the optic mode INS spectra of ice, Li and Ross proposed the 'two strengths of hydrogen bond' (LR) model [68]. They believed that the two molecular peaks are associated with different local dipoledipole configurations. They postulated that the relative intensities of the two optical bands are entirely dependent on the relative number of the two configurations. Moreover the different configurations are related to strong and weak H-bonds in the ice structure. For instance, in ice Ic (which has an identical INS spectrum to ice Ih), the protons are completely disordered. Hence, statistically, it has one C-configuration for every two D-configurations (see Fig. 9). Therefore, in ice Ic, there would be one weak H-bond for every two strong Hbonds. The situation is more or less the same for ice Ih, but here there are four configurations in the structure, which can be classed into two groups. Although we know type-C and type-D well (because they are shared by a number of phases, e.g. ice Ic, VII and VIII), the classification of the A- and B-types is difficult, because they only present in ice Ih and II. These strong and weak hydrogen bonds are randomly and isotropically distributed in the ice structure and their populations should correspond to the ratio of the integrated peak areas. The observed value, low-energy mode to highenergy mode, is about 1 92 which agrees well with the assumption of the protons in the structure is complete disordered. The difference between the force

523

constants used in these calculations, 2.1 and 1 eV/A2, is considerably larger than can be explained on the basis of electrostatic effects alone and the ratio, 2.1/1.1 = 1.9, compares poorly with that obtained from classic pair-wise potentials, 1.15-1.3 in Table 1. Unfortunately we have no means of estimating this ratio for polarisable potentials, because the many-body interactions in bulk ice are presently impossible to calculate. We could reasonably anticipate, however, that they are significantly larger than the maximum value achieved by the classic pairwise potentials based on the MD results in section 5.2. Our experimental and simulation results indicate that long-range interactions in ice are much weaker than we had imagined. Evidence from the INS spectrum of ice in small pores (radius -~10 - 30 A) is very similar to the bulk spectrum [54] and only a small difference is seen on the low energy side of the librational band at 68 meV. It is thought that this is due to the large numbers of water molecules close to the pore surfaces. Further, the similarity of the translational modes for the Ic, Ih and LDA ices indicates that the spectral features are not dependent on the long-range organisations of the lattice. Rather, the phonon frequencies are determined by nearest neighbour interactions. This implies that the role played by long-range interactions in water potentials may have much to do with the short intermolecular distances found in the structures produced by these potentials. If the difference between the forces produced by the distinct orientations of neighbours is only modest, then band splitting is unlikely to be observed, although band broadening may well be significant. Band splitting seems only to appear for the greatest anisotropy in the local force-field. To illustrate this effect we have varied the force constants among the different configurations. 6.1. The variation of the two force constants

A series of DOS calculations were performed with different ratios of K1 and K2 (other force constants were fixed at their Table 2 values). The results of the calculations are shown in Fig. 19; for a super-cell of 16 water molecules with four different ratios of the force constants: 1.8, 1.5, 1.3 and 1, while K2 fixed at 1.1 eV/~ 2. The features in the calculated g(co) for 1.8 shows excellent agreement with measured spectrum. The acoustic mode is well reproduced at 7 meV, it is both sharp and follows the Debye-model at low energies. It was the 20x20x20 q-points calculation that produced this ideal curvature and demonstrates the high quality of the BZ integration. The model also produces the correct positions for the optic modes and their triangular shapes. These triangular features are a direct result of the randomness of the strong and weak force constants used in the super-lattice. The phonons become localised giving modes that fill the gap between the two sharp features of earlier

524

model [49]. This proton disordering effect will be more clearly demonstrated in the next section 6.2.

(D)

(c)

E

I (B)

(A)

I

o.o

100.0

i

I

J

200.0

I

300.0

Energy transfer (cm-1) Fig. 19. Comparison of calculated PDOS for a super-cell of 16 water molecules with different ratio of the hydrogen bond force constants KI :K2 (the weak bond force constant K2 is fixed, having a value of 1.1 eV/flt2): the curve (A) is for the ratio of 1.8; (B) for 1.5; (C) for 1.3 and (D) for 1.0. The calculations also show that for a ratio of less than 1.5, the two optic bands begin to merge. Therefore, the INS data can only be reproduced when the ratio of the strong to weak force constants, among nearest neighbouring molecules, is greater than a critical value of 1.5. As we indicated in the earlier sections, the

525

classic pair potentials we have tested so far produce ratios are all less than the critical value. Indeed, if we assume that the orientational differences of the force constants come only from the charge terms, then, based on the classic dipole-dipole interactions it would be almost impossible to obtain a ratio greater than 1.5 for the pair-wise potentials. 6.2. Effects of Proton Disordering Ice Ih is a completely disordered proton system. In order to truly represent such a proton-disordered structure, an infinite lattice is required but is not realistically attainable. However, if the super-lattice is large enough to satisfy the following conditions it can be regarded as adequate. First, the averaged total dipole moment is near zero for the super-lattice used. Second, the calculated DOS has converged, where convergence implies that a further increase of the lattice cell size would not change the calculated results. The production of proton-disordered ice structures for LD and MD calculations is not trivial. If we assume that the protons in ice Ih or Ic structures are entirely disordered the strong and weak bond configurations have a ratio of 2 : 1. However, for small cells with periodic boundary conditions the ratio will vary topographically. A primary cell with 4 molecules (P63/mmc for oxygen network) has only two proton arrangements that obey the Bemal-Fowler ice rules [34]. These structural symmetries are Cmc2 and Cc. For 8 molecule orthorhombic cells, however, there are 17 proton arrangements and the strong- and weak- bonds can be mixed in an ordered manner [77]. For a lattice cell of 16 molecules, the possible proton arrangements are certainly over hundreds. To illustrate the size effect to the DOS, a series of LD calculations for a number of different size of ice Ih structures based on the same model were made as shown in Fig. 20. The curve (a) is the DOS for the Cmc2 structure. In this structure, all bonds in c-axis are strong and all bonds in the basal plane are weak. Therefore, the peak at 28 meV appears only in the integrated modes associated with vibrations in basal plane and the peak at 37 meV appears only on the c-axis of the hexagonal structure [78]. Although the peak positions are correct, the shape of the calculated spectrum does not agree with experimental data - the two peaks are very sharp and well separated from one another. This simulation result resembles the Renker's [49] and Criado et al's [67] results. When a large unit cell with 8 molecules was used, and with other properties remaining fixed, the improvement in the calculated result is shown in the curve (b). The strong and weak bonds are present in both crystal orientations. Another important improvement is that more phonon frequencies are found in the gap between the peaks, these will gradually build up into the shape seen in the measured spectrum. When the lattice cell is increased to 32 molecules, the

526

resulting spectrum, shown in the curve (d), begins to look very much like the measured spectrum, because the proton-configurations can be reasonably well mixed in all directions.

o

I'

or) iii

oo ii O >oo z iii

E)

d

|

I

0

. . . .

I

. . . .

1

t

~

~

,

!

. . . .

I,,,,1

10 20 30 40 ENERGY TRANSFER ( m e V )

50

Fig. 20. A plot shows a series of LD results using the different sizes of the ice Ih lattice to represents the proton disordering: (a) for a lattice cell with 4 molecules; (b) for 8 molecules; (c) for 16 molecules and (d) for 32 molecules. From this series of calculations, we have demonstrated that proton disordering in the lattice can generate the common triangular shapes seen for the two optic

527

peaks. The integrated intensities of the two peaks are directly proportional to the ratio of the numbers of the strong and weak bonds in the structure, in agreement with observation.

6.3. The anomalies of water and ice. The above conclusions were based on calculations and measurements taken from a variety of ices [68]. They offer the prospect of defining a potential for the water molecule that not only satisfactorily reproduces structural data but also generates an acceptable DOS. At present the LR model appears to offer considerable promise in this direction. The local structure of water is often considered to be ice-like and a good model for ice would be an obvious candidate for the structure and dynamics of water. Indeed there are some indications that the two peaks are present in the INS of liquid water, but shifted to lower energies, 24 meV and 32 meV [79]. Therefore, it is appropriate to briefly review the consequences that this model would have for the liquid state of water. Of course this process is by its very nature speculative but it does draw out the intriguing number of unusual properties of water that can be addressed through the Li/Ross model. Melting and boiling temperatures: In the LR model the strong H-bonds have a slightly greater bonding energy (i.e. more negative) than the weak ones. The vibrational frequency for the lower energy optic peak (the weak bonds) is about 24 meV, the thermal energy is very close the melting-point of ice 0~ (considering 300 K = 25 meV). At this temperature, the ice structure can not be sustained and weak bonds would be broken. The continuous ice structure would degenerate into large water clusters mainly connected by the strong H-bonds. The higher the temperature, the smaller the water clusters would become. When the temperature is high enough to break even the strong bonds, water molecules can evaporate. The vibrational frequency of the high-energy peak is 32 meV (-~ 380 K) which is almost equivalent to the boiling temperature of water, 100~ High heat capacity: Previously, the number of broken hydrogen bonds in liquid water was used to explain its high heat capacity. We believe that the progressive altering of the ratio of the strong to weak bonds may account for this high heat capacity. The strong and weak bonds act like an energy reservoir. Increasing the temperature would convert the lower energy states, of the strong bonds, to the higher energy states, of the weak bonds. High surface tension: Water molecules on the surface can readily orient themselves into the lowest energy configuration, which is the strongly bonded state. This phenomenon has been observed in the INS spectrum of vapour deposited ice [51]. Porosity in these ices is very high and a large number of water molecules are on surfaces. The INS spectrum shows a single, dominant,

528

peak at 37 meV at 10 K which is associated with the strong bonds. The measurements for water on surfaces of porous media such as silica gels and Vycor [54] show very similar spectra, again one peak at higher energy, 37 meV. This all points to a high population of strong bonds at the surface of water and this would inevitably lead to high surface tension. Polymorphism of ice: The model also casts a new light on our understanding of polymorphism in ice. The complex phase diagram of water has been explained in terms of the openness of the H-bond structure. However, tetrahedral structures are not restricted to ice (e.g. Si, Ge and diamond). In the LR model [68] ice Ih is under significant local stress, arising from the mixture of strong and weak bonds and their respective lengths. When external pressure is applied, the different configurations (or bond types) respond differently. In this picture, new phases appear as they are best able to relax these internal stresses. A broad range of different structures would be anticipated as well as metastable behaviour. Taking ice II as an example, the bonds in the hexagonal ring are all strong Hbonds, while the bonds between different hexagonal rings are all weak bonds. A rotation of the individual members of the rings allows energy relaxation as they switch from configuration type-C to configuration -D [58]. In ice VIII, because of the strong repulsion between the two interpenetrating sub-lattices, all the hydrogen bonds are considerably stretched with a O-O distance of 2.98 A. Weak bonds would require less energy to stretch and be energetically more stable than a structure of strong bonds. We believe that the proton-disordered forms of ice are frustrated systems. These "equilibrium" systems result from suitable mixtures of different bond types with different bond lengths. This idea is supported by the fact that many phases of ice can only be obtained by following specific paths in the T-P diagram. Stress energy is also able to account for the small energy differences between ice Ih and Ic. Quite simply Ih has an extra lattice parameter, c(-axis), which can be optimised in response to the surrounding stress and the total free energy of ice Ih is, therefore, lower than that of ice Ic. The total free energy of ice Ih can be reduced further when proton ordering increases. Such as in ice XI, where the c/a ratio decreases from that of ice Ih, 1.628, to 1.617, whilst the value in ice Ic is 1.632 [80]. Despite the large force constant difference between the strong and the weak bonds, the total energy difference between structures involving the two bonds may yet be rather small. There is, therefore, only slight stability to be gained from adopting an ordered structure for ice Ih(c) and much may be lost to entropy. Furthermore, this mixture of different configurations with their short and long bonds may also cause the positional disorder observed for oxygen

529

atoms [81] and why protons are sometimes found away from the O---O axis [82]. Internal stress could play a major role in the stability of ice structures.

7. DISCUSSION In this article, we have presented a series of LD and MD simulations for ice Ih using a variety of water potentials and the results were compared with INS measured DOS. Neutron measurements were shown to provide unique information on the fundamental intramolecular and intermolecular modes, some of which cannot be obtained from the standard IR and Raman techniques. A full knowledge of the intermolecular vibrations as modulated by the molecule's environment in the lattice systems is necessary for a complete analysis of the dynamics of these ice structures. Equipped with the precise knowledge of the structural information obtained by the diffraction measurements [81,82], one can model the system rigorously with suitable force fields or potential functions. The extensive simulation results show that classic pair-wise potentials were unsuccessful in reproducing the measured DOS for ice Ih. From the simulations, we conclude that two hydrogen bonding force constants are a basic requirement for reproducing the measured spectrum. If a water-water potential generates sufficiently large force constant differences for the different proton configurations (or the different relative dipole-dipole orientations in water or ice), it should produce the same effect as seen in the LR model. The anisotropic properties of the classic potentials are a result of charge interaction and this anisotropy should increase in the polarisable potentials and hence they produce a broad optic peak. This broad peak indicates that the orientational variation of the potential function has been increased considerably but it may still be less than the critical value of 1.5 as we indicated in the section 6.1. One would, therefore, expect that a better polarisable potential would, eventually, be able to reproduce the split optic peaks in the measured INS spectrum.

ACKNOWLEDGEMENTS The authors would like to thank the Engineering and Physical Science Research Council (UK) for financial support and the Rutherford-Appleton Laboratory for the use of neutron facilities. We would also like to think Mr. C.J. Burnham and S.L. Dong for providing a number of graphics, which they produced as part of their Ph.D. studies.

530

REFERENCES

1. R.A. Home (ed.), "Water and Aqueous Solutions", John Wiley & Sons, New York, 1972. 2. H.S. Frank, Water: A Comprehensive Treatise, F. Franks (ed.), Plenum Press, New York - London, 1972. 3. D. Eisenberg and W. Kauzmann (eds.), The Structure and Properties of Water, Oxford Univeristy Press, Oxford, 1965. 4. O. Matsuoka, E. Clementi and M. Yoshimine, J. Chem. Phys. 57(1976) 1351. 5. F.H. Stillinger and A. Rahman, J. Chem. Phys. 60 (1974) 1545. 6. W.L. Jorensen, J. Chandrasekhar, J.D. Maura, R. W. Impey and M.L. Klein, J. Chem. Phys. 79 (1983) 926. 7. N. Kumagai, K.K. Kawamura and T. Yokokawa, Mol. Sim., 12 (1994) 177. 8. G.W. Robinson, S.B. Zhu, S. Singh and M.W. Evans (eds.), Water in Biology, Chemistry and Physics, World Sciencitic, 1996. 9. P.V. Hobbs (ed.), Ice Physics, Clarendon Press, Oxford, 1974. 10. G.E. Walrafen, M.R. Fisher, M.S. Hokmabadi and W.H. Yang, J. Chem. Phys. 85 (1986) 6970 and 94 (1990) 2237. 11. Y. Marechal, J. Chem. Phys. 95 (1991) 5565. 12. J.E. Bertie and E. Whalley, J. Chem. Phys. 40 (1968) 1646. 13. B. Minceva-Sukarova, W.F. Sherman and G.R. Wilkinson, J. Phys. C 17 (1984) 5833. 14. J.C. Li, J.D. Londono, D.K. Ross, J.L. Finney, J. Tomkinson, and W.F. Sherman, J. Chem. Phys. 94 (1991) 6770. 15. J.C. Li, J.D. Londono, D.K. Ross, J.L. Finney, S.M. Bennington and A.D. Tayor, J. Phys.: Condens. Matter 4 (1992) 2109. 16. G.L. Squires (ed.), Introduction to thermal neutron Scattering, Dover Publications Inc., New York, 1996. 17. L. van Hove. Phys Rev 95 (1954) 249. 18. J. Tomkinson, M. Warner and A.D. Taylor, Mol. Phys. 51 (1984) 381. 19. H. Conroy J. Chem. Phys. 59 (1973) 3992. 20. G. Gilat (ed.), Methods in Computational Physics, Vol 15, Acadaemic Press, New York, 1976, p317. 21. A. Kolesnikov and J.C. Li, Physica B 234-236 (1997) 34. 22. J.C. Li, D.K. Ross, M.H.B. Hayes, W.F. Sherman and M. Adams, NATO ARW on Hydrogen bond Networks, M.C. Bellissent-Funel and J.C. Dore (eds.), Kluwer Academic Publishers, 1994, p139. 23. E.O. Brigham (ed.), The fast Fourier Transferm, Prentice-Hall, Inc. New York, 1974. 24. R.K. Pathria (ed.), Statistical Mechanics, Pergamon Press, Oxford, 1985.

531

25. M.P. Allen and D.J. Tildesley (eds.), Computer simulations of liquids, Oxford Science Publications, Clarendon Press, Oxford, 1987. 26. C.J. Burnham, Ph.D. thesis, University of Salford, Manchester, UK, 1998. 27. P.H. Berens, D.H.J. Mackay, G.M. White and K.R. Wilson, J. Chem. Phys. 79 (1983) 2375. 28. W.C. Rontgen, Ann. Phys. U. Chem. 45 (1892) 91. 29. J.A. Pople, Proc. Roy. Soc. A205 (1957) 163. 30. H.S. Frank and W.Y. Wen, Dies. Faraday Soc. 24 (1957) 133. 31. J.C. Li, J. Phys. Chem. 101 (1997) 6237. 32. J.L. Finney, J.E. Quinn and J.O. Baum, Water Science Review, F. Franks (ed.), Cambridge University Press, 1985, Vol 1, p93. 33. M.D. Morse and S.A. Rice, J. Chem. Phys. 76 (1982) 650. 34. J.D. Bernal and R.H. Fowler, J. Chem. Phys. 1 (1933) 515. 35. M. Sprik and M. L. Klein, J. Chem. Phys. 89 (1988) 7556. 36. L.X. Dang and T-M. Chang, J. Chem. Phys. 106 (1997) 8149. 37. U. Niesar, G. Corongiu, E. Clementi, G.R. Kneller and D.K. Bhattacharya, J. Phys. Chem. 94 (1990) 7949. 38. C. Millot and A. Stone, Molecular Physics, 77 (1992) 439. 39. K. Liu, J.D. Cruzan and R.J. Saykally, Science, 96 (1992) 1832. 40. K. Liu, M.G. Brown, R.J. Saykally, R.J., Gregory and D.C. Clary, Nature, 381 (1996) 501. 41. S.S. Xantheas, J.Chem. Phys. 100 (1994) 7523 and 102 (1995) 4505. 42. J.K. Gregory and D.C. Clary, J. Phys. Chem., 100 (1996) 18014. 43. P. Bemes, J.L. Finney, J.D. Nicholas and J.E. Quinn, Nature, 282 (1979) 77. 44. D.N. Bemardo, Y. Ding, K. Krogh-Jespersen and R.M. Levey, J. Chem. Phys. 98 (1994) 4180. 45. J.C.Dore, Water Science Review, F. Franks (ed.), Cambridge University Press, 1985, Vol 1, pl. 46. P.A. Egelstaff, J.A. Polo, J.H. Root and L.J. Hahn, Phys. Rev. Lett. 47 (1981) 1733. 47. A.K. Soper and M.G. Phillips, Chem. Phys. 107 (1986) 47. 48. H. Prask, H. Boutin and S. Yip, J. Chem. Phys. 48, (1968) 3367. 49. B. Renker, in Physics and Chemistry of Ice, E.Whalley, S.J. Jones, and L.W. Gold (eds.), University of Toronto Press, Toronto, 1973, p82. 50. B.C. Boland (ed.), ISIS Facilities, Rutherford Appleton Laboratory, 1990. 51. A. Kolesnikov, J.C. Li, S.L. Dong, I. Balley, W. Hahn. R. Eccelston and S. Parker, Phys. Rev. Lett. 79 (1997) 1869. 52. J.C. Li and M. Adams, EuroPhys. Lett. 34 (1996) 675. 53. S.F. Parker, K.E. Horton and J. Tomkinson, (eds.) The TFXA User-Guide, Rutherford Appleton Laboratory: Technical Report, 1995.

532

54. J.C. Li, D.K. Ross and M.H.B. Hayes, J. Mol. Struc. 322 (1994) 131. 55. J.C. Li, J. Chem, Phys. 105 (1996) 6733. 56. A.I. Kolesnikov, J.C. Li, D.K. Ross, O.I. Barkalov, V.V. Sinitsyn, E.G. Ponyatovsky, J. Phys. Lett. A 168 (1992) 308. 57. S.L. Dong, A. Kolesnikov and J.C. Li, J. Phys. Chem. 101 (1997) 6087. 58. S.L. Dong, Y. Wang, A.I. Kolesnikov and J.C. Li, J. Chem. Phys. 109 (1998) 235. 59. A. Kolesnikov and J.C. Li, Submitted to Phys. Rev. B. 60. P. Faure and A. Chosson, J. Glaciology 21 (1978) 85. 61. R.E. Shawyer and P. Dean, J. Phys. C5 (1972) 1028. 62. P.T.T. Wong and E. Whalley, J. Chem. Phys. 65 (1976) 829. 63. W.L. Jorensen, J. Am. Chem. Soc. 103 (1981) 335. 64. G. Nielson and S.A. Rice, J. Chem. Phys. 80 (1984) 4456 and G. Nielson, R.M. Townsend and S.A. Rice, J. Chem. Phys. 81 (1984) 5288. 65. H.J.C. Berendsen, J.P.M. Postma, W.F. van Gusteren and J. Hermans, Intermole. Forces, B Pullman (ed.), Reidel, Dordrecht, 1981, p331 and H.J. Berendsen, J.R. Grigera and T. Straatsma, J. Phys. Chem. 87 (1987) 6269. 66. M. Marchi, J.S. Tse and L. Klein, J. Chem. Phys. 85 (1986) 2414. 67. A. Criado, F.J. Bermejo, M. Carcia-Hemandez and J.L Martinez, Phys. Rev E. 47 (1992) 3516. 68. J.C. Li and D.K. Ross, Nature 365 (1993) 327. 69. M. Leslie, Daresbury Laboratory, Warrington, England, WA4 4AD. 70. S.L. Dong, Ph.D. thesis, UMIST, Manchester, UK, 1998. 71. W. Smith (ed.), MDCSPC4, Technical Memorandum DL/SCI/TM84T, Daresbury Laboratory, 1991. 72. H. Itoh, K. Kawamura, T. Hondoh and S. Mae, Physica B 219 (1996) 469. 73. F. Sciortino and G. Corongiu, J. Chem. Phys. 98 (1993) 5694. 74. C.J. Burnham, J.C. Li and M. Leslie, J. Phys. Chem. B, 101 (1997) 6192. 75. C.J. Burnham, J.C. Li, S.S. Xantheas and M. Leslie, J. Chem. Phys. (1998), 9 accepted. 76. A. Rahman, F.H. Stilinger and H. Lemberg, J. Chem. Phys. 63 (1975) 5223. 77. C.M.B. Line and R.W. Whitworth, J. Chem. Phys. 104 (1996) 10008. 78. J.C. Li and M. Leslie J. Phys Chem. 101 (1997) 6304. 79. S.H. Chen, Proc. of the NATO Adv. Study Inst. on Hydrogen-Bonded Liquids, J.C. Dore and J. Teixeira (eds.), 1989, p289. 80. R. Howe and R.W. Whitworth, J. Phys. Chem. Solids 50 (1989) 963. 81. W.F. Kuhs, J.L. Finney, C.Vettier and D.V. Bliss, J. Chem. Phys. 81 (1984) 3612. 82. W.F. Kuhs and M.S. Lehmann, Strcutr of ice Ih, in Water Science Review, F. Franks (ed.), Cambridge University Press, 1985, Vol 2, p 1.