Stress testing the ELBA water model

Abstract

The ELBA coarse-grained model describes a water molecule as a single-site Lennard-Jones particle embedded with a point dipole. ELBA was previously reported to capture several properties of real water with relatively high accuracy, while being up to two orders of magnitude more computationally efficient than atomistic models. Here, we ‘stress test’ the ELBA model by investigating the temperature and pressure dependences of two most important water properties, the liquid density and the self-diffusion coefficient. In particular, molecular dynamics simulations are performed spanning temperatures from 268 K up to 378 K and pressures from 1 atm up to 4000 atm. Comparisons are made with literature data from experiments and from simulations of traditional three-site atomistic models. Remarkably, the ELBA results show an overall similar (and sometimes higher) accuracy with respect to the atomistic data. We also calculate a number of additional thermodynamic properties at ambient conditions, namely isothermal compressibility, shear viscosity, isobaric heat capacity, thermal expansion coefficient and melting point. The accuracy of ELBA is relatively good compared to atomistic and other coarse-grained models.

Keywords::

• water
• molecular dynamics
• coarse-graining
• temperature/pressure dependence

1. Introduction

The modelling and simulation of water is an area of great interest for both academic research and industrial applications. Numerous water models, with various characteristics and capabilities, have been developed over several decades now.1234 A popular approach involves atomistic models where each atom in the water molecule is represented by a simulated site,567 although optimised models also exist which include one or two extra sites.8910 While the simulations of atomistic models can provide accurate and realistic results, they are also highly demanding of computational resources. To reduce the computational cost, simplified models can be developed by representing one or several water molecules with lower resolution single sites; in this approach, known as ‘coarse-graining’, some atomic-level details are sacrificed to increase simulation efficiency.1112131415

Recently, the ELBA coarse-grained model has been developed to attempt striking a new balance between physical accuracy and simulation efficiency.161718 In the ELBA model, the three atoms of a water molecule are reduced to a single particle, with an electrical point dipole affixed to its centre to capture the well-known dipolar nature of water. Previous molecular dynamics simulations have shown relatively good accuracy of the ELBA model under ambient conditions of temperature and pressure. In particular, ELBA reproduced several fundamental water properties as accurately as the best atomistic models, including density, potential energy, self-diffusion coefficient, heat of vaporisation, critical point, surface tension and the liquid–vapour interface structure.[18] A significant issue was however observed for the radial distribution function, which showed qualitative and quantitative differences compared with experimental and atomistic data.[18] Notably, ELBA can also be used to hydrate atomistic solutes, including small organic molecules and proteins; in particular, uniquely for a coarse-grained model, ELBA was shown to be directly compatible with common atomistic force fields, meaning that no additional or ad hoc scaling factors, intermediate regions or extra sites were required.[19] In terms of computational cost, ELBA proves to be between one to two orders of magnitude more efficient to simulate than traditional multi-site atomistic models.[18]

While ambient conditions are obviously of primary importance, there is also substantial interest in applying molecular simulations to investigate phenomena which take place under a variety of temperatures and pressures, such as protein folding processes [20,21] and water filtration.[22] Hence, it is important to assess how well a water model can predict real water behaviour at non-ambient conditions. In this work, the performance of the ELBA model (with no modification to the original parameters) is examined over a wide range of temperatures (268–378 K) and pressures (1–4000 atm). We investigate the temperature and pressure dependences of two fundamental water properties, the liquid phase density and self-diffusion, which ELBA reproduces accurately under ambient conditions.[18] Comparisons are made with corresponding experimental data and simulation results for the atomistic models SPC,[5] SPC/E [7] and TIP3P,[6] which represent the most widely used water models in molecular simulation. As extra findings, we also compute the pressure dependence of liquid density and self-diffusion coefficient at room temperature for, respectively, TIP3P and SPC. To our knowledge, these properties have not been previously published.

In addition, five thermodynamic properties at ambient conditions are reported here for ELBA for the first time, namely, isobaric heat capacity, thermal expansion coefficient, isothermal compressibility, melting point and shear viscosity. For these properties, literature data exist also from simulation of other coarse-grained models, which are thus included in the comparison.

2. Methods

2.1 ELBA coarse-grained model

The ELBA model describes a water molecule as a single Lennard-Jones particle embedded with a point dipole (Figure 1).[16,18] The potential energy for a pair of water particles i and j is the sum of a Lennard-Jones interaction term and a dipole interaction term :

(1)

Figure 1 (Colour online) Atomistic water molecule and the ELBA model. The left panel shows a water molecule at the atomic level, with a negative charge (‘ − ’ sign) on the oxygen atom and two positive charges (‘+’ signs) on the hydrogen atoms. The right panel depicts an ELBA water site, with the arrow representing an electrical point dipole.

Both terms are in the ‘shifted-force’ form.[23,24] whereby the potential energy and its derivative (the force) go to zero smoothly at the cut-off point. The shifted-force scheme removes cut-off-related artefacts and simulation stability problems,[23,24] which are especially severe for orientation-dependent potentials such as the point dipole potential.[25] For the Lennard-Jones part, the following expression proposed by Stoddard and Ford [26] is adopted:

(2)

in which σ and ϵ have the standard meaning,[23,24] r is the inter-particle distance and r_c is the cut-off radius. For the dipole interactions, the classic electrostatic model [23,27] has been modified [17,18] as:

(3)

where ϵ₀ is the vacuum permittivity, μ_i and μ_j are the point dipole vectors of the interacting pair, r and r are, respectively, the distance vector and its magnitude and r_c is the cut-off distance. Values of the potential parameters used here are the same as in previous works [18,19]: ϵ = 0.55 kcal mol ^− 1, , μ = 2.6 D and .

2.2 General simulation details

Molecular dynamics simulations were performed using the software LAMMPS (version 11 Nov 2013).[28,29] Input scripts are available on our group website.[30] The simulation systems contained 8000 water sites in a cubic region. Conventional periodic boundary conditions were adopted. Conditions of constant temperature and pressure were applied in most of the simulations; in these cases, the edge length of the starting simulation box was set to 6.2 nm, yielding an initial water density of approximately 1 g cm ^− 3 (consistent with the density of real water at standard ambient conditions). However, in some cases (detailed below), we fixed the volume, and hence also the system density, to a constant value. In all simulations, the temperature was controlled using the Langevin thermostat,[31] with a collision frequency of 1 ps ^− 1. The pressure was maintained (when needed) using the barostat by Berendsen et al., [32] with a damping time of 1 ps and a bulk modulus of 2.174 × 10⁴ atm. The integration time step for the ELBA simulations was 10 fs. At every time step, the net momentum of the mass centre of the entire system was removed to prevent drifting during the simulation. The interaction cut-off radius was 12 Å and no long range interactions (beyond the cut-off distance) were included for the ELBA model. All these settings are consistent with previous work. [18] Additional atomistic simulations were also run to study the pressure dependence of the TIP3P density and the SPC diffusion coefficient, as no corresponding results were found in the literature. For these simulations the time step was 2 fs. Bonds and angles were constrained using the SHAKE algorithm [33] with a relative tolerance of 10 ^− 4. The non-bonded cut-off distance was 10 Å for SPC [18] and 13 Å for TIP3P.[18,34] Long range electrostatic interactions were included with the particle-particle particle-mesh (PPPM) method [35] with a relative tolerance of 10 ^− 5.

Specific details of individual calculations are given in Section 2.3.

2.3 Details of individual calculations

To investigate density, self-diffusion, isothermal heat capacity, thermal expansion coefficient and isothermal compressibility, we ran series of 7 ns long simulations. The initial 2 ns was regarded as equilibration, while the subsequent 5 ns was considered as production; during production, the relevant properties were sampled every 0.1 ps. For each simulation, three independent repeats were run by assigning the initial velocities with different random seeds.

Calculations of density and self-diffusion coefficient followed standard procedures. [2,24] Regarding isothermal heat capacity, thermal expansion coefficient and isothermal compressibility, details are given in the following Sections 2.3.1, 2.3.2 and 2.3.3, respectively. The calculations of shear viscosity and melting point were carried out according to the protocols described below in Sections 2.3.4 and 2.3.5, respectively.

2.3.1 Isobaric heat capacity

The isobaric heat capacity C_p at 298 K was estimated as [36]:

(4)

where and are the total energies per molecule at temperature T₁ = 288 K and T₂ = 308 K. The second term at the right-hand side of the equation is the quantum contribution, which for classical molecular simulations is about  − 9.3 J mol ^− 1K ^− 1 under standard ambient conditions.363738

2.3.2 Thermal expansion coefficient

The thermal expansion coefficient α at 298 K was estimated using the finite-difference expression method [37,39]:

(5)

where V is the simulation box volume, and and are the densities at temperature T₁ (288 K) and T₂ (308 K), respectively.

2.3.3 Isothermal compressibility

The isothermal compressibility was calculated from two sets of constant-volume (and hence also constant-density) simulations, by applying the finite-difference equation as follows [36,37]:

(6)

where and are the pressures corresponding to the constant density values of, respectively, and . These density values were fixed at the outset by setting the system volume accordingly.

2.3.4 Shear viscosity

The shear viscosity was calculated according to the Green–Kubo formula,[40,41] in which the autocorrelation function of the stresses is integrated as where V is the system volume, is the Boltzmann constant, T is the system temperature and are the off-diagonal components of the stress tensor . The calculation was carried out over 1 ns production stages following 1 ns equilibration stages.[42] A total of six independent simulations were run by setting different initial velocities.

2.3.5 Melting point

The liquid–solid coexistence method [43,44] was followed in order to determine the melting point. The simulation system consisted of a cubic region containing 8000 water sites in total (as described above, in Section 2.2). The cube was divided into two equal regions, designated liquid and ice region, respectively. Then three separate simulation stages were performed. In the first stage, lasting 1 ns, the whole system was cooled down to 200 K to achieve complete freezing of both regions. In the second stage, lasting 1 ns, particles assigned to the ice region were immobilised (thus retaining their frozen state), while the temperature for the liquid region was raised to 350 K to ensure complete melting. Finally, in the third stage, lasting 10 ns, the whole system was equilibrated to a series of target temperatures, with all particles freely moving. During freezing, the thermostat subtracts energy from the system, therefore the total energy decreases. Conversely, during melting, the thermostat adds energy to the system, thus the total energy increases.[44] By examining the total energy evolution over time and the final state of these trial simulations, the melting point was gradually approached. In the simulations where the total energy neither increased nor decreased and both ice and liquid phases coexisted at a specific temperature, this temperature was considered to be the melting point. Alternatively, if no coexistence appeared, but the temperature difference between a completely melted and a completely frozen system was equal to 1 K, the melting point was estimated as the average between them (see Figure 2, and Section 1 in the Supplementary Material). In total, the three-stage procedure described was repeated five times with different initial conditions.

Figure 2 (Colour online) Schematic diagram of the melting point calculation procedure. See main text for a detailed description of the different stages.

3. Results and discussion

3.1 Temperature dependence

The density and the self-diffusion coefficient of the ELBA water have been calculated at temperatures ranging from 268 and 378 K, under a constant pressure of 1 atm. Results for both properties are compared with previously published results from experiment and from simulations of the widely used three-site atomistic water models SPC, SPC/E and TIP3P.645464748

Figure 3 shows the density data as a function of temperature. Compared with the experimental data, it is clear that the ELBA results exhibit a larger rate of change with temperature in comparison to real water. Specifically, the ELBA results underestimate the density above room temperature while overestimate it below room temperature. Note, however, that the discrepancies observed are not particularly severe; the largest error, at around 373 K, is ∼0.03 g cm ^− 3, which corresponds to a relative error of only ∼3% with respect to the experimental measurement. When compared to the data from the atomistic models, the ELBA results are remarkably accurate (where accurate is intended to mean close to the experimental value). In particular, ELBA proves more accurate than SPC for temperatures above ∼285 K. ELBA is also more accurate than TIP3P for temperatures above ∼290 K. However, ELBA is less accurate than SPC/E for all temperatures.

Figure 3 (Colour online) Density as a function of temperature at 1 atm. The standard deviation of the ELBA data is less than 0.00005 g cm^− 3 for all values; corresponding error bars are smaller than the size of the symbols. The experimental data are from the NIST database.[73] Atomistic data are from Paschek. [74]

The self-diffusion coefficient data as a function of temperature are displayed in Figure 4. It can be seen that the ELBA results exhibit a lower rate of change with temperature in comparison to real water. In terms of absolute values, the accuracy of ELBA is rather good in a fairly large region centred around room temperature (approximately from 270 to 320 K). In this region, ELBA is substantially more accurate than SPC and TIP3P, and largely as accurate as SPC/E. For temperatures above ∼320 K, the ELBA values increasingly diverge from the experimental curve, indicating that the mobility of the ELBA water is less sensitive than real water to the temperature increase.

Figure 4 (Colour online) Self-diffusion coefficient as a function of temperature at 1 atm. The standard deviation for the ELBA results is at most 0.05 × 10^− 9 m² s^− 1; corresponding error bars are of the order of the size of the symbols. The experimental data,[75] as well as results for SPC,[48] SPC/E [48] and TIP3P [47] are shown for comparison.

3.2 Pressure dependence

The pressure dependence of ELBA was investigated by applying external isotropic pressures from 1 to 4000 atm while maintaining a standard ambient temperature of 298 K. From Figure 5, it can be seen that the ELBA results are close to experiment only in the initial region around ambient condition, while they increasingly underestimate the experimental values with increase in pressure. This implies that ELBA is less compressible by high pressure than real water. Compared with the atomistic models, ELBA is markedly less accurate than both SPC/E and TIP3P, but it proves more accurate than SPC for pressures below 1000 atm.

Figure 5 (Colour online) Density as a function of pressure at 298 K of ELBA, together with experimental data [73] and simulation results for SPC,[76] SPC/E [76] and TIP3P (this work). The standard deviations for the ELBA and TIP3P results are at most 0.00005 and 0.00007 g cm^− 3, respectively; corresponding error bars are smaller than the size of the symbols.

Regarding the self-diffusion coefficient, Figure 6 shows that the accuracy of ELBA is rather high at ambient pressure, but it decreases with increase in pressure, as the experimental values are more and more underestimated. Nevertheless, ELBA proves markedly more accurate than both SPC and (especially) TIP3P, while the SPC/E data almost match the experimental values.

Figure 6 (Colour online) Self-diffusion coefficient as a function of pressure at 298 K. Results from experiment [77] and for SPC (this work), SPC/E [78] and TIP3P [79] are shown for comparison. The standard deviation for the ELBA and SPC values are at most 0.018 × 10^− 9 and 0.23 × 10^− 9 m² s^− 1, respectively; corresponding error bars are of the order of the size of the symbols.

3.3 Thermodynamic properties at ambient conditions

Results for five thermodynamic properties of the ELBA model at 298 K and 1 atm are presented in Table 1, together with literature data from experiments, three-site atomistic models and also other coarse-grained models. The atomistic models are SPC,[5] SPC/E [7] and TIP3P,[6] already used as standards for comparison in the temperature and pressure dependences reported in the previous section. Regarding the other coarse-grained models, comparable results have been reported for single-site models (like ELBA) and for ‘coarser’ descriptions, whereby several water molecules are reduced to a smaller number of coarse-grained particles. Considering single-site models, relevant data have been reported in the literature for the soft sticky dipole (SSD), the model of water (mW) and multiscale coarse-grained (MS-CG) models. SSD is a model comprising a Lennard-Jones sphere, a point dipole and an octupolar term used to capture hydrogen bonding. SSD has been used to study water in the liquid and ice phases,49505152 and as a solvent for different systems.5253545556 The mW model is characterised by a three-body term that captures tetrahedral coordination.[57] MS-CG is a model consisting of an effective potential obtained by the force-matching method.[58] The big multipole water (BMW) model [59] is a coarser model that reduces four water molecules to three charged sites, with the Born–Mayer–Huggins potential replacing the conventional Lennard-Jones potential. BMW has been applied to study the solvation and aggregation of peptides.[60] MARTINI is a popular model based on the Lennard-Jones potential, without explicit electrostatics, where each site represents four water molecules.[61] P-MARTINI [62] refers to a polarisable variant of MARTINI which comprises three sites, two of which are oppositely charged. The GROMOS model [36] maps five water molecules to two oppositely charged sites. This model has also been used as a solvent for proteins.[63] WT4 is a model that represents clusters of 11 water molecules by tetrahedral particles comprising four charged sites.[64] WT4 has been used as a solvent for ions and nucleic acids.[64]

Table 1 Thermodynamic properties at ambient conditions.^a

	Mapping	Cp (J mol^− 1 K^− 1)	α (10^− 4 K^− 1)	κT (10^− 5 atm^− 1)	Tm (K)	η (mPa s)
ELBA	1 → 1	64.37 ± 0.17	9.64 ± 0.03	2.914 ± 0.004	301.7 ± 9.8	0.952 ± 0.093
Experiment^b		75.3	2.56	4.53	273.15	0.8931
SPC^c	1 → 3	75.6	7.3, 10.6^d	4.74	190^e	0.49
SPC/E^b	1 → 3	86.6	5.6^c	4.67, 5.17^c	215^e	0.729
TIP3P^b	1 → 3	78.41	9.2^e, 14.4^d	5.81	146^e	0.321
SSD^f	1 → 1	120.5 ± 0.5	–	–	Above 375	–
mW^g	1 → 1	33	–	1.9	274.6	–
MS-CG^h	1 → 1	–	25.0	14.7	–	–
BMW^i	4 → 3	–	–	3.3	Below 200	1.01 to 1.62^j
P-MARTINI^k	4 → 3	–	–	–	280–285	–
MARTINI^1	4 → 1	–	–	2.6,^i 6.1	290 ± 5	–
GROMOS^m	5 → 2	80.7	23	8.4–13.8	–	3.72
WT4^n	11 → 4	43^o	11.6	2.46	–	–

Notes: Isobaric heat capacity: C_p. Thermal expansion coefficient: α. Isothermal compressibility: κ_T. Melting point: T_m. Shear viscosity: η. The ELBA results are reported as ‘average ±  standard deviation’. A ‘mapping’ entry m → n indicates m water molecules are modelled by n sites.

a Temperature of 298 K and pressure of 1 atm, unless otherwise stated.

b Ref. [1].

c Ref. [37].

d Ref. [45].

e Ref. [47].

f Ref. [50].

g Ref. [57], κ_T at 300 K.

h Ref. [58], at 300 K.

i Ref. [59], at 300 K.

j Ref. [70], at 300 K.

k Ref. [62].

l Ref. [71], at 300 K.

m Ref. [36].

n Ref. [64].

o Ref. [72], value for a WT4 and SPC mixture.

Regarding the isobaric heat capacity C_p, the ELBA result is lower than the experimental value. This could be expected since there is no explicit hydrogen-bonding network in the ELBA model. In fact, the high heat capacity of real liquid water is believed to be partly due to the presence of hydrogen bonds, because when water is heated, the corresponding energy is partly used to bend or break the hydrogen bonds, rather than being totally transferred to the kinetic energy of the water molecules.[65,66] Quantitatively, the difference between ELBA and the experimental value is ∼10 , corresponding to a relative error of 14.5%. Compared with the atomistic models, ELBA proves slightly more accurate than SPC/E, for which the relative error is 15%; however, SPC and TIP3P are in turn markedly more accurate than ELBA, with relative errors of 0.3% and 4%, respectively. As for the other coarse-grained models, it can be seen that mW is characterised by a C_p value that is less than half the experimental value. The WT4 model also underestimates , being 43% lower than the experimental measurement. Both the SSD and GROMOS models instead overestimate , although with rather different relative errors of, respectively, 60% and 7%.

Regarding the thermal expansion coefficient α, it can be seen that the ELBA value is not very accurate, as it is almost four times larger than the experimental value. However, in comparison with the atomistic models, the accuracy of ELBA is relatively good, being within the range reported for SPC and TIP3P (which have α values about three to five times larger than the experimental value). The SPC/E result is the most accurate, even though it is still over twice as large as the experimental result. It should be pointed out that these results are consistent with the temperature-dependent density data reported in Figure 3, where all densities from simulations show a larger rate of change with temperature with respect to the experimental data. Specifically, as the density of the models decreases faster than that experimentally observed with increasing temperature, the thermal expansion coefficient is expected to be overestimated. Compared with the other coarse-grained models, ELBA proves slightly more accurate than WT4 and substantially more accurate than GROMOS and MS-CG.

As for the isothermal compressibility , the ELBA model underestimates the experimental value by ∼36%, while the atomistic models all prove more accurate, with relative errors in the range from 4% to 28%. Compared with the other coarse-grained models, ELBA is instead more accurate than all but one model (BMW).

Regarding the melting temperature, we obtained for ELBA a value of 301.7 ± 4.4 K, which is 28 K higher than the experimental value. By considering only the magnitude of the error with respect to the experimental temperature, ELBA would seem to be more accurate than the atomistic models, which all underestimate the true melting temperature with larger absolute errors. However, it must be stressed that underestimating the melting temperature is a much less serious problem than overestimating it. In most applications, simulations are carried out at (or around) room temperature, where water must be in the liquid phase; in this context, models overestimating the melting temperature may instead freeze. In particular, considering those coarse-grained models overestimating the melting temperature (ELBA, MARTINI, P-MARTINI and SSD), the solid phase is predicted to be the most thermodynamically stable phase for temperatures close to room temperature (or even substantially higher in the case of SSD); as noted, this is in principle very problematic. In practice however, since freezing is an activated process,[67] unless a model is already in the solid phase, simulated liquid water can remain in the (supercooled) liquid phase for temperatures well below the predicted melting temperature. Regarding ELBA, we have run microsecond long simulations at decreasing temperatures, and we have observed spontaneous freezing only for temperatures below 250 K (see Section 2 in the Supplementary Material). Moreover, for any other system we have previously simulated, including the water–vacuum interface,[18] hydrated lipid bilayers in the lamellar and inverse hexagonal phases,[16,17] and hydrated small molecules and proteins,[19] the ELBA water has consistently remained in the liquid phase. To lower the melting temperature, ELBA could be reparametrised. In particular, decreasing the dipole magnitude should have the effect of decreasing the melting temperature, consistently with what has been reported for generic ‘Lennard-Jones plus point dipole’ models in reduced units.[68,69] However, changing the dipole will of course have repercussions not only on the melting point, but also on all other physical properties; a reparametrisation of the whole force field could therefore be needed, with no guarantee of eventually obtaining an overall better model.

Finally, regarding the shear viscosity , ELBA yielded a result only 7% larger than the experimental value; the corresponding accuracy is higher than that of the atomistic models (which have relative errors between 18% and 64%) as well as the other coarse-grained models BMW (relative error from 13% to 81%) and GROMOS (relative error of 316%). Such a high accuracy of ELBA in predicting the shear viscosity was not unexpected, as the shear viscosity is related to the self-diffusion coefficient (through the Stokes–Einstein expression), and ELBA was shown previously to reproduce closely the experimental self-diffusion coefficient.[18]

4. Conclusion

We have studied the temperature and pressure dependences of density and self-diffusion for the coarse-grained ELBA water model. Notably, the accuracy of ELBA was found to be overall comparable to that of the standard three-site atomistic models, in that the relative errors with respect to experiment obtained with ELBA were mostly within the range of errors characterising the atomistic models. In fact, in some cases ELBA was shown to be even more accurate than some of the atomistic models.

Moreover, we computed five properties at standard ambient conditions, namely isothermal compressibility, shear viscosity, isobaric heat capacity, thermal expansion coefficient and melting point, for which comparisons could also include other coarse-grained models. For the thermal expansion coefficient and the isothermal compressibility, ELBA proved less accurate than the atomistic models, but more accurate than the other coarse-grained models. Regarding the heat capacity, ELBA was shown to be as accurate as SPC/E and more accurate than most available coarse-grained models. For the shear viscosity, the ELBA prediction was even closer to experiment than the results reported for the atomistic models.

Together with the recent demonstration of the direct compatibility of ELBA with atomistic force fields for organic molecules and proteins,[19] the work presented here opens up a range of opportunities for mixed atomistic/coarse-grained simulations under non-ambient conditions of temperature and pressure.

Supplementary data

Supplemental data for this article can be accessed doi:10.1080/08927022.2015.104736710.1080/08927022.2015.1047367.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

WD would like to thank the China Scholarship Council for financial support. High-performance computational facilities were provided by the MidPlus regional centre of excellence [grant number EP/K000128/1].

Notes

^1. Email: [email protected]

^2. Email: [email protected]