Péter T. Kiss

A systematic development of a polarizable potential of water

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 (2005)], the polarizable GCPM [P. Paricaud, M. Predota, A. A. Chialvo, and P. T. Cummings, J. Chem. Phys. 122, 244511 (2005)], and our earlier BKd3 model [P. T. Kiss and A. Baranyai, J. Chem. Phys. 137, 084506 (2012)]. The best overall performance is shown by the new model.

Clusters of classical water models

The properties of clusters can be used as tests of models constructed for molecular simulation of water. We searched for configurations with minimal energies for a small number of molecules. We identified topologically different structures close to the absolute energy minimum of the system by calculating overlap integrals and enumerating hydrogen bonds. Starting from the dimer, we found increasing number of topologically different, low-energy arrangements for the trimer(3), the tetramer(6), the pentamer(6), and the hexamer(9). We studied simple models with polarizable point dipole. These were the BSV model [J. Brodholt et al., Mol. Phys. 86, 149 (1995)], the DC model [L. X. Dang and T. M. Chang, J. Chem. Phys. 106, 8149 (1997)], and the GCP model [P. Paricaud et al., J. Chem. Phys. 122, 244511 (2005)]. As an alternative the SWM4-DP and the SWM4-NDP charge-on-spring models [G. Lamoureux et al., Chem. Phys. Lett. 418, 245 (2006)] were also investigated. To study the impact of polarizability restricted to the plane of the molecule we carried out calculations for the SPC-FQ and TIP4P-FQ models, too [S. W. Rick et al., J. Chem. Phys. 101, 6141 (1994)]. In addition to them, justified by their widespread use even for near critical or surface behavior calculations, we identified clusters for five nonpolarizable models of ambient water, SPC/E [H. J. C. Berendsen et al., J. Phys. Chem. 91, 6269 (1987)], TIP4P [W. L. Jorgensen et al., J. Chem. Phys. 79, 926 (1983)], TIP4P-EW [H. W. Horn et al., J. Chem. Phys. 120, 9665 (2004)], and TIP4P/2005 [J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005)]. The fifth was a five-site model named TIP5P [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 (2000)]. To see the impact of the vibrations we studied the flexible SPC model. [K. Toukan and A. Rahman, Phys. Rev. B 31, 2643 (1985)]. We evaluated the results comparing them with experimental data and quantum chemical calculations. The position of the negative charge in the models plays a crucial role. In this respect models with SPC geometry provided structures different from the TIP4P-type potentials, including polarizable ones. The TIP4P variants form configurations similar to one another. Results for TIP4P-EW and for TIP4P/2005 were especially close to each other in every respect. This is also true for the BSV and the DC pair. The charge-on-spring models (SWM4-DP and SWM4-NDP) are also very similar to each other, despite the sign exchange of charges on the spring particle and the oxygen. The spherical polarization of water is crucial. Due to the planar polarization of the SPC-FQ and the TIP4P-FQ models, they prefer planar arrangements contrary to other polarizable models and quantum chemical calculations. The tetrahedral geometry of TIP5P stabilizes additional clusters with peculiar geometries and small O-O distances. Inclusion of vibrations causes only insignificant changes in the characteristic geometries but decreases the internal energy relative to its reference rigid version. Comparing with quantum mechanical calculations the GCP model provided the best overall results.

A transferable classical potential for the water molecule.

We developed a new model for the water molecule which contains only three Gaussian charges. Using the gas-phase geometry the dipole moment of the molecule matches, the quadrupole moment closely approximates the experimental values. The negative charge is connected by a harmonic spring to its gas-phase position. The polarized state is identified by the equality of the intermolecular electrostatic force and the spring force acting on the negative charge. In each timestep the instantaneous position of the massless negative charge is determined by iteration. Using the technique of Ewald summation, we derived expressions for the potential energy, the forces, and the pressure for Gaussian charges. The only properties to be fitted are the half-width values of the Gaussian charge distributions and the parameters of the nonelectrostatic repulsion-attraction potential. We determined the properties of gas-phase clusters up to six molecules, the internal energy and density of ambient water and hexagonal ice. We calculated the equilibrium density of ice VII as a function of pressure. As an additional test, we calculated the pair-correlation function, the isotherm compressibility, the heat capacity, and the self-diffusion coefficients for ambient water. As far as we know, this is the first classical model of water which is able to estimate both ends of the phase diagram, the high pressure ice VII, and the gas clusters of water with excellent accuracy.

A new polarizable force field for alkali and halide ions

We developed transferable potentials for alkali and halide ions which are consistent with our recent model of water [P. T. Kiss and A. Baranyai, J. Chem. Phys. 138, 204507 (2013)]. Following the approach used for the water potential, we applied Gaussian charge distributions, exponential repulsion, and r(-6) attraction. One of the two charges of the ions is fixed to the center of the particle, while the other is connected to this charge by a harmonic spring to express polarization. Polarizability is taken from quantum chemical calculations. The repulsion between different species is expressed by the combining rule of Kong [J. Chem. Phys. 59, 2464 (1972)]. Our primary target was the hydration free energy of ions which is correct within the error of calculations. We calculated water-ion clusters up to 6 water molecules, and, as a crosscheck, we determined the density and internal energy of alkali-halide crystals at ambient conditions with acceptable accuracy. The structure of hydrated ions was also discussed.

Sources of the deficiencies in the popular SPC/E and TIP3P models of water

Motivated by the results of Vega et al. [J. Phys. Condens. Matter 20, 153101 (2008)] about the phase diagram of water, and by the results of Kiss and Baranyai [J. Chem. Phys. 131, 204310 (2009)] about the properties of gas-phase clusters, we carried out a comparative study of the structure modeled by SPC∕E and TIP3P interactions in ambient liquid water. The gas-phase clusters of SPC∕E and TIP3P models show erroneous structures, while TIP4P-type models, either polarizable or not, provide qualitatively correct results. The trimers of SPC∕E and TIP3P are planar in gas phase, contrary to experimental and TIP4P-type models. The aim of this study was to see whether traces of these false geometries characteristic to SPC∕E and TIP3P in gas phase can also be found in the liquid phase. For this purpose we selected trimers formed by adjacent neighbors of water molecules in the liquid and calculated their geometrical features. We determined angles formed by the HO bonds of the molecules with OO vectors and with the normal vector of the OOO plane in the selected trimers. Our results showed that, despite high temperature, the SPC∕E and TIP3P water contains larger number of planar arrangements than other TIP4P-type models. Although structural differences presented in this study are small, they are accurately detectable. These results weaken the reliability of studies obtained by the SPC∕E or TIP3P models even in the liquid phase.

Density maximum and polarizable models of water

To estimate accurately the density of water over a wide range of temperatures with a density maximum at 4 °C is one of the most stringent tests of molecular models. The shape of the curve influences the ability to describe critical properties and to predict the freezing temperature. While it was demonstrated that with a proper parameter fit nonpolarizable models can approximate this behavior accurately, it is much more difficult to do this for polarizable models. We provide a short overview of ρ-T diagrams for existing models, then we give an explanation of this difficulty. We present a version of the BK model [A. Baranyai and P. T. Kiss, J. Chem. Phys. 133, 144109 (2010); and ibid. 135, 234110 (2011)] which is capable to predict the density of water over a wide range of temperature. The BK model uses the charge-on-spring method with three Gaussian charges. Since the experimental dipole moment and the geometry is fixed, and the quadrupole moment is approximated by a least mean square procedure, parameters of the repulsion and dispersive attraction forces remained as free tools to match experimental properties. Relying on a simplified but plausible justification, the new version of the model uses repulsion and attraction as functions of the induced dipole moment of the molecule. The repulsive force increases, while the attractive force decreases with the size of the molecular dipole moment. At the same time dipole moment dependent dispersion forces are taking part in the polarization of the molecule. This scheme iterates well and, in addition to a reasonable density-temperature function, creates dipole distributions with accurate estimation of the dielectric constant of the liquid.

Comparison of wind power estimates from the ECMWF reanalyses with direct turbine measurements

Reanalysis data are rarely used for wind power estimates because of the limited spatial and temporal resolution. Here we report on a detailed comparison of wind speed and electric power time series recorded at a continental location in Hungary and estimates provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-40 and Interim databases at nearby grid points. The results show that the temporal behavior is adequately represented in reanalysis records with damped magnitudes, as expected. However, characteristic shape differences in the wind speed histograms for turbine measurements and reanalysis hinder a perfect match of statistics. A satisfying agreement of histograms for measured and modeled output powers is achieved by scaling up surface wind speeds to have the same long time average value as for the turbine records. The presented calibration permits us to provide wind power estimates for large geographic areas, where the wind field is similarly coherent as around the test site.

Testing recent charge-on-spring type polarizable water models. I. Melting temperature and ice properties

We determined the freezing point of eight molecular models of water. All models use the charge-on-spring (COS) method to express polarization. The studied models were the COS/G2, COS/G3 [H. Yu and W. F. van Gunsteren, J. Chem. Phys. 121, 9549 (2004)], the COS/B2 [H. Yu, T. Hansson, and W. F. van Gunsteren, J. Chem. Phys. 118, 221 (2003)], the SWM4-DP [G. Lamoureux, A. D. MacKerell, Jr., and B. Roux, J. Chem. Phys. 119, 5185 (2003)], the SWM4-NDP [G. Lamoureux, E. Harder, I. V. Vorobyov, B. Roux, and A. D. MacKerell, Jr., Chem. Phys. Lett. 418, 245 (2006)], and three versions of our model, the BKd1, BKd2, and BKd3. The BKd1 is the original Gaussian model [P. T. Kiss, M. Darvas, A. Baranyai, and P. Jedlovszky, J. Chem. Phys. 136, 114706 (2012)] with constant polarization and with a simple exponential repulsion. The BKd2 applies field-dependent polarizability [A. Baranyai and P. T. Kiss, J. Chem. Phys. 135, 234110 (2011)], while the BKd3 model has variable size to approximate the temperature-density (T-ρ) curve of water [P. T. Kiss and A. Baranyai, J. Chem. Phys. 137, 084506 (2012)]. We used the thermodynamic integration (TI) and the Gibbs-Helmholtz equation to determine the equality of the free energy for liquid water and hexagonal ice (Ih) at 1 bar. We used the TIP4P and the SPC/E models as reference systems of the TI. The studied models severely underestimate the experimental melting point of ice Ih. The calculated freezing points of the models are the following: COS/G2, 215 K; COS/G3, 149 K; SWM4-DP, 186 K; BKd1, 207 K; BKd2, 213 K; BKd3, 233 K. The freezing temperature of the SWM4-NDP system is certainly below 120 K. It might even be that the water phase is more stable than the ice Ih at 1 bar in the full temperature range. The COS/B2 model melts below 100 K. The best result was obtained for the BKd3 model which indicates that correct description of the (T-ρ) curve improves the estimation of the freezing point. We also determined and compared the densities of high-pressure polymorphs of ice for these models.

Efficient Handling of Gaussian Charge Distributions: An Application to Polarizable Molecular Models.

We present a mesh-based Ewald summation method that is suitable for the calculation of the electrostatic interaction between Gaussian charge distributions, instead of point charges. As an application, we implemented the method in the Gromacs simulation package and tested it with a polarizable water model, showing that the interaction between Gaussian charge distributions can be computed with a small (10%) additional computational cost with respect to the point charge case. In addition, since the performance of polarizable models is strongly influenced by the number of iterations required for the self-consistent field solution, we tested also the Always Stable Predictor-Corrector (ASPC) method of Kolafa (Kolafa, J. J. Comp. Chem. 2003, 25, 335) as an alternative to the steepest descent (SD) based algorithm with predictor implemented in the Gromacs, and found that it speeds up the integration of the equations of motion by a factor of 1.6-2.0, depending on the target model.

Anomalous properties of water predicted by the BK3 model

Recently, we proposed a new model for water [P. T. Kiss and A. Baranyai, J. Chem. Phys. 138, 204507 (2013)]. We presented a detailed description of the development of this classical, polarizable model, and a large number of calculated properties. The model provided excellent estimates for ambient liquid properties and reasonably good results from high-pressure solids to gas-phase clusters. In this paper we present results of extensive calculations for temperature-dependent water anomalies in terms of the pressure. The calculated isobars of the temperature-density and the self-diffusion diagrams provide excellent estimates of the experimental values. The estimated compressibility isobars perfectly match the experimental ones if we shift our numbers by ∼10 K upwards. The calculated pressure-dependent viscosity values are excellent at higher temperatures and qualitatively correct at lower temperatures.