Hydrothermal properties of the COS/D2 water model: a polarizable charge-on-spring water model, at elevated temperatures and pressures

Abstract

Molecular simulations have been conducted to assess the pVT properties and static permittivity of the charge-on-spring polarizable water model COS/D2 at hydrothermal conditions from 300 to 450 °C and bulk densities of 0.001 to 1.0 g cm⁻³. The results indicate that the model performs well for reproducing volumetric and dielectric properties of real water. The liquid–vapor coexistence curve and critical point of COS/D2 water (ca. 334 °C, 21.7 MPa and 0.32 g cm⁻³) are also in reasonable agreement with those of real water (ca. 374 °C, 22.1 MPa and 0.322 g cm⁻³). We include comparisons to other polarizable (as well as non-polarizable) water models where available. The good performance of the COS/D2 model at high temperatures and pressures is noteworthy considering that the model was parameterized exclusively at room temperature and pressure. Moreover, these results imply that COS/D2 is a suitable water model for simulating solutions of non-electrolyte and electrolyte solutes at elevated temperatures and pressures.

---

I. Introduction

Hydrothermal fluid properties – physical and chemical properties of aqueous fluids and solutions at elevated temperatures and pressures – are key parameters in many geochemical and industrial processes and applications. Molecular simulations of fluids at elevated temperatures and pressures have become an important method for investigating such properties.^1,2 Commonly, molecular simulations allow insights into fluid properties at conditions that are difficult (or impossible) to achieve experimentally, and provide one of the only means of analyzing some structural (e.g. solvation) properties of materials at such conditions. For example, molecular simulations have been extensively used to investigate the structural properties and hydrogen-bonding of pure H₂O up to temperatures of several hundred degrees centigrade, and from vapor-like to liquid-like densities.¹ Similarly, molecular simulations have been used to characterize solvation, speciation and thermodynamic properties of aqueous non-electrolyte and electrolyte solutes near the critical point of H₂O^3–8 and up to extreme temperatures and pressures exceeding 1000 °C and 2 GPa.^9,10

To enable quantitative assessment of fluid (solvent and solute) properties from molecular simulations, a suitable water model is a key requirement. Despite decades of advances in terms of parameterization of water models, still the most commonly used models for simulations at hydrothermal conditions are the SPC¹¹ and SPC/E¹² models (both of which are three-site, rigid, non-polarizable models), owing to the simplicity of their design, their existing force-fields compatible with many biomolecules and electrolytes,¹³ and (especially in the case of SPC/E) their overall good performance in terms of hydrothermal properties of water. The TIP family of models¹⁴ are also popular in molecular dynamics simulation of water and solutions, and numerous variants of this family have emerged in recent years, each with advantages and disadvantages in terms of specific properties and ranges of applicability.¹⁵

Several approaches have been suggested and applied to account for polarizability in molecular simulations.^16–26 Quantum-chemical calculations of molecular polarizability provide the most rigorous approach, but are computationally expensive and may be prohibitive for complex systems.²⁰ Amongst classical methods for representing molecular systems, induced molecular point dipoles^27–34 or induced atomic dipoles^35–40 have been used to model polarizability. The atom-dipole model of Silberstein²¹ and Applequist et al.¹⁷ is based on isotropic polarization of individual atoms, which are combined assuming additivity in order to obtain molecular polarizability; anisotropy of the molecular polarizability results from the geometric constraints when summing the polarizabilities of the constituent atoms, according to their positions in the molecule. However, several studies have noted that the latter method fails to account for effects of the local chemical environment on polarizability,^20,24,26 and alternative methods with explicit treatment of local environment and structure have thus been developed (e.g., polarizability of specific bond types). Thus, Jensen et al.²⁰ noted that the isotropic component of molecular polarizability is implicitly additive, whereas the full molecular polarizability can be treated as an additive sum of atomic polarizabilities only if the details of local environment are fully accounted for. The method of Thole²⁴ is to represent bonded atoms according to “smeared” dipoles rather than isotropic point dipoles. This latter approach is adopted in many recent studies, including models for H₂O.^18,19 Other models use fluctuating charges to mimic polarizability^41–44 or Drude oscillators^45–55 (=charge-on-spring), wherein the position of the point charge is allowed to move in response to the electrostatic environment. Other recently developed polarizable water models include the GCP (Gaussian charge polarizable) model,⁵⁶ the SWM models,^48,49 and the BK-series models.^55,57,58 A full overview of polarizability in molecular mechanics was recently presented by Lopes et al.⁵⁹

Several of these models have been tested to various extents for simulations at hydrothermal conditions, especially along the liquid–vapor coexistence curve, although few data are available for the performance of polarizable water models at single-phase liquid, vapor or supercritical fluids conditions. In the present study, we test a newly developed polarizable water model (COS/D2 (ref. 60)) over a wide range of pressure–temperature conditions, with the aim of assessing whether this new model will be suitable for simulating electrolyte solutions in hydrothermal environments relevant to industrial and geological systems. As described below, the COS/D2 model is a (damped) charge-on-spring model designed to improve upon the simulated properties of previous charge-on-spring models at ambient conditions;⁶⁰ however, as yet the COS/D2 model has not been tested under hydrothermal conditions. The latter tests are the subject of the present study.

Of the models listed above, the SPC/E and TIP4P models in particular have been shown to provide good quantitative agreement with the properties of water (e.g., the dielectric properties⁶¹ and critical point⁶²) at elevated temperatures and pressures, and consequently these models are still widely used to represent aqueous fluids at hydrothermal conditions.¹ In part, the continued prevalence of these water models in high-temperature simulations likely reflects the fact that few other models have been featured in such comprehensive force-field parameterizations. For example, in the case of the SPC/E model, the GROMOS force field (used in the present study) contains numerous biomolecular residues, non-electrolytes and ions, allowing simulations of complex systems including dilute electrolytes as well as concentrated brines,¹³ which are of interest in industrial process applications as well as in geochemistry of Earth's interior. Moreover, Vega and Abascal⁶³ compared several (nonpolarizable) water models and found that the SPC/E model ranked second best (inferior only to the TIP4P/2005 model) in terms of reproducing the properties of H₂O over a wide range of physical conditions, despite the simplicity of the SPC/E model. Hence, in the present study we use the SPC/E model as a benchmark, for comparison of the newly developed COS/D2 model (described below).

Owing to the rigid, non-polarizable nature of the SPC and SPC/E water models, these models may not be appropriate for simulating some scenarios at high-temperatures (as well as ambient conditions), in which polarization plays a role.⁶³ For example, electronic polarization likely plays a key role in concentrated electrolyte solutions, owing to the inhomogeneous environment of H₂O molecules such fluids. In addition, neglect of molecular polarizability imposes several limitations in simulating H₂O at extreme conditions (i.e., high temperatures and/or pressures⁶³). These limitations may be compounded in the case of electrolyte solutions in which solutes contain a high charge density, such as divalent Ca²⁺ and/or SO₄²⁻. Consequently, ion association and clustering in solutions containing divalent ions can be dramatically overestimated even at ambient conditions, if a non-polarizable water model is used;⁶⁴ this problem can be overcome by using a polarizable water model for the solvent.⁶⁴ Thus, polarizable water models, coupled with ion force-field parameters optimized for the respective solvent, would be advantageous for simulating ion solvation and ion association/clustering (as well as other phenomena) at hydrothermal conditions. Simulations of concentrated electrolyte solutions using a polarizable water model at hydrothermal conditions will thus require (i) rigorous testing and detailed characterization of the polarizable solvent at elevated temperatures and pressures, and (ii) ion force-field parameterization optimized for the specific solvent. The latter of these (ion parameterization) itself requires knowledge of the solvent pressure–volume–temperature (pVT) and dielectric properties.^65,66 Therefore, in the present study we focus on obtaining a detailed characterization of the hydrothermal properties of one polarizable water model.

In recent years, several polarizable water models have emerged based on various approaches to simulating molecular polarizability.^56,57,60 In the present study, we focus on a newly developed polarizable model based on the “charge-on-spring” (Drude oscillator) design and termed the “COS” family of models.⁵¹ The COS water models are composed of a rigid three-point frame representing the atomic positions, similar to the SPC family of models, with the addition of a flexible point charge (the charge on spring = COS) virtual site. The three-point frame consists of atomic positions of hydrogen and oxygen, which are treated in the normal way (partial point charges located at the atomic center, Lennard-Jones interaction parameters). The additional flexible point charge allows polarization of the otherwise rigid molecule by displacing the point-charge position according to the electrostatic environment of the molecule. The COS formulation balances the computational efficiency of rigid models with added functionality of solvent polarizability. Several COS-type models have been developed.⁵¹ Amongst these, the COS/D and COS/D2 models employ a non-linear dependence of the induced dipole on the electric field for higher field strengths, which mimics the electronic polarization of water molecules and the saturation of the polarization at high field strengths.⁶⁰ In the case of COS/D2, the non-linear dependence of the induced dipole on the electric field at high field strengths provides improved performance compared to models employing a linear dependence (e.g., COS/G2 (ref. 60)). In contrast to earlier charge-on-spring-type water models, the COS/D2 model also includes a Lennard-Jones repulsion potential at the hydrogen atoms, to prevent excessive hydrogen bonding, which had resulted in unrealistically low diffusivity for the earlier COS/D model.⁶⁰ Recent tests involving several COS-type water models⁶⁰ showed that the COS/D2 model quantitatively reproduces the thermodynamic and transport properties of real water at ambient temperature and pressure. However, the reliability of this model (which was parameterized strictly based on properties at ambient temperature and pressure) at hydrothermal conditions is unknown. In the present study, we investigate the application of the COS/D2 water model⁶⁰ at temperatures from 300 to 450 °C and from vapor-like to liquid-like bulk densities (0.001 to 1.0 g cm⁻³). The pVT properties, static dielectric properties, and fluid structure of COS/D2 water at these conditions provide an assessment of the suitability of the COS/D2 model for hydrothermal simulations.

II. Methods

Molecular simulations were conducted using the GROMOS simulation package.^67–69 The cubic boxes contained 4500 COS/D2 molecules, and periodic boundary conditions were applied. Solvent bond lengths were constrained with the SHAKE algorithm⁷⁰ with a relative geometric tolerance of 10⁻⁴. The SHAKE algorithm⁷⁰ was used (rather than the alternative LINCS algorithm⁷¹) in order to maintain consistency with the original study on the COS/D2 model.⁶⁰ Simulations were conducted at constant temperature using weak coupling⁷² to a single heat bath with a coupling constant of 0.1 ps, and at constant density (i.e., in the NVT ensemble). The translational center of mass motion was subtracted each 1000 time steps. Electrostatic interactions were handled using a triple-range cut-off scheme. This cut-off scheme consisted of an inner cut-off radius of 0.8 nm, within which all atomic interactions were updated every time step. Between the inner cutoff radius (0.8 nm) and the outer cutoff radius of 1.4 nm, atomic interactions were update every five time steps. Outside of the outer cutoff radius (>1.4 nm), atomic interactions were not explicitly calculated, and were approximated using a Barker–Watts reaction-field.⁷³ For each simulated state point, the system was equilibrated at the given temperature and density for 100 ps, and the next 2 ns of simulation were used for analysis. For calculations of the solvent dielectric constant (described below), each 2 ns simulation was extended to 6–8 ns (depending on the time required for convergence of the static permittivity).

Initially (for equilibration and the first 2 ns of simulation), a reaction field dielectric permittivity⁷³ (accounting for electrostatic interactions beyond a cutoff distance of 1.4 nm) was set according to the IAPWS formulation.⁷⁴ To determine the static permittivity ε(0) of the solvent at the simulation state points, we used the fluctuation of the total box dipole moment M⃑ according to:⁷⁵

where ε_rf is the relative dielectric permittivity of the reaction field continuum, V is the volume of the box, k_B is the Boltzmann constant, T is the absolute temperature, and ε₀ is the dielectric permittivity of free space. To ensure satisfactory convergence of the calculated relative dielectric constant, the initial input reaction field permittivity was adjusted according to the relative dielectric constant obtained from the analysis of the simulation, and the simulation was then continued for an additional 2 ns. This process was repeated (iterating as required), until the static dielectric constant of the COS/D2 solvent was determined to within less than one unit in terms of relative permittivity. In general the static permittivity could be constrained within 3 and 4 iterations, such that the total simulation time for each state point was between 6 and 8 ns.

The pVT properties of the COS/D2 fluid were determined from the NVT simulations according to the average pressure at each given T and density. Because the simulations were done in the NVT ensemble, the simulated state points at varying temperatures represent the isochoric properties of the fluid, whereas simulated state points at varying density represent the isothermal properties. Isochoric and isothermal pressure variations were used to calculate the derivatives ∂P/∂T|_V and ∂P/∂V|_T, and ∂T/∂V|_P was assessed from these using the chain rule. These derivatives were used to assess the isothermal compressibility κ_T and isobaric thermal expansivity α_P of the fluid, according to the equations

andrespectively. The pVT data were also used to extract second virial coefficients, according to the limit

Liquid–vapor coexistence was determined by applying the Maxwell equal area construction⁷⁶ to the P–V data along each simulated isotherm. The resulting liquid and vapor densities were then used to estimate the critical temperature and density.⁶²

For comparisons between COS/D2 and experiment and an existing non-polarizable water model, the same properties were calculated at the equivalent state points using the IAPWS formulations^74,77 for real water, and additional GROMOS simulations were conducted at the same conditions using the SPC/E water model. Simulations using SPC/E H₂O¹² were conducted using the 54A8 force field.¹³ Real water pVT properties were determined using the IAPWS-95 equation of state,⁷⁷ and static relative dielectric constants were calculated using the IAPWS formulation.⁷⁴ We also include comparisons to available literature data. This includes recent comprehensive data on hydrothermal pVT properties of TIP4P/2005,⁷⁸ as well as data on liquid–vapor equilibria and critical properties for numerous water models. Unfortunately, for many water models, few data are available throughout ranges of hydrothermal pressure–temperature conditions. Recent studies have reported data on hydrothermal properties of polarizable water models,^79,80 although not to such wide ranges of pressure and temperature as explored in the present study (relevant to geologic and industrial processes).

Simulations were conducted at temperatures of 300 to 450 °C at 25 °C increments, at fluid densities from 0.001 to 1 kg m⁻³ (0.001 to 1 g cm⁻³), for both COS/D2 and SPC/E H₂O (Table 1).

Table 1 pVT and dielectric properties of COS/D2 H₂O, and comparison with SPC/E H₂O and experimental values^ab

ρ/g cm^−3	COS/D2				SPC/E				Exp.^c
	P/MPa	κT/MPa^−1	αP/°C^−1	ε(0)	P/MPa	κT/MPa^−1	αP/°C^−1	ε(0)	P/MPa	κT/MPa^−1	αP/°C^−1	ε(0)
a Entries marked by a dash (—) correspond to pVT points characterized by two-phase (liquid + vapor) coexistence.b The complete dataset is available in the ESI.c Calculated according to the IAPWS equations.^74,77
350 °C (≈Tc SPC/E)
0.001	0.3	3.4	6.9 × 10^−4	1	0.3	4.3	1.9 × 10^−3	1	0.2	6.0	1.6 × 10^−3	1
0.05	11.3	1.2 × 10^−1	3.8 × 10^−3	1.3	8.9	1.7 × 10^−1	6.5 × 10^−3	1.5	11.1	1.3 × 10^−1	4.7 × 10^−3	1.3
0.1	18.1	1.0 × 10^−1	8.1 × 10^−3	1.8	12.5	2.0 × 10^−1	1.6 × 10^−2	2.2	16.1	1.8 × 10^−1	1.6 × 10^−2	1.8
0.2	23.8	1.6 × 10^−1	2.8 × 10^−2	2.9	14.4	1.4	2.2 × 10^−1	3.8	16.7	—	—	—
0.3	25.6	2.3 × 10^−1	6.3 × 10^−2	4.3	14.5	—	—	5.8	16.7	—	—	—
0.4	27.6	1.2 × 10^−1	4.7 × 10^−2	6.2	15.1	1.8 × 10^−1	6.3 × 10^−2	8.0	16.7	—	—	—
0.5	31.9	2.6 × 10^−2	1.4 × 10^−2	8.5	19.2	2.0 × 10^−2	1.2 × 10^−2	11.2	16.7	—	—	—
0.6	46.7	6.7 × 10^−3	5.3 × 10^−3	10.9	35.5	6.1 × 10^−3	4.5 × 10^−3	13.8	20.6	9.9 × 10^−3	6.9 × 10^−3	13.9
0.7	87.4	2.3 × 10^−3	2.5 × 10^−3	13.3	79.0	2.2 × 10^−3	2.2 × 10^−3	17.3	54.6	2.6 × 10^−3	2.6 × 10^−3	17.6
0.8	183	9.5 × 10^−4	1.4 × 10^−3	18.0	179	9.4 × 10^−4	1.3 × 10^−3	19.8	142	9.9 × 10^−4	1.3 × 10^−3	21.7
0.9	365	4.6 × 10^−4	8.9 × 10^−4	19.8	355	4.8 × 10^−4	8.4 × 10^−4	22.9	320	4.7 × 10^−4	8.4 × 10^−4	25.9
1.0	683	2.5 × 10^−4	6.1 × 10^−4	25.2	649	2.8 × 10^−4	6.1 × 10^−4	26.3	628	2.6 × 10^−4	5.7 × 10^−4	30.6
375 °C (≈Tc real + 0.9 °C)
0.001	0.3	3.3	6.7 × 10^−4	1	0.3	4.1	1.8 × 10^−3	1	0.2	6.0	1.6 × 10^−3	1
0.05	12.2	1.0 × 10^−1	3.3 × 10^−3	1.3	9.8	1.4 × 10^−1	5.4 × 10^−3	1.5	11.8	1.1 × 10^−1	3.8 × 10^−3	1.3
0.1	20.1	8.1 × 10^−2	6.4 × 10^−3	1.8	14.5	1.3 × 10^−1	1.1 × 10^−2	2.1	17.9	1.2 × 10^−1	9.8 × 10^−3	1.8
0.2	28.4	8.7 × 10^−2	1.6 × 10^−2	2.9	18.5	2.3 × 10^−1	3.7 × 10^−2	3.7	22.0	4.2 × 10^−1	7.8 × 10^−2	3.1
0.3	32.7	8.5 × 10^−2	2.4 × 10^−2	4.4	20.6	1.8 × 10^−1	4.2 × 10^−2	5.7	22.3	7.0	1.9	4.9
0.4	37.4	4.7 × 10^−2	1.8 × 10^−2	6.5	24.1	5.1 × 10^−2	1.7 × 10^−2	7.9	22.5	5.6 × 10^−1	1.8 × 10^−1	7.2
0.5	45.7	1.6 × 10^−2	8.7 × 10^−3	8.6	32.5	1.5 × 10^−2	7.2 × 10^−3	10.3	24.7	3.7 × 10^−2	1.7 × 10^−2	10.1
0.6	66.3	5.2 × 10^−3	4.1 × 10^−3	11.2	54.6	4.9 × 10^−3	3.5 × 10^−3	13.2	37.7	7.0 × 10^−3	5.0 × 10^−3	13.3
0.7	116	2.0 × 10^−3	2.2 × 10^−3	14.0	108	1.9 × 10^−3	2.0 × 10^−3	17.2	79.2	2.2 × 10^−3	2.3 × 10^−3	16.9
0.8	221	8.8 × 10^−4	1.3 × 10^−3	16.8	212	8.9 × 10^−4	1.2 × 10^−3	19.2	177	9.2 × 10^−4	1.3 × 10^−3	20.7
0.9	414	4.3 × 10^−4	8.5 × 10^−4	20.4	398	4.6 × 10^−4	8.1 × 10^−4	22.3	364	4.5 × 10^−4	8.0 × 10^−4	24.7
1.0	744	2.4 × 10^−4	5.9 × 10^−4	23.1	705	2.6 × 10^−4	5.8 × 10^−4	24.5	684	2.5 × 10^−4	5.5 × 10^−4	29.2

---

Finally, one set of simulations was done as an initial test of ion solvation in COS/D2. The test simulations involved a single Na⁺ ion dissolved in a system of 4500 COS/D2 (or SPC/E), with a bulk density of 100 kg m⁻³ (0.1 g cm⁻³) at 375 °C. Owing to the reconnaissance nature of these latter tests, no attempt was made at this stage to optimize the ion parameters specifically for the COS/D2 solvent. Thus, the Na⁺ parameters used here are those optimized for the SPC/E solvent from Reif and Hünenberger.⁶⁵ Structural solvation properties were analyzed by computing the radial distribution functions and running coordination numbers of the ion–solvent (Na⁺–O) pair correlation.

III. Results and discussion

A first noteworthy comparison between simulations of COS/D2 and SPC/E is in terms of computational time required to run equivalent simulations. Naturally, because of the greater complexity of the COS/D2 model (three point charges, the charge-on-spring, and LJ-interaction parameters for the hydrogen atoms) compared to the SPC/E model (only three point charges; no charge-on-spring; LJ interaction for the oxygen atom only), the COS/D2 simulations are somewhat more computational intensive. The specific computational times will depend on the particular setup used, but we found that a rule-of-thumb for estimating computational times is that for equivalent systems (number of molecules, bulk densities, etc.) a simulation using COS/D2 takes approximately 4× the time required for an analogous SPC/E simulation. Thus, a simulation requiring ∼10 hours using SPC/E would take ∼40 hours for COS/D2, etc.

The pVT properties and related derivatives, as well as the relative static dielectric constants of COS/D2 water, as well as SPC/E and real water, at the state points investigated in the present study are listed in Table 1. Table 2 lists the approximate pVT conditions of the critical point for these two H₂O models and real H₂O, as well as for several other water models reported from literature.

Table 2 Critical properties of water models

	Tc/°C	Pc^b/MPa	ρc/g cm^−3
a Calculated according to the IAPWS-95 EOS.^77b Blank entries = no data available in literature.
COS/D2	334	21.7	0.32
SPC/E	350	14.5	0.30
Exp.^a	374.1	22.1	0.322
SPC^83,88	314	12.6	0.27
MSPC/E^89	337	—	0.287
SPC/P^34	278	—	0.34
SCPDP^34	265	—	0.32
KJ^34	412	—	0.34
PPC^90	333	—	0.30
TIP3P^63	305	12.6	0.27
TIP4P^63	315	14.9	0.31
TIP4P-2005 ^63	367	14.6	0.31
TIP5P^63	248	8.6	0.37
TIP4P-FQ^91	297	—	0.30
TIP4P/P^34	314	—	0.35
BK3 (ref. 57)	356	23.3	0.33
BKd3 (ref. 55)	361	24.6	0.32
GCPM^56	369	21.4	0.33

---

Volumetric properties of the H₂O are key parameters for simulating hydrothermal systems. In addition, for simulating concentrated electrolyte solutions, ion parameterization based on solvation free energies require knowledge of the solvent pVT properties as input parameters.^65,66 Fig. 1 shows H₂O isochores projected in pressure–temperature space, with the upper panels showing the full temperature–pressure range of this study (up to 1000 MPa = 10 kbar) and the lower panels focusing on the liquid–vapor coexistence curves. Isochoric pT properties are important for matching simulated properties with results of hydrothermal experiments, because many such experiments are conducted under isochoric conditions (e.g., synthetic fluid inclusions,⁸¹ or hydrothermal diamond anvil cells⁸²). Unfortunately, although isochoric properties are important for simulations designed to investigate experimental results, few published data on isochoric properties of other water models are available in literature for comparison. Fig. 1 shows comparison with data for the TIP4P/2005 model from literature.⁷⁸ The TIP4P/2005 model performs quite well for reproducing the experimental pVT trends, and at high temperatures and pressures the results from SPC/E and TIP4P/2005 appear fairly similar (Fig. 1, upper right panel), although in detail near the liquid–vapor coexistence curve the differences between these models become apparent (Fig. 1, lower right panel). In general, both SPC/E and COS/D2 models perform fairly well in terms of qualitative and quantitative reproduction of the H₂O isochoric pT properties (Fig. 1 and Table 1). At lower pressures, in the vicinity of the liquid–vapor coexistence curve, the COS/D2 model generally performs somewhat better than the SPC/E model (which underestimates pressures, including the vapor pressure), but at higher pressures the SPC/E model outperforms the COS/D2 model (which tends to overestimate pressure at a given temperature and density). For both models, absolute difference with respect to experimental isochoric pressures generally increase with increasing temperature to 450 °C, owing to the divergence of the isochores with increasing p and T. Nevertheless, the correspondence between model and experimental isochoric pressures is good – the maximum difference between the COS/D2 model versus real H₂O occurs at 1 g cm⁻³ and 450 °C, at which the COS/D2 model has a pressure of 927 MPa, versus a pressure of 855 MPa for real H₂O⁷⁷ (Fig. 1 and Table 1).

Fig. 1 Pressure–temperature projection showing isochores of real (left), COS/D2 (middle) and SPC/E plus TIP4P/2005 (right) H₂O. Lower panels show expanded view of the region proximal to the liquid–vapor curve, corresponding to the grey-shaded areas in the upper panels. The TIP4P/2005 data are from literature.⁷⁸

Fig. 2 shows the pVT properties of real and simulated H₂O projected in pressure–density space and showing isothermal properties as well as the liquid–vapor coexistence. Both the COS/D2 and SPC/E models recreate qualitatively the trends in isothermal pV properties, exhibiting the classical cubic form. Quantitative agreement between the models and real H₂O is fair. Fig. 2 also shows data for the TIP4P/2005 model from literature.⁷⁸ Note that the 452 °C isotherm for TIP4P/2005 nearly overlaps with the 425 °C isotherm for SPC/E, while the 352 °C isotherm for TIP4P/2005 nearly overlaps with the 325 °C isotherm for SPC/E (Fig. 2; right). Thus, these data indicate that at a given pressure and density, the temperature of TIP4P/2005 H₂O is approximately 25 °C lower than for SPC/E (also see Fig. 1). Fig. 3 shows a close-up view of the liquid–vapor coexistence curve from 250 °C to the critical point, for various water models. Fig. 1 also shows the liquid–vapor coexistence curves of real and modeled H₂O according to the IAPWS-95 equation of state⁷⁷ and the present simulations (see the lower panels of Fig. 1). For both real and modeled H₂O, the liquid–vapor coexistence curve terminates at the critical point of the fluid, each of which is listed in Table 2. Table 2 also lists published values of critical properties of other selected water models. The critical point of real H₂O occurs at 374.1 °C, 22 MPa and 0.322 g cm⁻³.⁷⁷ The critical point of SPC/E H₂O has been variously reported at either 350 or 367 °C.^2,62 In the present study, we find the critical point of SPC/E H₂O at 350 °C, 14.5 MPa and 0.3 g cm⁻³, thus the critical temperature of this model is at the lower end of the range reported in literature.² For the polarizable COS/D2 model, the critical point occurs at 334 °C, 21.7 MPa and 0.32 g cm⁻³ (Fig. 3 and Table 2). The critical pressure and density of COS/D2 are in good agreement with those of real H₂O, but the critical temperature underestimated by about 40 °C. Even so, the COS/D2 model provides a better predicted critical temperature than some common H₂O models including SPC H₂O (T_c = 314 °C),⁸³ TIP3P (T_c = 305 °C), TIP4P (T_c = 315 °C)¹⁵ and TIP5P (T_c = 248 °C) (Table 2). However, several other models yield significantly better estimates of the critical temperature, including the TIP4P/2005 (T_c = 367 °C), BK3 and BKd3 (T_c = 356 and 361 °C, respectively), and GCPM (T_c = 369 °C) models (Fig. 3 and Table 2). The COS/D2 also yields one of the better estimates for the critical pressure (21.7 MPa, compared to 22.1 MPa for real water), which is generally underestimated by the SPC/E model and TIP family of models (although is in good agreement with both the BK3/BKd3 and GCPM models; Table 2), and a good estimate for the critical density (0.32 g cm⁻³, compared to 0.322 g cm⁻³ for real water), which is generally in good agreement amongst most of the models listed in Table 2 (with the exceptions of TIP3P and TIP5P). Consequently, the COS/D2 model also yields fair reproduction of the vapor pressures and saturated fluid densities along the liquid–vapor curve from 300 °C to the critical point (Fig. 1–3). Note that the BK3 and BKd3 models, TIP4P/2005 model, and the GCPM model all yield excellent reproduction of the vapor pressure and liquid–vapor curve properties for real water, as well (Fig. 3).^15,56,57,63 In summary, in terms of liquid–vapor equilibria and critical properties, the COS/D2 model is generally around the middle amongst the water models included here, performing somewhat better than some models, but not as well as some others. It should be reemphasized that the fairly good performance of the COS/D2 model in terms of critical properties and liquid–vapor equilibria is despite the fact that the model was parameterized only according to conditions at ambient temperature and pressure.⁶⁰

Fig. 2 Pressure–density projection showing isotherms of real (left), COS/D2 (middle) and SPC/E plus TIP4P/2005 (right) H₂O. Liquid–vapor coexistence is represented by thick black curves. For COS/D2 and SPC/E H₂O, simulated state points are shown as open square symbols, and the metastable/unstable parts of the isotherms (within the liquid–vapor field) are colored grey. The TIP4P/2005 data are from literature.⁷⁸

Fig. 3 Liquid–vapor coexistence curves of several water models. The COS/D2 and SPC/E data are from the present study, whereas data for the other models are from literature.^{55,56,63,83,88,92,93} Note that the TIP5P liquid–vapor curve is not shown because the critical point for this model lies at 248 °C (i.e., below the lowest temperature shown; Table 2).

In addition to pVT properties, we investigated the derivative properties of isothermal compressibility and isobaric thermal expansivity (Table 1 and Fig. 4). Fig. 4 shows the variation of these properties with fluid density, along isotherms representing both the estimated critical temperature, and the estimate critical temperature plus 1 °C, for each of the COS/D2 and SPC/E models as well as real H₂O.⁷⁷ Each of the three fluids exhibits the divergence to extreme compressibility (κ_T) and thermal expansion coefficient (α_P) in the vicinity of the critical point (Fig. 4). Fig. 4 also highlights that the COS/D2 model provides better reproduction of the critical density of real water, compared to the SPC/E model (left panels), and somewhat better reproduction of the near-critical values of both κ_T and α_P (left and right panels). Again, both the SPC/E and COS/D2 models provide good quantitative reproduction of the near-critical derivative properties, which are key parameters for modeling thermodynamics of solvation of electrolytes and non-electrolytes.^4–6 The second virial coefficients of the COS/D2 water model from 300–450 °C show fairly good quantitative agreement with experiment (Fig. 5), providing marked improvement compared to the nonpolarizable models SPC, SPC/E and TIP4P/2005, although the polarizable BKd3 model appears to reproduce the second virial coefficient of real water better than COS/D2 (Fig. 5).

Fig. 4 Isothermal compressibility (κ_T; top) and isobaric thermal expansivity (α_P; bottom) of real, COS/D2 and SPC/E H₂O functions of fluid density along the isotherms of T_c (left) and T_c + 1 °C (right).

Fig. 5 Second virial coefficients of H₂O and water models. The data for SPC and SPC/E water are from literature.⁸⁸ Note that the second virial coefficient of the TIP4P/2005 model overlaps with that of SPC/E,⁵⁵ whereas that of the BKd3 model overlaps with that of real H₂O.⁵⁵

Realistic simulations of electrolyte solutions at hydrothermal conditions require good reproduction of the static permittivity (dielectric constant) of the solvent.¹ In addition, optimization of ion parameters for simulations of concentrated electrolyte solutions requires knowledge of the dielectric permittivity for input into solvation free energy calculations.^65,66 Although various water models have been characterized in detail in terms of dielectric properties at ambient conditions or as a function of temperature at 1 atm pressure, few data are available in literature for dielectric properties of liquids, vapors and supercritical fluids at elevated temperatures and pressures. Fig. 6 shows the density dependence of the dielectric constant of COS/D2 and SPC/E H₂O and comparison to real water.⁷⁴ As expected, both water models converge to unity at the zero-density limit, and show increasing static permittivity with increasing density up to 1 g cm⁻³. In general, the SPC/E provides a somewhat better reproduction of the real-water dielectric constant (Fig. 6). In comparison, the COS/D2 underestimates the dielectric constant by approximately one to three units over the range of intermediate densities. With increasing density towards 1 g cm⁻³ (and temperatures from 300 to 450 °C), the dielectric constant of COS/D2 H₂O seems to approach that of SPC/E H₂O, suggesting that the dielectric constant of COS/D2 may exceed that of SPC/E at ρ > 1 g cm⁻³. Nevertheless, the differences in dielectric constant between COS/D2, SPC/E and real H₂O are relatively small considering the variability in dielectric constant with changing temperature and density, suggesting that the COS/D2 model can be used for simulations of electrolyte solutions at elevated temperature and pressure, without taking flexible water models into account, which reproduce the static permittivity of real water even better.^84,85

Fig. 6 Density dependence of the relative static permittivity (dielectric constant) of real, COS/D2 and SPC/E H₂O at 300, 350, 400 and 450 °C.

Molecular simulations of liquids and solutions at hydrothermal conditions are also useful for assessing interatomic interactions and fluid structure. Therefore, we compared the structural properties of COS/D2 and SPC/E H₂O to test the suitability of the COS/D2 model for structural studies. Fig. 7 shows the radial distribution functions, g(r), of oxygen–oxygen, oxygen–hydrogen and hydrogen–hydrogen at 300 and 400 °C, at bulk densities of 0.2 to 1 g cm⁻³ at 0.2 g cm⁻³ increments. At both temperatures, and from 0.4 to 1 g cm⁻³, both models show very similar radial distribution functions (Fig. 7). At lower density (representing vapor-like conditions), the peaks for SPC/E H₂O are more pronounced than those of COS/D2 (Fig. 7), although this trend is reversed at liquid-like densities and peaks for COS/D2 are somewhat more pronounced than those for SPC/E. This latter trend (more pronounced peaks for COS/D2 versus SPC/E) is similar to observations for H₂O liquid at ambient conditions.⁶⁰ Nevertheless, the good agreement between the two models, especially at ρ > 0.4 g cm⁻³, suggests that COS/D2 shows overall similar fluid structure to SPC/E model H₂O at hydrothermal conditions.

Fig. 7 Comparison of oxygen–oxygen (top), oxygen–hydrogen (middle) and hydrogen–hydrogen (bottom) radial distribution functions of COS/D2 versus SPC/E H₂O at 300 (left) and 400 °C (right). Results based on neutron diffraction data for real H₂O⁸⁶ are plotted for comparison. Note that the fluid densities in the neutron diffraction experiments⁸⁶ do not correspond exactly to the densities presented here (Exp. densities are: 0.72, 0.78, 0.83, 0.89 and 0.92 g cm⁻³ at 300 °C; and 0.58, 0.63, 0.73 and 0.87 g cm⁻³ at 400 °C (ref. 86)), and are superimposed on the nearest density data from our simulations.

In Fig. 7, we also compare the simulated g(r) results with data obtained from neutron diffraction experiments on H₂O at 300 and 400 °C.⁸⁶ It should be noted that structure factors obtained from neutron diffraction experiments require model calculations to extract radial distribution functions, and as such, radial distribution functions extracted from diffraction data are not purely experimental in nature. For example, the data shown in Fig. 7 were obtained from diffraction data using an empirical potential structure refinement technique, which used the SPC/E potential as a starting point reference.⁸⁶ Nevertheless, this comparison illuminates some of the key features of the fluid structure, which are consistently displayed by the COS/D2 (and SPC/E) model. Firstly, note that the first sharp peak in the oxygen–hydrogen (r ∼ 0.09 nm) and hydrogen–hydrogen (r ∼ 0.15 nm) radial distribution functions from the diffraction data (Fig. 7) correspond to intramolecular correlations, which are omitted from the calculations from molecular simulations (Fig. 7). The hydrogen–hydrogen data shows good consistency between neutron-diffraction and simulation results in terms of first and second intermolecular maxima (Fig. 7). The oxygen–hydrogen correlation shows good consistency in terms of the second intermolecular maximum, but the first peak is significantly subdued in the neutron-diffraction results compared to simulations (e.g., at 400 °C, becoming a shoulder rather than a real peak). Soper⁸⁶ reported that the decrease/disappearance of the first intermolecular oxygen–hydrogen peak with increasing temperature and/or decreasing density corresponds to a decrease in hydrogen bonding. Again, COS/D2 (and SPC/E) do not show such a decrease in the first oxygen–hydrogen maximum. However, it should also be noted that other neutron diffraction studies⁸⁷ have reported oxygen–hydrogen radial distribution functions exhibiting an obvious first inter-molecular peak, similar to those presented here. Additional study is needed, to clarify the oxygen–hydrogen correlation in H₂O at these conditions. Finally, as also discussed by Soper,⁸⁶ the first maximum in the simulated oxygen–oxygen radial distribution function is somewhat more pronounced than that obtained from neutron diffraction data. In addition, at 300 °C the second maximum in the oxygen–oxygen g(r) is somewhat broader and located at higher radii in the simulation results, compared to the diffraction data. This latter phenomenon suggests that the oxygen–oxygen correlation is weighted towards short-range interactions in the simulations, such that the longer-range interactions may be underestimated.⁸⁶ Soper⁸⁶ described this as a general drawback of empirical potentials for H₂O, and suggested that incorporating molecular polarization may compensate for the observed discrepancy. It is thus interesting to note that the agreement in between the neutron diffraction data and the simulations seems to improve somewhat at higher temperature (400 °C; Fig. 7).

Comprehensive simulations of concentrated electrolyte solutions (“brines”) using the COS/D2 model will require optimization of ion parameters specific for this solvent; as such, we have not attempted detailed simulations of that nature in the present study. However, as an initial test of ion–solvent interaction in hydrothermal COS/D2, we have simulated Na⁺ solvation at dilute conditions (∼0.01 molal), at 0.1 g cm⁻³ and 375 °C (Fig. 8). It should be noted that the ion parameters used in the present study were optimized for the SPC/E solvent,⁶⁵ and not for COS/D2. Nevertheless, the radial distribution functions and running coordination numbers (Fig. 8) provide a crude comparison to assess ion solvation in COS/D2 at dilute conditions. Firstly, the excellent agreement between solvation structures of Na⁺ in SPC/E versus COS/D2 is notable (Fig. 8), despite the somewhat different Lennard-Jones parameters and partial charges on atoms between these two models.⁶⁰ Secondly, in the comparison of running coordination numbers (density-normalized, integrated g(r)), both models yield approximately the same value at the first minimum in the radial distribution function (∼5), which represents the conventional hydration number. This value of ∼5 for the Na⁺ hydration number consistent with previous simulations of dilute NaCl at similar conditions.⁷ At longer radial distances (beyond the first minimum in g(r)) there is a slight tendency towards higher coordination of SPC/E compared to COS/D2. This latter phenomenon implies that the density perturbation of the solvent around ions⁵ extends further in SPC/E H₂O compared to polarizable COS/D2. This is the expected result, because polarization of H₂O in the vicinity of the ion enhances its solvation.

Fig. 8 Solvation structure of dilute (∼0.01 molal) Na⁺ in COS/D2 and SPC/E, 375 °C and bulk density 0.1 g cm⁻³. Black curves represent the sodium–oxygen radial distribution function (g(r)) while grey curves represent the density-normalized integrated radial distribution function (=running coordination number). Detail of the first minimum and second maximum of g(r) shown in the inset.

IV. Conclusions

The polarizable (charge-on-spring) water model COS/D2 provides good quantitative agreement with the pVT and dielectric properties of H₂O from 300 to 450 °C and over a wide range of fluid densities (from vapor-like to liquid-like), as well as fair reproduction of the liquid–vapor curve and critical point of H₂O. We should emphasize that several other recent (nonpolarizable and polarizable) models exhibit comparable or better reproduction of the critical properties of water (Table 2), although these models have generally not been characterized throughout the full range of conditions studied here (including single-phase liquid, vapor and supercritical fluid conditions). In particular, the good performance of the SPC/E model is once again demonstrated by the present results, as also previously emphasized by Vega and Abascal.⁶³ The good performance of SPC/E is noteworthy especially because this relatively simple model continues to outperform several more recent and “sophisticated” models for a variety of hydrothermal properties. Nevertheless, some of these alternative models appear quite promising. The SPC/E model, like other non-polarizable models, is likely to miss key aspects of structure and speciation in concentrated electrolyte solutions at hydrothermal conditions. The COS/D2 model is suitable for scenarios in which accounting for polarizability is a key objective, including at hydrothermal conditions. We do not suggest that the COS/D2 model is necessarily the “best” model for any particular problem or objective. Advantages of the COS/D2 model (beyond its good reproduction of the physical properties of H₂O) include its compatibility with the GROMOS forcefield, and its relative simplicity – as well as its concomitantly fast computation time even for rather large simulations. The good agreement of the COS/D2 model with experiment is especially noteworthy considering that the model was parameterized only according to properties at ambient temperature and pressure.⁶⁰ Thus, the model performs well at conditions far from those of the parameterization, providing indirect affirmation that the charge-on-spring formulation is well suited for simulating H₂O at such conditions. Furthermore, the results of this study indicate that the COS/D2 water model can be used reliably for simulations of aqueous fluids and solutions at high temperatures and pressures (provided ion parameters optimized for use with this water model). The present study represents a step towards including solvent polarizability in simulations of aqueous solutions at hydrothermal conditions; we expect that the COS/D2 model will be particularly useful for modeling solutions of concentrated and/or highly charged electrolytes, in which dielectric properties and ion solvation are key parameters.

Acknowledgements

The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007–2013) under grant agreement no. 331171, and from the Swiss National Science Foundation (SNF) under SNF grant #200021_132320.

References

1. A. G. Kalinichev, Rev. Mineral. Geochem., 2001, 42, 83–130.
2. Y. Marcus, Supercritical Water: A Green Solvent: Properties and Uses, John Wiley & Sons, Inc., 2012.
3. J. Reimer, M. Steele-MacInnis, J. M. Wambach and F. Vogel, J. Phys. Chem. B, 2015, 119, 9847–9857.
4. A. A. Chialvo, P. T. Cummings, J. M. Simonson and R. E. Mesmer, J. Chem. Phys., 1999, 110, 1075–1086.
5. A. A. Chialvo, P. T. Cummings, J. M. Simonson and R. E. Mesmer, J. Chem. Phys., 1999, 110, 1064–1074.
6. A. A. Chialvo, P. G. Kusalik, P. T. Cummings and J. M. Simonson, J. Chem. Phys., 2001, 114, 3575–3585.
7. T. Driesner, T. M. Seward and I. G. Tironi, Geochim. Cosmochim. Acta, 1998, 62, 3095–3107.
8. K. Yui, M. Sakuma and T. Funazukuri, Fluid Phase Equilib., 2010, 297, 227–235.
9. J. Brodholt and B. Wood, Am. Mineral., 1993, 78, 558–564.
10. C. L. Lin and R. H. Wood, J. Phys. Chem., 1996, 100, 16339–16409.
11. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren and J. Hermans, in Intermolecular Forces, ed. B. Pullman, Riedel, Dordrecht, 1981, pp. 331–342.
12. H. J. C. Berendsen, J. R. Grigera and J. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
13. M. M. Reif, P. H. Hünenberger and C. Oostenbrink, J. Chem. Theory Comput., 2012, 8, 3705–3723.
14. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
15. C. Vega, J. L. F. Abascal and I. Nezbeda, J. Chem. Phys., 2006, 125, 034503.
16. J. Applequist, Acc. Chem. Res., 1977, 10, 79–85.
17. J. Applequist, J. R. Carl and K.-K. Fung, J. Am. Chem. Soc., 1972, 94, 2952–2960.
18. C. J. Burnham, J. Li, S. S. Xantheas and M. Leslie, J. Chem. Phys., 1999, 110, 4566–4581.
19. C. J. Burnham and S. S. Xantheas, J. Chem. Phys., 2002, 116, 1500–1510.
20. L. Jensen, P.-O. Åstrand, A. Osted, J. Kongsted and K. V. Mikkelsen, J. Chem. Phys., 2002, 116, 4001–4010.
21. L. Silberstein, Philos. Mag., 1917, 33, 92–128.
22. L. Silberstein, Philos. Mag., 1917, 33, 215–222.
23. L. Silberstein, Philos. Mag., 1917, 33, 521–533.
24. B. T. Thole, Chem. Phys., 1981, 59, 341–350.
25. P. T. van Duijnen, A. H. de Vries, M. Swart and F. Grozema, J. Chem. Phys., 2002, 117, 8442–8453.
26. P. T. van Duijnen and M. Swart, J. Phys. Chem. A, 1998, 102, 2399–2407.
27. U. Niesar, G. Corongiu, E. Clementi, G. R. Kneller and D. K. Bhattacharya, J. Chem. Phys., 1990, 94, 7949–7956.
28. R. Kozack and P. Jordan, J. Chem. Phys., 1992, 96, 3120.
29. G. Corongiu, Int. J. Quantum Chem., 1992, 42, 1209.
30. J. Brodholt, M. Sampoli and R. Vallauri, Mol. Phys., 1995, 86, 149.
31. A. Chialvo and P. Cummings, Fluid Phase Equilib., 1998, 150, 73.
32. A. Chialvo and P. Cummings, J. Chem. Phys., 1996, 105, 8274.
33. L. Dang and T. Chang, J. Chem. Phys., 1997, 106, 8149.
34. K. Kiyohara, K. E. Gubbins and Z. A. Panagiotopoulos, Mol. Phys., 1998, 94, 803.
35. P. Cieplak, P. Kollman and T. Lybrand, J. Chem. Phys., 1990, 92, 6755.
36. C. J. Burnham, J. Li, S. S. Xantheas and M. Leslie, J. Chem. Phys., 1999, 110, 4566.
37. C. J. Burnham and S. S. Xantheas, J. Chem. Phys., 2002, 116, 1500.
38. T. Ishiyama and A. Morita, J. Phys. Chem. C, 2007, 111, 721.
39. P. Ren and J. W. Ponder, J. Phys. Chem. B, 2003, 107, 5933.
40. P. Tröster, K. Lorenzen, M. Schwörer and P. Tavan, J. Phys. Chem. B, 2013, 117, 9486–9500.
41. M. Sprik and M. Klein, J. Chem. Phys., 1988, 89, 7556.
42. S. Rick, S. Stuart and B. Berne, J. Chem. Phys., 1994, 101, 6141.
43. B. Chen, J. Xing and J. Siepmann, J. Phys. Chem. B, 2000, 104, 2391.
44. L. P. Wang, J. Chen and T. van Voorhis, J. Chem. Theory Comput., 2012, 9, 452.
45. I. Svishchev, P. Kusalik, J. Wang and R. Boyd, J. Chem. Phys., 1996, 105, 4742.
46. P. J. van Maaren and D. van der Spoel, J. Phys. Chem. B, 2001, 105, 2618.
47. H. Yu, T. Hansson and W. van Gunsteren, J. Chem. Phys., 2003, 118, 221.
48. G. Lamoureux, E. Harder, I. V. Vorobyov, B. Roux and A. D. MacKerrell, Chem. Phys. Lett., 2006, 418, 245–249.
49. G. Lamoureux, A. D. MacKerrell and B. Roux, J. Chem. Phys., 2003, 119, 5185–5197.
50. P. Paricaud, M. Předota, A. Chialvo and P. Cummings, J. Chem. Phys., 2005, 122, 244511.
51. A.-P. E. Kunz and W. F. van Gunsteren, J. Phys. Chem. A, 2009, 113, 11570–11579.
52. A. Baranyai and P. Kiss, J. Chem. Phys., 2010, 133, 144109.
53. L. Viererblová and J. Kolafa, Phys. Chem. Chem. Phys., 2011, 13, 19925.
54. P. T. Kiss and A. Baranyai, J. Chem. Phys., 2012, 137, 084506.
55. P. T. Kiss and A. Baranyai, J. Chem. Phys., 2012, 137, 194103.
56. P. Paricaud, M. Predota, A. A. Chialvo and P. T. Cummings, J. Chem. Phys., 2005, 122, 244511.
57. P. T. Kiss and A. Baranyai, J. Chem. Phys., 2013, 138, 204507.
58. P. T. Kiss and A. Baranyai, J. Chem. Phys., 2014, 140, 154505.
59. P. Lopes, B. Roux and A. MacKerell, Theor. Chem. Acc., 2009, 124, 11.
60. S. Bachmann and W. F. van Gunsteren, J. Chem. Phys., 2014, 141, 22D515.
61. E. Wasserman, B. Wood and J. Brodholt, Geochim. Cosmochim. Acta, 1995, 59, 1–6.
62. Y. Guissani and B. Guillot, J. Chem. Phys., 1993, 98, 8221–8235.
63. C. Vega and J. L. F. Abascal, Phys. Chem. Chem. Phys., 2011, 13, 19663–19688.
64. E. Wernersson and P. Jungwirth, J. Chem. Theory Comput., 2010, 6, 3233–3240.
65. M. M. Reif and P. H. Hünenberger, J. Chem. Phys., 2011, 134, 144104.
66. M. M. Reif and P. H. Hünenberger, J. Chem. Phys., 2011, 134, 144103.
67. A.-P. E. Kunz, J. R. Allison, D. P. Geerke, B. A. C. Horta, P. H. Hünenberger, S. Riniker, N. Schmid and W. F. van Gunsteren, J. Comput. Chem., 2012, 33, 340–353.
68. N. Schmid, C. D. Christ, M. Christen, A. P. Eichenberger and W. F. van Gunsteren, Comput. Phys. Commun., 2012, 183, 890–903.
69. W. F. van Gunsteren, The GROMOS software for (bio)molecular simulation.
70. J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. Phys., 1977, 23, 327–341.
71. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
72. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. di Noa and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
73. J. A. Barker and R. O. Watts, Mol. Phys., 1973, 26, 789–792.
74. D. P. Fernández, A. R. H. Goodwin, E. W. Lemmon, J. M. H. Levelt Sengers and R. C. Williams, J. Phys. Chem. Ref. Data, 1997, 26, 1125–1166.
75. M. Neumann, Mol. Phys., 1983, 50, 841–858.
76. J. C. Maxwell, Nature, 1875, 11, 357–359.
77. W. Wagner and A. Pruß, J. Phys. Chem. Ref. Data, 2002, 31, 387–535.
78. F. Römer, A. Lervik and F. Bresme, J. Chem. Phys., 2012, 137, 074503.
79. C. Pinilla, A. H. Irani, N. Seriani and S. Scandolo, J. Chem. Phys., 2012, 136, 114511.
80. T. M. Yigzawe and R. J. Sadus, J. Chem. Phys., 2013, 138, 044503.
81. S. M. Sterner and R. J. Bodnar, Geochim. Cosmochim. Acta, 1984, 48, 2659–2668.
82. W. A. Bassett, A. H. Shen, M. Bucknum and I. M. Chou, Rev. Sci. Instrum., 1993, 64, 2340–2345.
83. J. J. de Pablo, J. M. Prausnitz, H. J. Strauch and P. T. Cummings, J. Chem. Phys., 1990, 93, 7355–7359.
84. G. Raabe and R. J. Sadus, J. Chem. Phys., 2011, 134, 234501.
85. Y. Wu, H. L. Tepper and G. A. Voth, J. Chem. Phys., 2006, 124, 024503.
86. A. K. Soper, Chem. Phys., 2000, 258, 121–137.
87. T. Tassaing, M.-C. Bellissent-Funel, B. Guillot and Y. Guissani, Europhys. Lett., 1998, 42, 265–270.
88. G. C. Boulougouris, G. E. Ioannis and D. N. Theodorou, J. Phys. Chem. B, 1998, 102, 1029–1035.
89. A. Z. Panagiotopoulos, Mol. Phys., 1987, 61, 813.
90. I. M. Svishchev and T. M. Hayward, J. Chem. Phys., 1999, 111, 9034.
91. E. M. Yezdimer and P. T. Cummings, Mol. Phys., 1999, 97, 993.
92. J. L. F. Abascal and C. Vega, J. Chem. Phys., 2005, 123, 234505.
93. R. Sakamaki, A. K. Sum, T. Narumi and K. Yasuoka, J. Chem. Phys., 2011, 134, 124708.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra13495a

---