Disordered Atom Molecular Potential for Water Parameterized against

Article pubs.acs.org/JPCB

Disordered Atom Molecular Potential for Water Parameterized against Neutron Diffraction Data. Application to the Structure of Ice Ih Alan K. Soper* ISIS Facility, Rutherford Appleton Laboratory, Harwell Science and Innovation Campus, Didcot, Oxfordshire OX11 0QX, U.K. ABSTRACT: A disordered atom molecular potential (DAMP) for water is described that accurately accounts for the observed neutron interference differential scattering cross sections for light water, heavy water, and two different mixtures of these liquids (x = 0.5 and x = 0.64, where x is the mole fraction of light water in the mixtures) at T = 283 K. This potential, when used in a NVT Monte Carlo computer simulation, produces an intermolecular pressure of ∼0 kbar and a configurational energy of approximately −50 kJ/mol, close to the values found in the ambient liquid at this temperature. The same potential is used as the reference potential in an empirical potential structure refinement of ice diffraction data at T = 258 K measured at the same time as the water data and under the same conditions. Particularly intriguing is the finding that the O···O−H angle in ice, which would be 0° for a linear hydrogen bond, is actually more disordered in ice than in the liquid. A rationalization of these findings is presented. It remains to be seen whether this potential has any value other than simply as a description of the ambient liquid structure.

■

INTRODUCTION Of all the liquids ever studied, water must have, rightly or wrongly, pride of place for being the most researched liquid of all time. This interest arises from a number of factors, only one of which is the fact that water is essential to life on Earth and possibly to life in the universe.1 Arguably as important is the fact that water is extremely hard to describe in the notation that works for most other small molecule liquids:2−4 on heating from 273 K it has a density maximum and compressibility minimum; the diffusion constant initially increases when low temperature water is pressurized; its most common solid form, ice Ih, is less dense than the liquid at the same temperature; when supercooled, the specific heat increases markedly, suggesting a possible hiatus at a lower, unreachable, temperature. Hence, water presents an intriguing set of properties that have produced significant literature documenting attempts to understand them over many years. Essential to understanding these properties are the molecular level forces between and within water molecules. Starting from the initial Bernal and Fowler attempt to describe these forces in terms of the repulsive and van der Waals forces common to all models of liquids plus a set of effective Coulomb charges placed on the atoms within the molecule,5 a number of such potentials have evolved,6−10 most of which are still in wide use in computer simulations of water and aqueous systems. Some versions of these potentials incorporate, via path integral molecular dynamics, the correct quantum nature of the proton in water,11−14 thereby capturing quite accurately the measured intramolecular as well as intermolecular interactions. The parameters in these potentials were mostly chosen by fitting © XXXX American Chemical Society

them to known water and ice properties, such as thermodynamic and structural quantities. However, dissatisfied with the ad hoc nature of these parameters, a number of other approaches were tried based on, for example, density functional theory combined with ab initio molecular dynamics.15 Another lucrative thread is the use of effective potentials that are derived from quantum mechanical calculations of the interactions of clusters of water molecules,16−18 then applied to the liquid and solid states of water: these potentials are not fit to thermodynamic and structural properties as such and often include many-body interactions so that in principle they can be applied over a wide range of state conditions. These later methods apparently capture water properties quite accurately. One of the tests of the accuracy of any particular choice of molecular force field for water is its ability to reproduce the measured structure of water as manifest by the O−O, O−H, and H−H radial distribution functions. Working within the approximation that these functions are only weakly dependent on isotope, it is in principle possible to extract them from neutron diffraction experiments on heavy water, D2O, light water, H2O, and mixtures of these liquids or else in combination with X-ray diffraction experiments. This is because the hydrogen scattering length for neutrons changes with isotope, while the X-ray scattering from water is dominated by Special Issue: Branka M. Ladanyi Festschrift Received: September 30, 2014 Revised: November 7, 2014

A

dx.doi.org/10.1021/jp509909w | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

from an alloy of titanium and zirconium that has almost zero coherent scattering for neutrons. The samples were 1 mm thick, and the front and back wall thicknesses were each 1 mm. The cans were mounted on a sample changer with a water bath temperature control that was able to cycle between the various samples at each temperature and also to control the temperature within 1 K. Sufficient space above each sample was provided to allow the liquid to expand into when it froze. The advantage of freezing the sample in situ is that it means the amount of material in the neutron beam is known, based only on the density of ice and the thickness of the container. An alternative approach based on making ice and powdering it outside the diffractometer, then placing it in the beam, is unsatisfactory because it introduces the uncertainty of how much material is actually in the neutron beam due to packing effects, making absolute normalization of the scattering more problematic. A disadvantage of freezing the sample in situ is that the ice sample will likely have significant preferred orientation, making the obtaining of truly powder averaged data difficult. To partially circumvent this latter problem, the sample was repeatedly thawed and refrozen, after checking at each cycle that all evidence of crystallinity had disappeared before starting the next freezing cycle. The result was two sets of scattering data, one in the liquid state at 283 K and the other in the crystalline state at 258 K. Four mixtures of heavy and light water were studied, namely, pure D2O, pure H2O, and two mixtures of these two liquids, with x = 0.5 (here called ”HDO”) and x = 0.64 (here called ”Null”), where x is the mole fraction of light water in the mixture. The latter case corresponds to the situation where the hydrogen atoms have zero coherent scattering because the neutron scattering length for a proton is of opposite sign to that for a deuteron. Hence, within the approximation that heavy water and light water have the same structure, the ”Null” sample is in principle sensitive to only the O−O pair distribution function. All the scattering data were corrected for detector deadtime, background, attenuation, multiple scattering, and container scattering and put on an absolute scale of differential scattering cross section (units of barns sr−1 atom−1; 1 barn = 10−24 cm2) by comparison with scattering from a 3 mm thick vanadium plate, which, like the samples, spanned the beam. The single atom and inelastic scattering were subtracted from the data by the iterative procedure described in ref 26. For the ice data, although the instrumental resolution on SANDALS is not good enough to perform a full Bragg peak crystal refinement, the observable Bragg peaks in each sample of ice were successfully indexed against the known structure of ice at this temperature.31,32 Setting Up a Reference Potential for Water. The majority of simple point charge water potentials use a LennardJones potential to describe the short-range repulsive and longer range dispersion (attractive) forces between molecules in addition, of course, to the effective Coulomb charges placed on particular sites. In the early days of computer simulation, there was considerable dispute between the height and position of the first peak in the O−O radial distribution function gOO(r) as obtained by simulation compared to that obtained in X-ray scattering experiments.33 The experimentalists believed the simulated peak was too high, while the simulators believed the experimental peak was too low. In fact, the much more recent X-ray work27 shows that this peak clearly is not as high and occurs at a larger radius, ∼2.8 Å, than the early simulations

the O−O correlation. In practice, however, the outcomes of these experiments have apparently differed quite widely,19−25 which happens primarily because of the large incoherent and inelastic scattering of neutrons by protons. Unfortunately there is currently no exact way that these effects can be removed from the data, leading to discrepancies in the way these effects are dealt with prior to structural analysis. Fortunately a recently introduced iterative procedure for subtracting inelasticity effects26 appears to have removed some of this uncertainty: comparison of the O−O radial distribution function obtained in that work is closely similar to that reported independently from X-ray diffraction.27 In addition, the methods used to interpret the scattering data vary quite widely, from direct Fourier transform of the data, as in refs 21 and 25, to reverse Monte Carlo (RMC), as in refs 22 and 24, and empirical potential structure refinement (EPSR), as in refs 23 and 26. The EPSR approach attempts to build a degree of reality into the modeling process: molecules are defined by realistic harmonic potentials, and the starting potential (called the “reference” potential) of an EPSR simulation for water can be one of the site−site potentials that have been described above. In fact there would be nothing in principle to stop any of the previously described potentials and methods being used as the reference potential in an EPSR simulation. Concepts like hydrogen bonding, van der Waals forces, atomic overlap forces, and so on are therefore built into the simulation at the outset, and it does not rely solely on the data to supply this information except with regard to any features of the reference potential that are clearly incompatible with the data. The data generate a perturbation to that reference potential (called the “empirical” potential because it is derived from the data) only to the extent necessary. Of course if there is a feature of the reference potential that prevents the simulation from fully reproducing the data, even with the empirical potential present, this means there is either a fault with the reference potential or a fault with the data. The success of this approach is confirmed because it will be noted that reference interatomic potentials, or their equivalents, are now appearing in RMC analyses.28−30 To date, not one of the simple point charge site−site potentials is able to describe existing neutron and X-ray scattering data from water at a level of accuracy compatible with the measuring accuracy, even with the uncertainties associated with data artifacts. This comment is less true of the more recent potentials that include molecular flexibility, quantum effects, and three-body forces.13,18 Even so, the question arises whether there is any site−site reference potential for water that accurately describes all the scattering data for a range of heavy and light water mixtures without the need for an accompanying empirical potential. The rest of this paper demonstrates that there is at least one such potential, and there may be a range of such potentials, depending on exactly how they are defined. In addition, by use of some new neutron scattering data from mixtures of heavy and light water and ice above and below the freezing point, the hydrogen bond structure in ice just below the melting point is compared with the same structure in the liquid just above the freezing point.

■

METHODS Neutron Scattering Experiment. The experiment described here was performed on the small angle neutron diffractometer for amorphous and liquid samples (SANDALS) at the ISIS pulsed neutron facility, U.K. Mixtures of heavy and light water were contained in standard flat plate cans made B



Article

suggested. In fact the first ab initio simulations of water15 showed that the first peak in this function could indeed be quite low in height and at a larger radius; it all depended on the softness of the repulsive core. Recent work34 fully supports these assertions. Therefore, for the present purpose, the use of a Lennard-Jones core potential was abandoned in favor of a softer, Gaussian potential placed at r = 0, which however lacks an attractive dispersion force. The steepness of the repulsive core was then controlled via the amplitude and width of this potential. Using a finite core potential introduces a possibility for errors to occur when running a Monte Carlo simulation, since the Coulomb potential for point charges diverges at r = 0, and this can cause atomic overlap between sites of opposite charge. To prevent this from happening, the point charges conventionally used are replaced here by a charge cloud distributed uniformly inside a sphere of specified radius. When the two spheres do not overlap, the interaction potential with charges qa, qb follows exactly Coulomb’s law, qq (Coul) (r ) = a b Uab 4π ϵ0r (1)

I2(r , b , x) =

1 2 ⎡ 2 2⎤ qaqb ⎢ 3 ra − 5 rb − r ⎥ = ⎥ 4π ϵ0 ⎢⎢ 2ra 3 ⎥⎦ ⎣

(

)

(2)

N

which can be derived by straightforward application of the electrostatic properties of uniform charged spheres. A more complicated but still analytic expression applies when (ra − rb) < r < (ra + rb): qq (Coul) Uab (r ) = a b [I1(r , ra , rb , x) − I2(r , rb , x)]−ra −rb r 4π ϵ0 (3)

U (intra) =

6 ⎤ 3 ⎡⎢ c xi⎥ 3 3 ∑ i 8a b r ⎢⎣ i = 1 ⎥⎦

∑

c1 = rb2(3a 2 − r 2)

c3 =

r(3(a 2 + b2) − r 2) 3

c4 = −

3r 5

c6 =

1 6

(6)

mα(i)mβ(i) mα(i) + mβ(i)

(7)

which replaces the need to specify a distinct Debye−Waller factor for each term in eq 6. C is a constant that controls the width of the intramolecular correlation functions and is normally set to 65kT/Å to give a distribution of distances in these correlation functions that is consistent with what is observed experimentally. In addition to these main features of the proposed reference potential, some further, mostly small amplitude, exponential and Gaussian terms are added to the O−O, O−H, and H−H intermolecular potential functions to improve the overall agreement between simulation and experiment. Each exponential and Gaussian contribution to the pair interaction potential is characterized by a position p, a width w, and an amplitude A, with a set of values for each of the O−O, O−H, and H−H site−site pair potentials. The exponential potential is given by

(3(a 2 − r 2) + b2) 4

c5 =

dα(i)β(i)

(4)

wα(i)β(i) = C

3b2(a 2 − r 2) 2.0

(dα(i)β(i) − rα(i)β(i))2

where the first summation proceeds over all the N molecules in the simulation box and the second proceeds over all the distinct pairs of atoms, α(i),β(i), within each molecule, for which interatomic distances dα(i)β(i) are specified. rα(i)β(i) is the actual distance between these atoms found in each instance of the molecule. The weighting function for each term in the summation is given by

with

c2 =

∑

wα(i)β(i)

i = 1 α(i), β(i) > α(i)

where I1(r , a , b , x) =

(5)

Hence the present approach is analogous, but different in detail, to the Gaussian charge model already used with considerable success.34 Instead of performing the full Ewald summation to correct for long-range effects, the local molecular field (LMF) theory35 is invoked to truncate the Coulomb force field with a complementary error function, using a field width of σ = 3.5 Å: Tc(r) = erfc(r/σ). No further mean field correction was made for the slowly varying part of the residual potential. This approach has been shown to give excellent results in the case of modeling water with Coulomb force fields.36 A characteristic feature of EPSR simulations with molecules is that these simulations must capture correctly the quantum mechanical zero point disorder of the real molecule. This is because the scattering data at large wave vector (Q) decay more rapidly than occurs in classical force models. While the ideal here would be to perform a full path integral quantum simulation,13 a pragmatic approach is to allow the atoms within a molecule to distribute according to a harmonic potential, but treat the molecule as rigid for most simulation moves, with only the occasional “shake” so that each molecule adopts a (slightly) different conformation from time to time as the simulation proceeds, which however remains, on average, consistent with the specified molecular structure. The intramolecular potential is therefore defined by the expression

but when the spheres overlap, the interaction potential remains finite, becoming quadratic with r near r = 0. For example, the potential energy of two overlapping uniform charged spheres of radii ra, rb with ra ≥ rb is given, for r ≤ ra − rb, by (Coul) (r ) Uab

3 ⎡ 2 x3 ⎤ bx+ ⎥ 3⎢ 3⎦ 4rb ⎣

and C



Article

⎛ (p (1) − r ) ⎞ (expon) (r ) = A ab(1)⎜⎜ ab Uab − 1⎟⎟ ⎝ wab(1) ⎠ ⎡ (p (1) − r ) ⎤ ⎥ exp⎢ − ab wab(1) ⎥⎦ ⎢⎣

intermolecular site−site potentials. The values of Coulomb charges and the positions, widths, and amplitudes of the exponential and Gaussian terms were adjusted by trial and error to give the best possible fit to the scattering data while also constraining the pressure to near 0 kbar and the configurational energy to near −50 kJ/mol. To summarize this section, the total intermolecular potential between atom types a and b used in this work is given by

(8)

while the Gaussian terms are of the form (Gauss) (r ) = Uab

⎡ (p (i) − r )2 ⎤ ab ⎥ 2wab(i) ⎥⎦ ⎢⎣

(expon) (inter) (Coul) Uab (r ) = Uab (r ) Tc(r ) + [Uab (r )

∑ Aab(i) exp⎢− i=2

(Gauss) + Uab (r )]T (r )

(9)

The Morse-like exponential potential (eq 8) is introduced at low r to slightly enhance the short-range interaction while preventing atomic overlap. In practice this term is only used for the O−H interaction. In principle the non-Coulomb part of the potential is truncated by a function of the form T(r) = 1 for r < r1, T(r) = 0.5[1 + cos π((r − r1)/(r2 − r1))] for r1 < r < r2, and T(r) = 0 for r > r2 where r1 = 9.0 Å and r2 = 17.5 Å.37 In practice, however, all the exponentials and Gaussians used have positions much smaller than r1 = 9.0 Å so that the effect of this truncation function is negligible. Table 1 lists the individual atom parameters used in the definition of the electrostatic and intramolecular potentials.

(10)

where the definition of the relevant terms has been given in the previous paragraphs. Figure 1 shows the site−site potentials for

Table 1. Individual Atom Parameters Used To Define the Intra- and Intermolecular Interaction Potentials atomic charge [e] charge radius [Å] mass [amu]

O

H

−0.8 1.0 16.0

+0.4 0.5 2.0

Figure 1. Site−site water potentials used as the reference potential in the present work.

water at 283 K that result from this parametrization. Note that, as anticipated, the potentials remain finite at r = 0, and with the exception of the exponential term in the O−H potential and Gaussian term at r = 0 in the O−O potential, the other Gaussian terms are too small to be visible on the scale of this graph. Running the Computer Simulations. For the liquid water simulation at 283 K, 1944 water molecules were placed in a cubic simulation box of dimension 38.779 827 Å, giving an atomic number density of 0.1 atom/Å3. The usual periodic boundary conditions and minimum image convention were applied. For the ice simulation at 258 K, there is a hexagonal box of dimension a = 40.635 910 Å and c = 44.105 156 Å, these numbers being derived from the known lattice constants32 at this temperature, with nine unit cells along each of the a and b axes and six unit cells along the c axis. In this case the center of mass of each water molecule was “tethered” to its lattice positions using a harmonic force field similar to eq 6, with the strength of the tethering constant set to give an O−O radial distribution function reasonably consistent with the measured “Null” ice data, which is sensitive only to this distribution. This tethering potential is needed to prevent the ice lattice from disordering too much because of the random nature of the Monte Carlo moves. The reference potential was set to be the same as the water reference potential as above, but by itself, this potential was not sufficient to drive the structure to give the best possible fit to the scattering data, at least in the time scale available to perform the simulations. As a result, the empirical potential was allowed to increase from zero in this case until

Table 2. Bond Lengths Used To Define the Structure of the Water Molecule atom pair

bond length [Å]

O−H H−H

0.976 1.550

Table 2 lists the set of intramolecular bond lengths used to describe the water molecule, and Table 3 lists the set of exponential and Gaussian parameters for each of the three Table 3. Parameters for the Exponential and Gaussian Terms in the Site−Site Intermolecular Potentialsa atom pair

position [Å]

width [Å]

amplitude [kJ/mol]

O−O

0.00* 0.00 3.30 4.30 1.48* 1.86 2.70 4.40 0.00* 0.00 2.90

0.00* 0.83 0.40 0.55 0.40* 0.32 0.30 0.50 0.00* 0.51 0.30

0.0* +2250.0 −3.0 −3.4 +3.0* −5.6 +2.7 +0.4 0.0* +500.0 +1.0

O−H

H−H

a

The parameters of the exponential terms are marked by asterisks (∗). D



Article

further increases in its amplitude did not produce any improvement in the fit to the data. At the same time the Bragg intensities were calculated exactly out to Q = 3.5 Å−1, using the method described in ref 38. Beyond that the differential cross section was estimated by direct Fourier transform of the simulated radial distribution functions. Each Bragg peak was broadened by a Gaussian resolution function appropriate to SANDALS.

■

RESULTS Water at 283 K. Figure 2 shows the simulated total interference differential cross sections using the reference

Figure 3. Fourier transforms (eq 11) of the simulations and data for water at 283 K in Figure 2, using the same Qmin and Qmax for both. The respective curves are shifted upward for clarity. The inset shows the Fourier transform of the simulation and data for the “Null” water sample; this sample is dominated by the O−O pair distribution function.

Finally the site−site radial distribution functions from these simulations and data are shown in Figure 4. Good overlap

Figure 2. Simulated total interference differential cross sections for H2O (bottom), “Null”, HDO, and D2O at 283 K (lines), using the intermolecular potentials shown in Figure 1, compared to the new neutron diffraction data (points). The respective curves are shifted upward for clarity.

potential shown in Figure 1. It can be seen that the simulation closely matches the data over a broad range of Q, both qualitatively and quantitatively. Inevitably there are small discrepancies in places: these could probably be removed with further refinement of the potential parameters, but the present potential is sufficiently accurate to show that the scattering data can be simulated with an effective pairwise additive potential. The most difficult area is at low Q where the inelasticity effects manifest themselves most strongly in energy dispersive data.39 It is not impossible that the inelasticity effects are still not fully subtracted in this region, introducing some uncertainty in the low Q behavior of the interference differential cross sections. Comparing the same data and simulation in r space is also instructive (Figure 3). Here the Fourier transform G(r) is defined as G (r ) =

1 2π 2ρ

∫Q

Q max min

Q 2 I (Q )

sin Qr dQ Qr

Figure 4. Simulated site−site radial distribution functions for water at 283 K (lines) compared to estimates of these functions from the scattering data (points) using the method described in ref 37 with a feedback factor f = 0.95. The dashed lines show the intramolecular contributions to the H−H and O−H distributions.

between simulation and data is seen for all three functions. In particular it is seen that there is significant penetration of the H−H intermolecular functions into the region where the corresponding intramolecular function still has intensity. This feature is apparently required by the present data in the region 1 < r < 2 Å (see Figure 3). Indeed, if anything, the simulation has slightly overestimated the position of the first intermolecular O−H peak, suggesting there might be some penetration of that function into the intramolecular region as well. It is important to recognize the distinction between the total G(r) shown in Figure 3 and the site−site distributions shown in Figure 4. The total G(r), eq 11, is obtained by Fourier transforming both the data and simulation over the same Q range, with the object of canceling out features that arise from the limitation of the data to finite range of Q. Hence, if the simulation fits the data perfectly, the lines and dots in Figure 3 would coincide. The fact that they do not fully overlap means

(11)

with I(Q) being the differential cross sections shown in Figure 2 and ρ the atomic number density. The same Q range has been used for both simulation and data so that the G(r) functions from each source are directly comparable. Once again the close agreement between simulation and data is clear. In particular it can be seen from the inset of this figure that the simulation of the “Null” sample, which is essentially the O−O pair correlation function, is extremely close to the data, within the usual Fourier truncation oscillations. E



Article

that there are residual discrepancies between data and fit highlighted that may not be evident in the Q-space presentation. On the other hand the solid lines in Figure 4 are the simulated functions and should have no associated truncation effects, whereas the dots are based on Fourier transforms of the data so that even if a perfect fit to the scattering data was achieved, the simulation (lines) and the estimated functions (dots) in this figure might not necessarily agree. Ice at 258 K. The sequence of Figures 2−4 is now repeated but for the structure refinement of the ice data at 258 K, Figures 5−7. Given the substantial increase in the amount of

Figure 6. Fourier transforms (eq 11) of the simulations and data for ice at 258 K in Figure 5, using the same Qmin and Qmax for both. The respective curves are shifted upward for clarity. The inset shows the Fourier transform of the simulation and data for the “Null” ice sample; this sample is dominated by the O−O pair distribution function.

Figure 5. Simulated total interference differential cross sections for H2O (bottom), “Null”, HDO, and D2O ice at 258 K (lines), compared to recent neutron diffraction data (points). The respective curves are shifted upward for clarity. The inset shows the same data but in a narrower Q region, highlighting the region of the Bragg peaks.

structural detail in these ice data from the Bragg peaks, the fits to the scattering data, Figure 5, are not as good as for water, but most of the Bragg peaks are captured qualitatively and many of them quantitatively. Although the sample was cycled repeatedly between liquid and crystal, it is still possible that preferred orientation exists in the crystal data due to the thin, flat plate nature of the container. Fortunately in r space, Figure 6, preferred orientation or texture is less of an issue, since local order is less affected by crystallite orientation. Here again it is found that the simulation captures to a reasonable degree the quantitative and qualitative features of the diffraction data. The worst discrepancy occurs in the region 1 Å < r < 2 Å and appears to be associated with the O−H radial distribution function, as can be seen by inspecting Figure 7 in this region. From this figure it might be concluded that the simulation has underestimated the O−H pair distribution function in this region. However, inspection of the running coordination numbers, Nab(r), which is the total number of atoms of type b within a distance r of atoms of type a,

∫0

Nab(r ) = 4πρcb

r

r′2 gab(r′) dr′

Figure 7. Simulated site−site radial distribution functions for ice at 258 K (lines) compared to estimates of these functions from the scattering data (points).

Figure 8. Simulated site−site radial distribution functions for ice at 258 K (lines) and running coordination numbers from eq 12 (points). Each curve has been shifted vertically by 2 units for clarity.

(12)

where cb is the atomic fraction of type b atoms, shows, Figure 8, that at the first minimum in gOH(r) near r = 2.6 Å, the OH coordination number has already reached 2, which is the maximum allowed in a hydrogen bonded system, so it is not clear where the extra intensity at lower distances implied by

these data would come from. This discrepancy remains unresolved at this point, since for liquid water, the reference potential on its own had little difficulty putting the required F



Article

and for many of them the simulated intensity is in quantitative agreement, particularly the changes observed when going from D2O to H2O, Figure 5. For other peaks the simulated intensities based on the crystallographic structure differ from the data, an effect that is tentatively assigned to residual texture or preferred orientation in the samples. In real space, Figure 6, these problems do not show up so readily, since they probably affect the long-range structure more than in the local region. Even so there does appear a mismatch between simulation and experiment in the region 1 Å < r < 2 Å where the D2O data need more intensity from the simulation while the H2O data require less scattering density. In fact both issues could be resolved by introducing more intensity in the O−H distribution in this region, since this term has a negative cross section weight for H 2 O. However, as already discussed, the coordination number of the first O−H peak is already 2, Figure 8; that is, it is at its limit for fully hydrogen-bonded ice, so it is unclear where this extra intensity would come from, assuming the water molecules remain undissociated. One possibility is that in ice the water molecule intramolecular bond angle is larger than in water, allowing for a larger range of intramolecular distances and so allowing the required penetration into the intramolecular region. So far various possibilities along these lines have not been tested. In fact a notable feature of this model of ice is its high degree of disorder. This is seen in the fact that between the first and second peaks in the O−O distribution the density does not go to zero as would be expected for an ideal crystal. A similar comment applies to the first two peaks in the O−H distribution, implying there is considerable degree of orientational disorder. The degree of disorder is controlled by the tethering potential in this EPSR model, and it is certainly possible to generate a more ordered structure by making this tethering potential stronger. However, doing so makes the simulated Bragg peaks in Figure 5 become too large compared to the data, so it seems the disorder in the model is an effect of the data rather than of the simulation protocol. It is interesting to compare the local order found here in ice Ih with that found in water. Useful distributions to quantify this order are the triple atom included angle distributions O···O···O and O···O−H, where ··· signifies a noncovalent “bond” between distinct molecules and − signifies an intramolecular covalent bond. For a perfect tetrahedral structure the O···O···O distribution would show a sharp peak near θ = 109.47°, and for linear hydrogen bonds the O···O−H distribution would have two peaks, one at 0° and one at ∼104°, the latter due to the second hydrogen on each water molecule. For the present purpose two oxygen atoms are regarded as bonded if their separation is 3.2 Å or less, this distance corresponding the point in the water gOO(r) function where the coordination number reaches 4. In fact, for ice, a pronounced peak near θ = 110° is indeed observed in O···O···O, Figure 9a, but it has a significant broadness covering a base of more than 120o. Hence in this ice model, which appears to be consistent with the crystallography, there is an underlying disordered component. For water the same peak is much more diffuse and has a maximum closer to 100°, and there is a significant contribution near 60°, which corresponds to non-hydrogen-bonded close-packed configurations. Perhaps less intuitive are the O···O−H included angle distributions, Figure 9b. In both water and ice the distributions peak at θ = 0°, corresponding to linear hydrogen bonds being the most probable, but the deviations from linearity for this

intensity in this region (Figures 3 and 4). It could be conjectured that it is associated with the method of “tethering” the water molecules to their crystallographic positions adopted in this work, but so far that assertion remains untested. Figure 8 also highlights the fact that in real ice near its melting point, the ideal lattice of each water molecule being surrounded by exactly four neighbors in a tetrahedral manner may not be realized. The requirement for a coordination number of 4 will only be met on average. Instead there will be a distribution of distances and some overlap between the first and second coordination shells. This is seen because none of the O−O, O−H, or H−H running coordination numbers have a flat plateau with zero gradient at any distance.

■

DISCUSSION DAMP Water Potential. The preceding results demonstrate that a simple pairwise additive intermolecular potential, coupled with a realistic model of the intramolecular disorder, can be developed that gives a quantitative account of observed neutron diffraction data from mixtures of heavy and light water. There is no suggestion that this is the only such potential or even the best that could be found, but it does show that such a potential can exist. The current version of this disordered atom molecular potential (DAMP) lacks the obvious London dispersion force40 which is common to all materials, a weakness that will need to be rectified if it were to be used more widely. Indeed the need for small negative amplitude Gaussian components in the O−O interaction potential in the region 3 Å < r < 5 Å, Table 3, is a likely consequence of this omission. The potential deliberately avoids the use of a 1/rn power law repulsive potential at short distances because such potentials tend to give a too sharp first peak in the O−O radial distribution function, but such a potential is needed if divergences associated with the use of Coulomb point charges are to be avoided. They are avoided in the present work by regarding the atomic charge as distributed uniformly over a spherical volume. It is possible the use of the Buckingham exp-6 potential41 as done in ref 34 could be used in place of the Gaussian repulsive core used here, but the divergence of the dispersive force at short distances would also need to be suppressed in that case. Hence further work to improve the parametrization of this potential is anticipated. The site−site radial distribution functions that derive from this potential and the new diffraction data follow closely those found in other recent experiments.27,26 For the H−H function in particular, and to some extent for the O−H function as well, there appears to be some penetration of the intermolecular distributions at low r into the outer part of the intramolecular region. This is necessitated by trying to get the best possible fit in the region 1 Å < r < 3 Å for all four sets of diffraction data, Figure 3. Although it is tempting to suggest this could simply be an artifact of the procedure used to remove the inelastic and single atom scattering from the scattering data, in fact such effects tend to be slowly varying in r space, and since, in the present case, the area of the O−H intramolecular peak is correctly captured by the data for both heavy and light water, it is difficult to see how the data can be correct in that region of r but incorrect at nearby r values. Ice Structure Compared to Water Structure. The present attempt to build a structure of ice consistent with total scattering data from mixtures of heavy and light water frozen in situ is partly successful. The known Bragg peak positions from crystallography are well reproduced by the data, G



Article

Scheme 1. Schematic of Hydrogen Bonding in Ice (a) and Water (b)a

Figure 9. Simulated O···O···O (a) and O···O−H (b) included angle distribution functions, P(θ), for ice at 258 K (red) and water at 283 K (green). The trivial sin θ dependence in these distributions that would occur if the angles were randomly distributed has been divided out.

peak are larger in ice than in water. On the other hand, the second peak near θ ≈ 110° is stronger in ice than in water. Hence it seems that any given hydrogen bond in ice may be more disordered than in water but that the hydrogen bonding overall is stronger in ice, since the second peak in this distribution is markedly stronger in ice. A possible rationalization of this result can be made with reference to following schematic. In Scheme 1a a central water molecule, C, is surrounded by four water molecules A, B, D, and E in a near tetrahedral arrangement, forming four hydrogen bonds in total. If the bond A···C−HA is linear (0°), then the corresponding bond angle A···C−HB is unlikely to be at the tetrahedral angle because the intramolecular angle is less than the required 109.47° for this to happen. More importantly the hydrogen atoms will be subject to the ever present zero point disorder, so the chances of the two O−H bonds on molecule C lining up in this way are small. Hence there will always be some width to the O···O−H angles, but since each molecule is held in a tetrahedral lattice by its neighbors, the widths of the distributions of the A···C−HA angles (centered on 0°) and A··· C−HB angles (centered on 109.47°) are likely to be similar. In water on the other hand, Scheme 1b, any given water molecule, C, is unlikely to be bonded to four other water molecules with equally strong hydrogen bonds. If the A···C bond is the strongest and most linear, then it is likely the others will be weaker and may be even broken. The likelihood of the hydrogen atom HB being in a near-linear bond to another molecule at the same time as HA is much lower than in (a); hence, it does not suffer the same degree of restriction as it would in ice. The same argument would apply to the B···C−HA angle if the B···C bond were linear. Another factor to be considered is that according to Figures 4 and 7, the hydrogen atom can approach a neighboring water molecule more closely in water compared to ice: this would also likely result in a more linear hydrogen bond.

a

Solid lines represent stronger hydrogen bonds, and dashed lines represent weaker hydrogen bonds. The red and blue circles represent oxygen, while the light gray circles represent hydrogen.

There appear to have been surprisingly few total scattering studies of ice of the kind described here, and with which these results can be compared. Nield and co-workers studied the diffuse scattering in ice Ih single crystals and modeled the data with RMC.42−44 Although some doubt is expressed about the consistency of the outcomes of their analysis, nevertheless the values quoted for the O−O root-mean-square displacement and for the average O···O−H angles are not inconsistent with the distributions shown here. In another study, Teimleitner and others45 performed a total neutron scattering study on heavy water−ice (with much better instrumental resolution than was available in the present study but over a much smaller range of Q) and modeled the total scattering, including the Bragg peaks, with the RMCPOW method.46 Although the chosen temperatures, 120 and 200 K, are lower than the present data and the corresponding peaks in the site−site radial distribution functions correspondingly sharper, the O···O···O and O···O− H angle distributions appear to be quite compatible with those shown in Figure 9. Finally an EPSR study of water over a range of temperatures and of ice at 220 K has recently been submitted,47 where the ice data had been presented and analyzed much earlier23 using a structural model that did not enforce the periodic lattice of ice Ih. The outcomes of this more recent analysis of the earlier data on ice are closely similar to those shown here using the much more recent scattering data. H



■

Article

(9) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (10) Abascal, J. L. F.; Sanz, E.; Fernández, R. G.; Vega, C. A Potential Model for the Study of Ices and Amorphous Water: TIP4P/Ice. J. Chem. Phys. 2005, 122, 234511. (11) Kuharski, R. A.; Rossky, P. J. A Quantum Mechanical Study of Structure in Liquid H2O and D2O. J. Chem. Phys. 1985, 82, 5164− 5177. (12) Paesani, F.; Voth, G. A. The Properties of Water: Insights from Quantum Simulations. J. Phys. Chem. B 2009, 113, 5702−5719. (13) Habershon, S.; Markland, T. E.; Manolopoulos, D. E. Competing Quantum Effects in the Dynamics of a Flexible Water Model. J. Chem. Phys. 2009, 131, 024501. (14) Habershon, S.; Manolopoulos, D. E. Free Energy Calculations for a Flexible Water Model. Phys. Chem. Chem. Phys. 2011, 13, 19714− 19727. (15) Laasonen, K.; Sprik, M.; Parrinello, M.; Car, R. Ab Initio Liquid Water. J. Chem. Phys. 1993, 99, 9080−9089. (16) Fanourgakis, G. S.; Xantheas, S. S. Development of Transferable Interaction Potentials for Water. V. Extension of the Flexible, Polarizable, Thole-Type Model Potential (TTM3-F, V. 3.0) To Describe the Vibrational Spectra of Water Clusters and Liquid Water. J. Chem. Phys. 2008, 128, 074506. (17) Babin, V.; Medders, G. R.; Paesani, F. Toward a Universal Water Model: First Principles Simulations from the Dimer to the Liquid Phase. J. Phys. Chem. Lett. 2012, 3, 3765−3769. (18) Medders, G. R.; Babin, V.; Paesani, F. Development of a “FirstPrinciples” Water Potential with Flexible Monomers. III. Liquid Phase Properties. J. Chem. Theory Comput. 2014, 10, 2906−2910. (19) Narten, A. H.; Thiessen, W. E.; Blum, L. Atom Pair Distribution Functions of Liquid Water at 25°C from Neutron Diffraction. Science 1982, 217, 1033−1034. (20) Soper, A. K.; Silver, R. N. Hydrogen-Hydrogen Pair Correlation Function in Liquid Water. Phys. Rev. Lett. 1982, 49, 471−474. (21) Soper, A. K.; Phillips, M. G. A New Determination of the Structure of Water at 25°C. Chem. Phys. 1986, 107, 47−60. (22) Pusztai, L. Partial Pair Correlation Functions of Liquid Water. Phys. Rev. B 1999, 60, 11851. (23) Soper, A. K. The Radial Distribution Functions of Water and Ice from 220 to 673 K and at Pressures up to 400 MPa. Chem. Phys. 2000, 258, 121−137. (24) Temleitner, L.; Pusztai, L.; Schweika, W. The Structure of Liquid Water by Polarized Neutron Diffraction and Reverse Monte Carlo Modelling. J. Phys.: Condens. Matter 2007, 19, 335207. (25) Zeidler, A.; Salmon, P. S.; Fischer, H. E.; Neuefeind, J. C.; Simonson, J. M.; Lemmel, H.; Rauch, H.; Markland, T. E. Oxygen as a Site Specific Probe of the Structure of Water and Oxide Materials. Phys. Rev. Lett. 2011, 107, 145501. (26) Soper, A. K. The Radial Distribution Functions of Water As Derived from Radiation Total Scattering Experiments: Is There Anything We Can Say for Sure? ISRN Phys. Chem. 2013, 2013, 279463. (27) Skinner, L. B.; Huang, C.; Schlesinger, D.; Pettersson, L. G. M.; Nilsson, A.; Benmore, C. J. Benchmark Oxygen-Oxygen PairDistribution Function of Ambient Water from X-ray Diffraction Measurements with a Wide Q-Range. J. Chem. Phys. 2013, 138, 074506. (28) Tucker, M. G.; Keen, D. A.; Dove, M. T.; Goodwin, A. L.; Hui, Q. RMCProfile: Reverse Monte Carlo for Polycrystalline Materials. J. Phys.: Condens. Matter 2007, 19, 335218. (29) Funnell, N. P.; Dove, M. T.; Goodwin, A. L.; Parsons, S.; Tucker, M. G. Local Structure Correlations in Plastic Cyclohexane a Reverse Monte Carlo Study. J. Phys.: Condens. Matter 2013, 25, 454204. (30) Steinczinger, Z.; Pusztai, L. Comparison of the TIP4P-2005, SWM4-DP and BK3 Interaction Potentials of Liquid Water with Respect to their Consistency with Neutron and X-Ray Diffraction Data of Pure Water. 2013, arXiv:1312.4557.

CONCLUSION The preceding work has described a neutron total scattering experiment on water at 283 K and ice at 258 K, using H2O, D2O, and two mixtures of these liquids. The data were analyzed using the EPSR technique, but for the water data a unique disordered atom molecular potential, DAMP, was set up that was able to describe the experimental data without the need for an additional empirical potential. The ice data were analyzed using the same approach, using the potential created in the water study as the starting point. A model based on the known crystal structure was then perturbed using the EPSR technique to give the best possible fit to the data. Some discrepancies between data and simulation remained which currently are not fully understood because they would require breaking the long accepted wisdom that in ice each water molecule has four neighbors in a near tetrahedral arrangement. Some intriguing comparisons of the nonlinearity of hydrogen bonds in ice and water were made, which can be heuristically reconciled by reference to simple structural ideas. The use of EPSR to refine a crystal structure is a relatively new idea and one that presents practical difficulties in execution. The height of each Bragg peak is sensitive to the degree of disorder present in the model, and it is easy to melt the lattice in a Monte Carlo simulation, hence the use of tethering constraints on the molecular positions in the present work (but not on the molecular orientations). On the other hand total scattering from a crystal, involving as it does the study of the pair correlation function for a material instead of the more common single particle density function, can give a different structural insight on the local order from that obtained by conventional crystallographic techniques.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The author thanks Sam Callear for help with running the SANDALS experiments.

■

REFERENCES

(1) Finney, J. L. Water? What’s So Special About It? Philos. Trans. R. Soc. London, Ser. B 2004, 359, 1145−1165. (2) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of Repulsive Forces in Determining the Equilibrium Structure of Simple Liquids. J. Chem. Phys. 1971, 54, 5237−5247. (3) Hansen, J. P.; MacDonald, I. R. Theory of Simple Liquids; Academic Press: London, 1986. (4) Gray, C. G.; Gubbins, K. E. Theory of Molecular Fluids. Vol. 1: Fundamentals; Oxford University Press: Oxford, U.K., 1984. (5) Bernal, J. D.; Fowler, R. H. A Theory of Water and Ionic Solution, with Particular Reference to Hydrogen and Hydroxyl Ions. J. Chem. Phys. 1933, 1, 515−548. (6) Stillinger, F. H.; Rahman, A. Improved Simulation of Liquid Water by Molecular Dynamics. J. Chem. Phys. 1974, 60, 1545−1557. (7) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (8) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. I



Article

(31) Rottger, K.; Endriss, A.; Ihringer, J.; Doyle, S.; Kuhs, W. F. Lattice Constants and Thermal Expansion of H2O and D2O Ice Ih between 10 and 265 K. Acta Crystallogr., Sect. B: Struct. Sci. 1994, 50, 644−648. (32) Rottger, K.; Endriss, A.; Ihringer, J.; Doyle, S.; Kuhs, W. F. Lattice Constants and Thermal Expansion of H2O and D2O Ice Ih between 10 and 265 K. Addendum. Acta Crystallogr., Sect. B: Struct. Sci. 2012, 68, 91−91. (33) Narten, A. H.; Levy, H. A. Liquid Water: Molecular Correlation Functions from X-ray Diffraction. J. Chem. Phys. 1971, 55, 2263−2269. (34) Paricaud, P.; Předota, M.; Chialvo, A. A.; Cummings, P. T. From Dimer to Condensed Phases at Extreme Conditions: Accurate Predictions of the Properties of Water by a Gaussian Charge Polarizable Model. J. Chem. Phys. 2005, 122, 244511. (35) Chen, Y.-G.; Weeks, J. D. Local Molecular Field Theory for Effective Attractions between Like Charged Objects in Systems with Strong Coulomb Interactions. Proc. Nat. Acad. Sci. U.S.A. 2006, 103, 7560−7565. (36) Rodgers, J. M.; Hu, Z.; Weeks, J. D. On the Efficient and Accurate Short-Ranged Simulations of Uniform Polar Molecular Liquids. Mol. Phys. 2011, 109, 1195−1211. (37) Soper, A. K. Partial Structure Factors from Disordered Materials Diffraction Data: An Approach Using Empirical Potential Structure Refinement. Phys. Rev. B 2005, 72, 104204. (38) Soper, A. K.; Page, K.; Llobet, A. Empirical Potential Structure Refinement of Semi-Crystalline Polymer Systems: Polytetrafluoroethylene and Polychlorotrifluoroethylene. J. Phys.: Condens. Matter 2013, 25, 454219. (39) Soper, A. K. Inelasticity Corrections for Time-of-Flight and Fixed Wavelength Neutron Diffraction Experiments. Mol. Phys. 2009, 107, 1667−1684. (40) London, F. The General Theory of Molecular Forces. Trans. Faraday Soc. 1937, 33, 8b−26. (41) Buckingham, R. A. The Classical Equation of State of Gaseous Helium, Neon and Argon. Proc. R. Soc. A 1938, 264−283. (42) Nield, V. M.; Whitworth, R. W. The Structure of Ice Ih from Analysis of Single-Crystal Neutron Diffuse Scattering. J. Phys.: Condens. Matter 1995, 7, 8259−8271. (43) Beverley, M. N.; Nield, V. M. Extensive Tests on the Application of Reverse Monte Carlo Modelling to Single-Crystal Neutron Diffuse Scattering from Ice Ih. J. Phys.: Condens. Matter 1997, 9, 5145−5156. (44) Beverley, M. N.; Nield, V. M. Analysis of Single-Crystal Neutron Diffuse Scattering from Ice Ih. J. Phys. Chem. B 1997, 101, 6188−6191. (45) Temleitner, L.; Pusztai, L. Investigation of Structural Disorder in Ice Ih Using Neutron Diffraction and Reverse Monte Carlo Modelling. In Physics and Chemistry of Ice; Kuhs, W. F., Ed.; RSC Publishing: Cambridge, U.K., 2007 (46) Mellergard, A.; McGreevy, R. L. Reverse Monte Carlo Modelling of Neutron Powder Diffraction Data. Acta Crystallogr., Sect. A: Found. Crystallogr. 1999, 55, 783−789. (47) Soper, A. K. Water and Ice Structure in the Range 220 - 365K from Radiation Total Scattering Experiments. 2014, arXiv:1411.1322.

J


Disordered Atom Molecular Potential for Water Parameterized against

Recommend Documents