Incorporating Phase-Dependent Polarizability in Non-Additive Electrostatic Models for Molecular Dynamics Simulations of the Aqueous Liquid-Vapor Interface

Abstract

We discuss a new classical water force field that explicitly accounts for differences in polarizability between liquid and vapor phases. The TIP4P-QDP (4-point transferable intermolecular potential with charge dependent-polarizability) force field is a modification of the original TIP4P-FQ fluctuating charge water force field of Rick et al.(1) that self-consistently adjusts its atomic hardness parameters via a scaling function dependent on the M-site charge. The electronegativity (χ) parameters are also scaled in order to reproduce condensed-phase dipole moments of comparable magnitude to TIP4P-FQ. TIP4P-QDP is parameterized to reproduce experimental gas-phase and select condensed-phase properties. The TIP4P-QDP water model possesses a gas phase polarizability of 1.40 Å(3) and gas-phase dipole moment of 1.85 Debye, in excellent agreement with experiment and high-level ab initio predictions. The liquid density of TIP4P-QDP is 0.9954(±0.0002) g/cm(3) at 298 K and 1 atmosphere, and the enthalpy of vaporization is 10.55(±0.12) kcal/mol. Other condensed-phase properties such as the isobaric heat capacity, isothermal compressibility, and diffusion constant are also calculated within reasonable accuracy of experiment and consistent with predictions of other current state-of-the-art water force fields. The average molecular dipole moment of TIP4P-QDP in the condensed phase is 2.641(±0.001) Debye, approximately 0.02 Debye higher than TIP4P-FQ and within the range of values currently surmised for the bulk liquid. The dielectric constant, ε = 85.8 ± 1.0, is 10% higher than experiment. This is reasoned to be due to the increase in the condensed phase dipole moment over TIP4P-FQ, which estimates ε remarkably well. Radial distribution functions for TIP4P-QDP and TIP4P-FQ show similar features, with TIP4P-QDP showing slightly reduced peak heights and subtle shifts towards larger distance interactions. Since the greatest effects of the phase-dependent polarizability are anticipated in regions with both liquid and vapor character, interfacial simulations of TIP4P-QDP were performed and compared to TIP4P-FQ, a static polarizability analog. Despite similar features in density profiles such as the position of the GDS and interfacial width, enhanced dipole moments are observed for the TIP4P-QDP interface and onset of the vapor phase. Water orientational profiles show an increased preference (over TIP4P-FQ) in the orientation of the permanent dipole vector of the molecule within the interface; an enhanced z-induced dipole moment directly results from this preference. Hydrogen bond formation is lower, on average, in the bulk for TIP4P-QDP than TIP4P-FQ. However, the average number of hydrogen bonds formed by TIP4P-QDP in the interface exceeds that of TIP4P-FQ, and observed hydrogen bond networks extend further into the gaseous region. The TIP4P-QDP interfacial potential, calculated to be -11.98(±0.08) kcal/mol, is less favorable than that for TIP4P-FQ by approximately 2% as a result of a diminished quadrupole contribution. Surface tension is calculated within a 1.3% reduction from the experimental value. Results reported demonstrate TIP4P-QDP as a model comparable to the popular TIP4P-FQ while accounting for a physical effect previously neglected by other water models. Further refinements to this model, as well as future applications are discussed.

Figure 1

(a) Molecular polarizability (in ų) as a function of Q_M. (b) Distribution of molecular polarizabilities within the condensed phase. Polarizabiities were calculated using Equation 7 and the charges from simulation.

Figure 2

Dipole moment distributions for TIP4P-QDP, TIP4P-FQ, and QDP-P1.

Figure 3

Radial distribution functions for TIP4P-QDP, TIP4P-FQ, and QDP-P1.

Figure 4

Interfacial profiles as a function of z-position relative to the center of mass for (a) the density and (b) the dipole moment of TIP4P-QDP, TIP4P-FQ, and QDP-P1. A 0.2 g/cm³ offset was applied to the density profiles to distinguish unique features.

Figure 5

Hydrogen bond profile. (a) The average number of hydrogen bonds as a function of z-position for TIP4P-QDP, TIP4P-FQ, and QDP-P1. (b) The probability of hydrogen bond formation as a function of z-position as calculated from the ratio of hydrogen bonds to coordination number. A definition of O–O distance less than 3.5 Å and H–O—O angle less than 30° were used as the hydrogen bond criteria.

Figure 6

(a) Profile of 〈cos θ〉, where θ is angle formed between the permanent dipole vector of the water molecule and the fixed z-axis. (b) Profile of 〈cos ϕ〉, where ϕ is the angle formed between the molecular plane of water and the fixed z-axis. (c) Profile of the z-induced dipole moment. All z-values are relative to the center of mass of the system. The GDS is represented on each plot as a dashed vertical line.

Figure 7

Total interfacial potential as a function of z relative to the center of mass of the water slab.

Figure 8

Contribution of each term in equation A18 to the total scaling function, g(Q_M). The top panel features the first term of equation A18 (solid line) as it compares to the total scaling function (dashed line) as a function of Q_M. The lower panel features the second (solid line) and third (dashed line) terms of equation A18. Terms 2 and 3 are approximately equal in magnitude while opposite in sign, which results in the first term’s dominance in g(Q_M).