^{1,a)}, D. J. Anick

^{2}, P. K. Mankoo

^{3,b)}and G. F. Reiter

^{1}

### Abstract

We present an empirical flexible and polarizable water model which gives an improved description of the position, momentum, and dynamical (spectroscopic) distributions of H nuclei in water. We use path integral molecular dynamics techniques in order to obtain momentum and position distributions and an approximate solution to the Schrödinger equation to obtain the infrared (IR) spectrum. We show that when the calculated distributions are compared to experiment the existing empirical models tend to overestimate the stiffness of the H nuclei involved in H bonds. Also, these models vastly underestimate the enormous increase in the integrated IR intensity observed in the bulk over the gas-phase value. We demonstrate that the over-rigidity of the OH stretch and the underestimation of intensity are connected to the failure of existing models to reproduce the correct monomerpolarizabilitysurface. A new model, TTM4-F, is parametrized against electronic structure results in order to better reproduce the polarizabilitysurface. It is found that TTM4-F gives a superior description of the observed spectroscopy, showing both the correct redshift and a much improved intensity. TTM4-F also has a somewhat improved dielectric constant and OH distribution function. It also gives an improved match to the experimental momentum distribution, although some discrepancies remain.

This work was supported by the DOE, Office of Basic Energy Sciences under Contract No. DE-FG02-03ER46078. The authors thank T. Keyes of the Department of Chemistry, Boston University for discussions regarding spectral calculations, F. Paesani of the Department of Chemistry, University of Utah for discussions regarding calculations of the dielectric constant, and B. R. Johnson of the Department of Chemistry, Rice University for discussions regarding solving the rotational Schrödinger equation.

I. INTRODUCTION

II. ANALYSIS OF ELECTROSTATICS

III. MODELS

A. Water models used in this study

B. Parametrization of TTM4-F

IV. ELECTRONIC STRUCTURE

V. DISTRIBUTION FUNCTIONS

VI. DIELECTRIC CONSTANT

VII. SPECTROSCOPY

VIII. SIMULATION DETAILS

IX. RESULTS

A. Monomer properties

B. Cluster properties

C. Bulk properties

X. CONCLUSIONS

A. A comment on the TTM3-F model of Fanourgakis and Xantheas

### Key Topics

- Infrared spectra
- 38.0
- Polarizability
- 34.0
- Polymers
- 29.0
- Hydrogen bonding
- 27.0
- Molecular dynamics
- 24.0

## Figures

water cluster used as a model of the solvated proton. The central proton along the vertex of the three fused pentagons is scanned.

water cluster used as a model of the solvated proton. The central proton along the vertex of the three fused pentagons is scanned.

Derivative of monomer dipole component parallel to OH bond with respect to OH distance as a function of external field parallel to OH bond. A comparison between the electronic structure results of Hermansson (Ref. 31) and models is shown. The units are consistent with Fig. 6 of Hermansson. (Multiply by to convert a.u./Å to D/Å.) Also, the same geometry (104.5°, ) was used as in Hermansson’s paper.

Derivative of monomer dipole component parallel to OH bond with respect to OH distance as a function of external field parallel to OH bond. A comparison between the electronic structure results of Hermansson (Ref. 31) and models is shown. The units are consistent with Fig. 6 of Hermansson. (Multiply by to convert a.u./Å to D/Å.) Also, the same geometry (104.5°, ) was used as in Hermansson’s paper.

Vibrational component of the calculated gas-phase IR spectrum at for various models. Both classical and quantum results are shown. The quantum results were obtained from solving the three-dimensional vibrational Schrödinger equation. Experimental data are only provided in the stretch region. Only the experimentally observed antisymmetric peak is clearly visible as the symmetric peak has an intensity 20 times smaller. To make for an easier comparison with later figures, the intensities are calculated using a bulk liquid density of (noninteracting) monomers. Also, the experimental and quantum results are shown using a histogram of a nonzero width. The experimental spectrum was calculated using integrated IR intensities of 0.226 and quoted by Whalley and Klug (Ref. 78) and using frequencies of 3657 and .

Vibrational component of the calculated gas-phase IR spectrum at for various models. Both classical and quantum results are shown. The quantum results were obtained from solving the three-dimensional vibrational Schrödinger equation. Experimental data are only provided in the stretch region. Only the experimentally observed antisymmetric peak is clearly visible as the symmetric peak has an intensity 20 times smaller. To make for an easier comparison with later figures, the intensities are calculated using a bulk liquid density of (noninteracting) monomers. Also, the experimental and quantum results are shown using a histogram of a nonzero width. The experimental spectrum was calculated using integrated IR intensities of 0.226 and quoted by Whalley and Klug (Ref. 78) and using frequencies of 3657 and .

Top left: Comparison of proton scan for monomer and solvated proton from Fig. 1 using the surface. Also shown are the eigenvalues and eigenfunctions from solving the 1D Schrödinger equation for a proton in each well. The numbers show the transition frequency in wavenumbers. Top middle to bottom right: Comparison of proton scan for solvated proton from Fig. 1 using model and surfaces. Also shown are the eigenvalues and eigenfunctions from solving the 1D Schrödinger equation for a proton in each well. The numbers show the transition frequency in wavenumbers.

Top left: Comparison of proton scan for monomer and solvated proton from Fig. 1 using the surface. Also shown are the eigenvalues and eigenfunctions from solving the 1D Schrödinger equation for a proton in each well. The numbers show the transition frequency in wavenumbers. Top middle to bottom right: Comparison of proton scan for solvated proton from Fig. 1 using model and surfaces. Also shown are the eigenvalues and eigenfunctions from solving the 1D Schrödinger equation for a proton in each well. The numbers show the transition frequency in wavenumbers.

Top: Bulk liquid OO and OH distribution functions. Comparison of experimental data of Soper *et al.* (Ref. 1) with classical and path integral molecular dynamics results using the TTM4-F model. Bottom: Comparison of experimental (Ref. 1) bulk liquid OO and OH distribution functions and path integral simulation results using various models.

Top: Bulk liquid OO and OH distribution functions. Comparison of experimental data of Soper *et al.* (Ref. 1) with classical and path integral molecular dynamics results using the TTM4-F model. Bottom: Comparison of experimental (Ref. 1) bulk liquid OO and OH distribution functions and path integral simulation results using various models.

Dielectric constant calculated using path integral molecular dynamics for the RWK2 model. Results are shown for one (classical), two, four, and eight replicas.

Dielectric constant calculated using path integral molecular dynamics for the RWK2 model. Results are shown for one (classical), two, four, and eight replicas.

Comparison of experimental and calculated IR spectra. Experimental data for the liquid from Bertie and Lan (Ref. 95). Experimental data for ice Ih taken from refractive index measurements of Clapp *et al.* (Ref. 96) at . Calculated results are shown from both classical (left) and wavefunction calculations (right).

Comparison of experimental and calculated IR spectra. Experimental data for the liquid from Bertie and Lan (Ref. 95). Experimental data for ice Ih taken from refractive index measurements of Clapp *et al.* (Ref. 96) at . Calculated results are shown from both classical (left) and wavefunction calculations (right).

Experimental and calculated proton momentum distribution functions in bulk ice Ih. Results for TIP3P-F, RWK2, and TTM2-F are taken from Burnham *et al.* (Ref. 29).

Experimental and calculated proton momentum distribution functions in bulk ice Ih. Results for TIP3P-F, RWK2, and TTM2-F are taken from Burnham *et al.* (Ref. 29).

Comparison of calculated momentum distributions from TTM4-F using two different calculation methods: Path integral and solving the 1D Schrödinger equation.

Comparison of calculated momentum distributions from TTM4-F using two different calculation methods: Path integral and solving the 1D Schrödinger equation.

Calculated momentum distribution along the stretch direction for a given proton from solution of the 1D Schrödinger equation along the stretch direction vs H-bond parameter for TTM4-F for a variety of water phases. Also shown are experimental data points for the momentum distribution, using calculated averages from the simulation for the H-bond parameter.

Calculated momentum distribution along the stretch direction for a given proton from solution of the 1D Schrödinger equation along the stretch direction vs H-bond parameter for TTM4-F for a variety of water phases. Also shown are experimental data points for the momentum distribution, using calculated averages from the simulation for the H-bond parameter.

KE vs from solution of the 1D Schrödinger equation along the stretch direction for TTM4-F for a variety of water phases. Also shown are experimental data points using estimates for the KE from observed momentum distributions and frequencies from observed IR spectra in the OH stretch region. A value of for ice VI is taken from the data of Bertie *et al.* (Ref. 97) at . We were unable to find experimental IR data for the particular supercritical water thermodynamic point.

KE vs from solution of the 1D Schrödinger equation along the stretch direction for TTM4-F for a variety of water phases. Also shown are experimental data points using estimates for the KE from observed momentum distributions and frequencies from observed IR spectra in the OH stretch region. A value of for ice VI is taken from the data of Bertie *et al.* (Ref. 97) at . We were unable to find experimental IR data for the particular supercritical water thermodynamic point.

## Tables

Model parameters. Coefficients [of Eq. (13)] are in units of (kcal/mol).

Model parameters. Coefficients [of Eq. (13)] are in units of (kcal/mol).

Electrostatic properties of the gas-phase monomer. the axis points out of the plane of the molecule, the axis is along the molecular bisection, and the axis is perpendicular to both and and lies in the molecular plane. Data in parentheses refer to coupled cluster results using the basis set comprising of 157 Gaussian-type functions from Avila (Ref. 86).

Electrostatic properties of the gas-phase monomer. the axis points out of the plane of the molecule, the axis is along the molecular bisection, and the axis is perpendicular to both and and lies in the molecular plane. Data in parentheses refer to coupled cluster results using the basis set comprising of 157 Gaussian-type functions from Avila (Ref. 86).

Comparison of various bulk ice , ambient liquid , and gas-phase properties from various water models with experimental data (where available).

Comparison of various bulk ice , ambient liquid , and gas-phase properties from various water models with experimental data (where available).

Intramolecular geometries for various water models in the gas phase , ambient liquid water , and ice using both classical and path integral simulation methods.

Intramolecular geometries for various water models in the gas phase , ambient liquid water , and ice using both classical and path integral simulation methods.

Binding energies for water clusters. MP2/CBS results from Xantheas *et al.* (Ref. 40).

Binding energies for water clusters. MP2/CBS results from Xantheas *et al.* (Ref. 40).

Article metrics loading...

Full text loading...

Commenting has been disabled for this content