Thermodynamic anomalies in a lattice model of water

We investigate a lattice-fluid model of water, defined on a three-dimensional body centered cubic lattice. Model molecules possess a tetrahedral symmetry, with four equivalent bonding arms, aiming to mimic the formation of hydrogen bonds. The model is similar to the one proposed by Roberts and Debenedetti [J. Chem. Phys. 105, 658 (1996)], simplified in that no distinction between bond "donors" and "acceptors" is imposed. Bond formation depends both on orientation and local density. In the ground state, we show that two different ordered (ice) phases are allowed. At finite temperature, we analyze homogeneous phases only, working out phase diagram, response functions, the temperature of maximum density locus, and the Kauzmann line. We make use of a generalized first-order approximation on a tetrahedral cluster. In the liquid phase, the model exhibits several anomalous properties observed in real water. In the low temperature region (supercooled liquid), there are evidences of a second critical point and, for some range of parameter values, this scenario is compatible with the existence of a reentrant spinodal.     

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

We present a new parametrization of the flexible, polarizable Thole-type model for water [J. Chem. Phys. 116, 5115 (2002); J. Phys. Chem. A 110, 4100 (2006)], with emphasis in describing the vibrational spectra of both water clusters and liquid water. The new model is able to produce results of similar quality with the previous versions for the structures and energetics of water clusters as well as structural and thermodynamic properties of liquid water evaluated with classical and converged quantum statistical mechanical atomistic simulations. At the same time it yields accurate redshifts for the OH vibrational stretches of both water clusters and liquid water.     

Polarizable and flexible model for ethanol

A polarizable, flexible model for ethanol is obtained based on an extensive series of B3LYP/6-311++G(d,p) calculations and molecular dynamics simulations. The ethanol model includes electric-field dependence in both the atomic charges and the intramolecular degrees of freedom. Field-dependent intramolecular potentials have been attempted only once previously, for OH and HH stretches in water [P. Cicu et al., J. Chem. Phys. 112, 8267 (2000)]. The torsional potential involving the hydrogen-bonding hydrogen in ethanol is found to be particularly field sensitive. The methodology for developing field-dependent potentials can be readily generalized to other molecules and is discussed in detail. Molecular dynamics simulations of bulk ethanol are performed and the results are assessed based on comparisons with the self-diffusion coefficient [N. Karger et al., J. Chem. Phys. 93, 3437 (1990)], dielectric constant [J. T. Kindt and C. A. Schmuttenmaer, J. Phys. Chem. 100, 10373 (1996)], enthalpy of vaporization [R. C. Wilhoit and B. J. Zwolinski, J. Phys. Chem. Ref. Data, Suppl. 2, 2 (1973)], and experimental interatomic distributions [C. J. Benmore and Y. L. Loh, J. Chem. Phys. 112, 5877 (2000)]. The simultaneous variation of the atomic charges and the intramolecular potentials requires modified equations of motion and a multiple time step algorithm has been implemented to solve these equations. The article concludes with a discussion of the bulk structure and properties with an emphasis on the hydrogen bonding network.     

Slow dynamics in a primitive tetrahedral network model

We report extensive Monte Carlo and event-driven molecular dynamics simulations of the fluid and liquid phase of a primitive model for silica recently introduced by Ford et al. [J. Chem. Phys. 121, 8415 (2004)]. We evaluate the isodiffusivity lines in the temperature-density plane to provide an indication of the shape of the glass transition line. Except for large densities, arrest is driven by the onset of the tetrahedral bonding pattern and the resulting dynamics is strong in Angell's classification scheme [J. Non-Cryst. Solids 131-133, 13 (1991)]. We compare structural and dynamic properties with corresponding results of two recently studied primitive models of network forming liquids-a primitive model for water and an angular-constraint-free model of four-coordinated particles-to pin down the role of the geometric constraints associated with bonding. Eventually we discuss the similarities between "glass" formation in network forming liquids and "gel" formation in colloidal dispersions of patchy particles.     

Computer simulation study of metastable ice VII and amorphous phases obtained by its melting

Molecular dynamics simulations of metastable ice VII and cubic ice Ic are carried out in order to examine (1) the ability of commonly used water interaction potentials to reproduce the properties of ices, and (2) the possibility of generating low-density amorphous (LDA) structures by heating ice VII, which is known to transform to LDA at approximately 135 K at normal pressure [S. Klotz, J. M. Besson, G. Hamel, R. J. Nelmes, J. S. Loveday, and W. G. Marshall, Nature (London) 398, 681 (1999)]. We test four simple empirical interaction potentials of water: TIP4P [W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983)], SPC/E [H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. B 91, 6269 (1987)], TIP5P [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 (2000)], and ST2 [F. H. Stillinger and A. Rahman, J. Chem. Phys. 60, 1545 (1974)]. We have found that TIP5P ice VII melts at 210 K, TIP4P at 90 K, and SPC/E at 70 K. Only TIP5P water after transition has a structure similar to that of LDA. TIP4P and SPC/E have almost identical structures, dissimilar to any known water or amorphous phases, but upon heating both slowly evolve towards LDA-like structure. ST2 ice VII is remarkably stable up to 430 K. TIP4P and SPC/E predict correctly the cubic ice collapse into a high-density amorphous ice (HDA) at approximately 1 GPa whereas TIP5P remains stable up to approximately 5 GPa. The densities of the simulated ice phases differ significantly, depending on the potential used, and are generally higher than experimental values. The importance of proper treatment of long-range electrostatic interactions is also discussed.     

Computational study of structural and dynamical properties of formamide-water mixtures

A molecular dynamics simulation study of structural and dynamical properties in liquid mixtures of formamide and water is presented. Site-site radial pair distribution functions, local mole fractions, pair energy distributions, and tetrahedral orientational order are the quantities analyzed to investigate the local structure in the simulated mixtures, along with a review of the intermolecular structure in terms of the distribution of hydrogen bonds. Our results indicate that there is a substitution of formamide molecules by water in the hydrogen bonds and a formation of a common hydrogen bond network. By analyzing the extent of tetrahedral order in the liquid as a function of composition, it is observed that whereas the tetrahedral network of liquid water is progressively lost by increasing the formamide concentration, the water structure within the first coordination shell is preserved and somewhat enhanced. The hydrogen-bond mean lifetimes were estimated by performing a time integration of the autocorrelation functions of bond occupation numbers. The lifetimes associated with hydrogen bonds between water, formamide, and interspecies pairs are found to increase with increasing formamide concentration. The lifetimes of the water hydrogen bonds show the largest variations, supporting the picture of an enhancement of the water structure among the nearest neighbors within the first coordination shell. We have used two different force field models for water, SPC/E [J. C. Berendsen et al., J. Phys. Chem. 91, 6269 (1987)] and TIP4P/2005 [J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005)]. Our results for structural and dynamical properties yield very small differences between those models, the TIP4P/2005 predicting a slightly more structured liquid and, consequently, exhibiting a slightly slower translational and librational dynamics.     

Thermodynamic properties in the molecular dynamics ensemble applied to the gaussian core model fluid

The thermodynamic properties of pressure, energy, isothermal pressure coefficient, thermal expansion coefficient, isothermal and adiabatic compressibilities, isobaric and isochoric heat capacities, Joule-Thomson coefficient, and speed of sound are considered in a classical molecular dynamics ensemble. These properties were obtained using the treatment of Lustig [J. Chem. Phys. 100, 3048 (1994)] and Meier and Kabelac [J. Chem. Phys. 124, 064104 (2006)], whereby thermodynamic state variables are expressible in terms of phase-space functions determined directly from molecular dynamics simulations. The complete thermodynamic information about an equilibrium system can be obtained from this general formalism. We apply this method to the gaussian core model fluid because the complex phase behavior of this simple model provides a severe test for this treatment. Waterlike and other anomalies are observed for some of the thermodynamic properties of the gaussian core model fluid.     

Cluster-variation approximation for a network-forming lattice-fluid model

We consider a three-dimensional lattice model of a network-forming fluid, which has been recently investigated by Girardi et al. by means of Monte Carlo simulations [J. Chem. Phys. 126, 064503 (2007)], with the aim of describing water anomalies. We develop an approximate semianalytical calculation, based on a cluster-variation technique, which turns out to reproduce almost quantitatively different thermodynamic properties and phase transitions determined by the Monte Carlo method. Nevertheless, our calculation points out the existence of two different phases characterized by long-range orientational order, and of critical transitions between them and to a high-temperature orientationally disordered phase. Also, the existence of such critical lines allows us to explain certain "kinks" in the isotherms and isobars determined by the Monte Carlo analysis. The picture of the phase diagram becomes much more complex and richer, though unfortunately less suitable to describe real water.     

Low-lying isomers and finite temperature behavior of (H2O)6 -

(H2O)(6) (-) appears as a "magic" number water cluster in (H2O)(n) (-) mass spectra. The structure of the (H2O)(6) (-) isomer dominating the experimental population has been established only recently [N. I. Hammer et al., J. Phys. Chem. A 109, 7896 (2005)], and the most noteworthy characteristic of this isomer is the localization of the excess electron in the vicinity of a double-acceptor monomer. In the present work, we use a quantum Drude model to characterize the low-energy isomers and the finite temperature properties of (H2O)(6) (-). Comparison with ab initio calculations shows that the use of a water model employing distributed polarizabilities and distributed repulsive sites is necessary to correctly reproduce the energy ordering of the low-lying isomers. Both the simulations and the ab initio calculations predict that there are several isomers of (H2O)(6) (-) significantly lower in energy than the experimentally observed species, suggesting that the experimental distribution is far from equilibrium.     

Analysis of experimental data for the nucleation rate of water droplets

A formula for the stationary nucleation rate J is proposed and used for analysis of experimental data for the dependence of J on the supersaturation ratio S in isothermal homogeneous nucleation of water droplets in vapors. It is found that the experimental data are described quite successfully by the proposed formula which is based on (i) the Gibbs presentation of the nucleation work in terms of overpressure, (ii) the Girshick-Chiu [J. Chem. Phys. 93, 1273 (1990); 94, 826 (1991)] self-consistency correction to the equilibrium cluster size distribution, and (iii) the Reguera-Rubi [J. Chem. Phys. 115, 7100 (2001)] kinetic accounting of the nucleus translational-rotational motion. The formula, like that of Wolk and Strey [J. Phys. Chem. B 105, 11683 (2001)], could be used as a semiempirical relation describing the J(S) dependence for nucleation in vapors of single-component droplets or crystals of substances with insufficiently well known equations of state.     

Diffusional anomaly and network dynamics in liquid silica

The present study applies the power spectral analysis technique to understand the diffusional anomaly in liquid silica, modeled using the Beest-Kramer-van Santen (BKS) potential. Molecular-dynamics simulations have been carried out to show that power spectrum of tagged particle potential energy of silica shows a regime with 1f(alpha) dependence on frequency f which is the characteristic signature of multiple time scale behaviour in networks. As demonstrated earlier in the case of water [J. Chem. Phys. 122, 104507 (2005)], the variations in the mobility associated with the diffusional anomaly are mirrored in the scaling exponent alpha associated with this multiple time scale behavior. Our results indicate that in the anomalous regime, as the local tetrahedral order decreases with temperature or pressure, the coupling of local modes to network reorganizations increases and so does the diffusivity. This symmetry-dependence of the vibrational couplings is responsible for the connection between the structural and diffusional anomalies.     

Local hardness equalization: exploiting the ambiguity

In the density-functional theory of chemical reactivity, the local hardness is known to be an ambiguous concept. The mathematical structure associated with this problematic situation is elaborated and three common definitions for the local hardness are critically examined: the frontier local hardness [S. K. Ghosh, Chem. Phys. Lett. 172, 77 (1990)], the total local hardness [S. K. Ghosh and M. Berkowitz, J. Chem. Phys. 83, 2976 (1985)], and the unconstrained local hardness [P. W. Ayers and R. G. Parr, J. Am. Chem. Soc. 122, 2010 (2000)]. The frontier local hardness has particularly nice properties: (a) it has smaller norm than most, if not all, other choices of the local hardness and (b) it is "unbiased" in an information-theoretic sense. For the ground electronic state of a molecular system, the frontier local hardness is equal to the global hardness. For an electronic system in its ground state, both the chemical potential and the frontier local hardness are equalized. The frontier local hardness equalization principle provides a computational approach for designing reagents with desirable chemical reactivity profiles.     

Communication: rationale for a new class of double-hybrid approximations in density-functional theory

We provide a rationale for a new class of double-hybrid approximations introduced by Bre?mond and Adamo [J. Chem. Phys. 135, 024106 (2011)] which combine an exchange-correlation density functional with Hartree-Fock exchange weighted by λ and second-order M?ller-Plesset (MP2) correlation weighted by λ(3). We show that this double-hybrid model can be understood in the context of the density-scaled double-hybrid model proposed by Sharkas et al. [J. Chem. Phys. 134, 064113 (2011)], as approximating the density-scaled correlation functional E(c)[n(1/λ)] by a linear function of λ, interpolating between MP2 at λ = 0 and a density-functional approximation at λ = 1. Numerical results obtained with the Perdew-Burke-Ernzerhof density functional confirms the relevance of this double-hybrid model.     

Contribution of water dimer absorption to the millimeter and far infrared atmospheric water continuum

We present a rigorous calculation of the contribution of water dimers to the absorption coefficient alpha(nu,T) in the millimeter and far infrared domains, over a wide range (276-310 K) of temperatures. This calculation relies on the explicit consideration of all possible transitions within the entire rovibrational bound state manifold of the dimer. The water dimer is described by the flexible 12-dimensional potential energy surface previously fitted to far IR transitions [C. Leforestier et al., J. Chem. Phys. 117, 8710 (2002)], and which was recently further validated by the good agreement obtained for the calculated equilibrium constant Kp(T) with experimental data [Y. Scribano et al., J. Phys. Chem. A. 110, 5411 (2006)]. Transition dipole matrix elements were computed between all rovibrational states up to an excitation energy of 750 cm(-1), and J=K=5 rotational quantum numbers. It was shown by explicit calculations that these matrix elements could be extrapolated to much higher J values (J=30). Transitions to vibrational states located higher in energy were obtained from interpolation of computed matrix elements between a set of initial states spanning the 0-750 cm(-1) range and all vibrational states up to the dissociation limit (approximately 1200 cm(-1)). We compare our calculations with available experimental measurements of the water continuum absorption in the considered range. It appears that water dimers account for an important fraction of the observed continuum absorption in the millimeter region (0-10 cm(-1)). As frequency increases, their relative contribution decreases, becoming small (approximately 3%) at the highest frequency considered nu=944 cm(-1).     

Nonadiabatic coupling vectors for excited states within time-dependent density functional theory in the Tamm-Dancoff approximation and beyond

Recently, we have proposed a scheme for the calculation of nonadiabatic couplings and nonadiabatic coupling vectors within linear response time-dependent density functional theory using a set of auxiliary many-electron wavefunctions [I. Tavernelli, E. Tapavicza, and U. Rothlisberger, J. Chem. Phys. 130, 124107 (2009)]. As demonstrated in a later work [I. Tavernelli, B. F. E. Curchod, and U. Rothlisberger, J. Chem. Phys. 131, 196101 (2009)], this approach is rigorous in the case of the calculation of nonadiabatic couplings between the ground state and any excited state. In this work, we extend this formalism to the case of coupling between pairs of singly excited states with the same spin multiplicity. After proving the correctness of our formalism using the electronic oscillator approach by Mukamel and co-workers [S. Tretiak and S. Mukamel, Chem. Rev. (Washington, D.C.) 102, 3171 (2002)], we tested the method on a model system, namely, protonated formaldimine, for which we computed S(1)/S(2) nonadiabatic coupling vectors and compared them with results from high level (MR-CISD) electronic structure calculations.     

Temperature dependence of homogeneous nucleation rates for water: near equivalence of the empirical fit of Wölk and Strey, and the scaled nucleation model

It is pointed out that the temperature fitting function of W?lk and Strey [J. Phys. Chem. 105, 11683 (2001)], recently shown to convert the Becker-D?ring [Ann. Phys. (Leipzig) 24, 719 (1935)] nucleation rate into an expression in agreement with much of the experimental water nucleation rate data, also converts the Becker-D?ring rate into a form nearly equivalent with the scaled nucleation rate model, J(scaled)=J(oc) exp[-16piOmega(3)(T(c)T-1)(3)3(ln S)(2)]. In the latter expression J(oc) is the inverse thermal wavelength cubed/sec, evaluated at T(c).     

A justification for a nonadiabatic surface hopping Herman-Kluk semiclassical initial value representation of the time evolution operator

A justification is given for the validity of a nonadiabatic surface hopping Herman-Kluk (HK) semiclassical initial value representation (SC-IVR) method. The method is based on a propagator that combines the single surface HK SC-IVR method [J. Chem. Phys. 84, 326 (1986)] and Herman's nonadiabatic semiclassical surface hopping theory [J. Chem. Phys. 103, 8081 (1995)], which was originally developed using the primitive semiclassical Van Vleck propagator. We show that the nonadiabatic HK SC-IVR propagator satisfies the time-dependent Schrodinger equation to the first order of variant Planck's over 2pi and the error is O(variant Planck's over 2pi(2)). As a required lemma, we show that the stationary phase approximation, under current assumptions, has an error term variant Planck's over 2pi(1) order higher than the leading term. Our derivation suggests some changes to the previous development, and it is shown that the numerical accuracy in applications to Tully's three model systems in low energies is improved.     

Monte Carlo simulations of critical cluster sizes and nucleation rates of water

We have calculated the critical cluster sizes and homogeneous nucleation rates of water at temperatures and vapor densities corresponding to experiments by Wolk and Strey [J. Phys. Chem B 105, 11683 (2001)]. The calculations have been done with an expanded version of a Monte Carlo method originally developed by Vehkamaki and Ford [J. Chem. Phys. 112, 4193 (2000)]. Their method calculates the statistical growth and decay probabilities of molecular clusters. We have derived a connection between these probabilities and kinetic condensation and evaporation rates, and introduce a new way for the calculation of the work of formation of clusters. Three different interaction potential models of water have been used in the simulations. These include the unpolarizable SPC/E [J. Phys. Chem. 91, 6269 (1987)] and TIP4P [J. Chem. Phys. 79, 926 (1983)] models and a polarizable model by Guillot and Guissani [J. Chem. Phys. 114, 6720 (2001)]. We show that TIP4P produces critical cluster sizes and a temperature and vapor density dependence for the nucleation rate that agree well with the experimental data, although the magnitude of nucleation rate is constantly overestimated by a factor of 2 x 10(4). Guissani and Guillot's model is somewhat less successful, but both the TIP4P and Guillot and Guissani models are able to reproduce a much better experimental temperature dependency of the nucleation rate than the classical nucleation theory. Using SPC/E results in dramatically too small critical clusters and high nucleation rates. The water models give different average binding energies for clusters. We show that stronger binding between cluster molecules suppresses the decay probability of a cluster, while the growth probability is not affected. This explains the differences in results from different water models.     

Water models based on a single potential energy surface and different molecular degrees of freedom

Up to now it has not been possible to neatly assess whether a deficient performance of a model is due to poor parametrization of the force field or the lack of inclusion of enough molecular properties. This work compares several molecular models in the framework of the same force field, which was designed to include many-body nonadditive effects: (a) a polarizable and flexible molecule with constraints that account for the quantal nature of the vibration [B. Hess, H. Saint-Martin, and H. J. C. Berendsen, J. Chem. Phys. 116, 9602 (2002), H. Saint-Martin, B. Hess, and H. J. C. Berendsen, J. Chem. Phys. 120, 11133 (2004)], (b) a polarizable and classically flexible molecule [H. Saint-Martin, J. Hernandez-Cobos, M. I. Bernal-Uruchurtu, I. Ortega-Blake, and H. J. C. Berendsen, J. Chem. Phys. 113, 10899 (2000)], (c) a polarizable and rigid molecule, and finally (d) a nonpolarizable and rigid molecule. The goal is to determine how significant the different molecular properties are. The results indicate that all factors--nonadditivity, polarizability, and intramolecular flexibility--are important. Still, approximations can be made in order to diminish the computational cost of the simulations with a small decrease in the accuracy of the predictions, provided that those approximations are counterbalanced by the proper inclusion of an effective molecular property, that is, an average molecular geometry or an average dipole. Hence instead of building an effective force field by parametrizing it in order to reproduce the properties of a specific phase, a building approach is proposed that is based on adequately restricting the molecular flexibility and/or polarizability of a model potential fitted to unimolecular properties, pair interactions, and many-body nonadditive contributions. In this manner, the same parental model can be used to simulate the same substance under a wide range of thermodynamic conditions. An additional advantage of this approach is that, as the force field improves by the quality of the molecular calculations, all levels of modeling can be improved.