^{1}and András Baranyai

^{1}

### Abstract

Based on extensive studies of existing potentials we propose a new molecular model for water. The new model is rigid and contains three Gaussian charges. Contrary to other models, all charges take part in the polarization of the molecule. They are connected by harmonic springs to their gas-phase positions: the negative charge to a prescribed point on the main axis of the molecule; the positive charges to the hydrogens. The mechanical equilibrium between the electrostatic forces and the spring forces determines the polarization of the molecule which is established by iteration at every timestep. The model gives excellent estimates for ambient liquid properties and reasonably good results from high-pressure solids to gas-phase clusters. We present a detailed description of the development of this model and a large number of calculated properties compared to the estimates of the nonpolarizable TIP4P/2005 [J. L. F. Abascal and C. Vega, J. Chem. Phys.123, 234505 (Year: 2005)10.1063/1.2121687], the polarizable GCPM [P. Paricaud, M. Predota, A. A. Chialvo, and P. T. Cummings, J. Chem. Phys.122, 244511 (Year: 2005)10.1063/1.1940033], and our earlier BKd3 model [P. T. Kiss and A. Baranyai, J. Chem. Phys.137, 084506 (Year: 2012)10.1063/1.4746419]. The best overall performance is shown by the new model.

A.B. gratefully acknowledges the support of OTKA Grant No. K84382.

I. INTRODUCTION

II. THE MODEL

A. Geometry

B. Charge distribution

C. Nonelectrostatic interactions

D. Polarization

III. SIMULATION DETAILS

A. Handling of the charge-charge interactions

B. Integration of the equations of motion

IV. PARAMETRIZATION OF THE MODEL

V. EVALUATION OF THE RESULTS

A. Water at ambient conditions

B. Properties in terms of the temperature

C. Vapor-liquid equilibrium and critical properties

D. Gas phase clusters

E. Second virial coefficient

F. Melting properties

G. Ice phases

VI. CONCLUSIONS

### Key Topics

- Polarization
- 33.0
- Ice
- 27.0
- Polarizability
- 22.0
- Thermodynamic properties
- 16.0
- Dielectric constant
- 15.0

## Figures

Schematic depiction of the BK3 model.

Schematic depiction of the BK3 model.

System size-dependence of the self-diffusion coefficient is shown. The simulated systems contained 250, 432, 600, 1000, and 4000 molecules. The circles and diamonds refer to the uncorrected and corrected values.

System size-dependence of the self-diffusion coefficient is shown. The simulated systems contained 250, 432, 600, 1000, and 4000 molecules. The circles and diamonds refer to the uncorrected and corrected values.

(a) The convergence of the dielectric constant at different temperatures. (b) The temperature dependence of the dielectric constant of the BK3 model. The error bars were estimated from five runs for each temperature.

(a) The convergence of the dielectric constant at different temperatures. (b) The temperature dependence of the dielectric constant of the BK3 model. The error bars were estimated from five runs for each temperature.

Oxygen–oxygen (top), oxygen–hydrogen (middle), and hydrogen–hydrogen (bottom) pair-correlation functions at ambient conditions, 298 K and 1 bar. The experimental data are from Ref. 91 .

Oxygen–oxygen (top), oxygen–hydrogen (middle), and hydrogen–hydrogen (bottom) pair-correlation functions at ambient conditions, 298 K and 1 bar. The experimental data are from Ref. 91 .

The density of the liquid water at 1 bar as a function of temperature.

The density of the liquid water at 1 bar as a function of temperature.

(a) Self-diffusion coefficients as function of temperature. Finite size corrections are included only for BK3 model. (b) Shear viscosity of the BK3 model as a function of temperature. The error bars were estimated from using different off-diagonal components for the calculation.

(a) Self-diffusion coefficients as function of temperature. Finite size corrections are included only for BK3 model. (b) Shear viscosity of the BK3 model as a function of temperature. The error bars were estimated from using different off-diagonal components for the calculation.

The temperature dependence of the isobaric heat capacity at 1 bar. The experimental line was shifted up by 10 J K−1 mol−1 for better comparison.

The temperature dependence of the isobaric heat capacity at 1 bar. The experimental line was shifted up by 10 J K−1 mol−1 for better comparison.

Thermal expansion coefficient as a function of temperature at 1 bar. The crossing of the curves with the T axis gives the T md .

Thermal expansion coefficient as a function of temperature at 1 bar. The crossing of the curves with the T axis gives the T md .

The temperature dependence of the compressibility at 1 bar. The marks are the simulated values which have relatively large statistical noise. The lines are fitted fourth order polynomials.

The temperature dependence of the compressibility at 1 bar. The marks are the simulated values which have relatively large statistical noise. The lines are fitted fourth order polynomials.

Vapor-liquid coexistence curves. The crosses indicate the corresponding critical points.

Vapor-liquid coexistence curves. The crosses indicate the corresponding critical points.

(a) Equilibrium vapor pressures as a function of temperature. The crosses indicate the critical points. (b) Same as in (a) but using the corresponding reduced temperature (T/T c ) as independent variable.

(a) Equilibrium vapor pressures as a function of temperature. The crosses indicate the critical points. (b) Same as in (a) but using the corresponding reduced temperature (T/T c ) as independent variable.

Surface tension as a function of temperature. With the exception of the BKd3 model the tail correction of the dispersion forced is included.

Surface tension as a function of temperature. With the exception of the BKd3 model the tail correction of the dispersion forced is included.

Vaporization enthalpies as a function of temperature. The marks indicate simulated points and the lines are help for the eye.

Vaporization enthalpies as a function of temperature. The marks indicate simulated points and the lines are help for the eye.

Temperature dependence of the second virial coefficient. In the inset the high temperature region is enlarged. The meaning of the colors is the same as in other figures: green-GCPM; blue-TIP4P/2005; red- BK3; pink-BKd3.

Temperature dependence of the second virial coefficient. In the inset the high temperature region is enlarged. The meaning of the colors is the same as in other figures: green-GCPM; blue-TIP4P/2005; red- BK3; pink-BKd3.

The density of the ice VII polymorph as a function of pressure at 300 K. The experimental data were taken from Ref. 23 . The curves of the models were calculated by us.

The density of the ice VII polymorph as a function of pressure at 300 K. The experimental data were taken from Ref. 23 . The curves of the models were calculated by us.

## Tables

Parameters of the BK3 model (see text for details.)

Parameters of the BK3 model (see text for details.)

Properties of ambient water (298 K, 1 bar): internal energy U, density ρ, self-diffusion coefficient D, shear viscosity η, heat capacity C p , compressibility κ T , thermal expansion coefficient α p , dielectric constant ɛ, and average dipole moment ⟨μ⟩. The temperature of maximum density (T md ) is also shown. The results of the BK3 model were calculated in this work, the numbers in parentheses are the estimated errors in the last digit. Data for the TIP4P/2005, GCPM, and BKd3 models were taken from Refs. 6, 7, and 13 , and 74 . Experimental data were taken from Ref. 2 .

Properties of ambient water (298 K, 1 bar): internal energy U, density ρ, self-diffusion coefficient D, shear viscosity η, heat capacity C p , compressibility κ T , thermal expansion coefficient α p , dielectric constant ɛ, and average dipole moment ⟨μ⟩. The temperature of maximum density (T md ) is also shown. The results of the BK3 model were calculated in this work, the numbers in parentheses are the estimated errors in the last digit. Data for the TIP4P/2005, GCPM, and BKd3 models were taken from Refs. 6, 7, and 13 , and 74 . Experimental data were taken from Ref. 2 .

Critical temperature, (T c ), pressure, (P c ), and density, (ρ c ), for different models. The results for the BK3 model were calculated in this work, the results of TIP4P/2005, GCPM, and BKd3 models were taken from Refs. 17 and 7 , and 13 , respectively. Experimental data were taken from Ref. 2 .

Energies (E) and average oxygen–oxygen distances (d OO) of small water clusters for different models. The energies are in kJ/mol and the distances are in Å. In the case of the dimer the characteristic angle, θ (the angle between the OO line and the symmetry axis of the acceptor molecule), and the total dipole moment, (μ tot), are also shown. QC refers to the highest level quantum chemical calculations (Refs. 36 and 84 ). Experimental data were taken from Refs. 35 and 85 .

Energies (E) and average oxygen–oxygen distances (d OO) of small water clusters for different models. The energies are in kJ/mol and the distances are in Å. In the case of the dimer the characteristic angle, θ (the angle between the OO line and the symmetry axis of the acceptor molecule), and the total dipole moment, (μ tot), are also shown. QC refers to the highest level quantum chemical calculations (Refs. 36 and 84 ). Experimental data were taken from Refs. 35 and 85 .

Melting properties of ice Ih at 1 bar for different models. T m is the melting temperature, ρ l and ρ s are the densities of the equilibrium liquid and solid phases, H l and H s are the corresponding enthalpies, ΔH is the melting enthalpy, ΔS it the melting entropy, and dp/dT is the slope of the ice Ih-water equilibrium line at 1 bar. The results of the BK3 model were calculated in this work, data for the TIP4P/2005 and BKd3 models were taken from Refs. 6 and 14 . Experimental data were taken from Ref. 2 .

Melting properties of ice Ih at 1 bar for different models. T m is the melting temperature, ρ l and ρ s are the densities of the equilibrium liquid and solid phases, H l and H s are the corresponding enthalpies, ΔH is the melting enthalpy, ΔS it the melting entropy, and dp/dT is the slope of the ice Ih-water equilibrium line at 1 bar. The results of the BK3 model were calculated in this work, data for the TIP4P/2005 and BKd3 models were taken from Refs. 6 and 14 . Experimental data were taken from Ref. 2 .

Densities of different ice polymorphs are shown in g/cm3. The experimental data were taken from Ref. 23 . The results of the models are calculated by us. The error bar of the calculated densities is ±0.002 g/cm3.

Densities of different ice polymorphs are shown in g/cm3. The experimental data were taken from Ref. 23 . The results of the models are calculated by us. The error bar of the calculated densities is ±0.002 g/cm3.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content