András Baranyai

In classical statistical mechanics the trajectory in phase space represents the propagation of a classical Hamiltonian system. While trajectories play a key role in chaotic system theory, exploitation of a single trajectory has yet to be considered. This work shows that for ergodic dynamical systems the dynamical temperature can be derived using phase space trajectories.

We carried out extensive calculations of liquid water at different temperatures and pressures using the BK3 model suggested recently [P. T. Kiss and A. Baranyai, J. Chem. Phys. 138, 204507 (2013)]. In particular, we were interested in undercooled regions to observe the propensity of water to form tetrahedral coordination of closest neighbors around a central molecule. We compared the found tetrahedral order with the number of hydrogen bonds and with the partial pair correlation functions unfolded as distributions of the closest, the second closest, etc. neighbors. We found that contrary to the number of hydrogen bonds, tetrahedrality changes substantially with state variables. Not only the number of tetrahedral arrangements increases with lowering the pressure, the density, and the temperature but the domain size of connecting tetrahedral structures as well. The difference in tetrahedrality is very pronounced between the two sides of the Widom line and even more so between the low density amorphous (LDA) and high density amorphous (HDA) phases. We observed that in liquid water and in HDA, the 5th water molecule, contrary to ice and LDA, is positioned between the first and the second coordination shell. We found no convincing evidence of structural heterogeneity or regions referring to structural transition.

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.

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⁻⁶ 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.

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.

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

We studied the vapor-liquid coexistence region of seven 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)10.1063/1.1805516], the SWM4-DP [G. Lamoureux, A. D. MacKerell, Jr., and B. Roux, J. Chem. Phys. 119, 5185 (2003)10.1063/1.1598191], the SWM4-NDP [G. Lamoureux, E. Harder, I. V. Vorobyov, B. Roux, and A. D. MacKerell, Jr., Chem. Phys. Lett. 418, 245 (2006)10.1016/j.cplett.2005.10.135], 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)10.1063/1.3692602] 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)10.1063/1.3670962], 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, 194102 (2012)10.1063/1.4767063]. We calculated the second virial coefficient, the heat of vaporization, equilibrium vapor pressure, the vapor-liquid coexistence curve, and the surface tension in terms of the temperature. We determined and compared the critical temperatures, densities, and pressures of the models. We concluded that the high temperature slope of the (T-ρ) curve accurately predicts the critical temperature. We found that Gaussian charge distributions have clear advantages over the point charges describing the critical region. It is impossible to describe the vapor-liquid coexistence properties consistently with nonpolarizable models, even if their critical temperature is correct.

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)10.1063/1.1805516], the COS/B2 [H. Yu, T. Hansson, and W. F. van Gunsteren, J. Chem. Phys. 118, 221 (2003)10.1063/1.1523915], the SWM4-DP [G. Lamoureux, A. D. MacKerell, Jr., and B. Roux, J. Chem. Phys. 119, 5185 (2003)10.1063/1.1598191], the SWM4-NDP [G. Lamoureux, E. Harder, I. V. Vorobyov, B. Roux, and A. D. MacKerell, Jr., Chem. Phys. Lett. 418, 245 (2006)10.1016/j.cplett.2005.10.135], 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)10.1063/1.3692602] 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)10.1063/1.3670962], 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)10.1063/1.4746419]. 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.