Liquid-vapor coexistence properties of Water, modeled by the SPC/E Force Field [1], obtained by grand-canonical Wang-Landau/Transition-matrix Monte Carlo and histogram re-weighting. The implementation of the SPC/E model is described in the SPC/E reference configuration calculation described elsewhere on this website, except for three variations described below. Mean values of the saturation pressure, density, and activity (chemical potential- see below) for each phase are reported.
METHOD | Grand-canonical hybrid Wang-Landau/Transition-matrix Monte Carlo and histogram re-weighting [2-4, 10-14] |
Fluid | Water |
Model | SPC/E [1] |
V | 8000 Å^{3} |
TRUNCATION | |
Lennard-Jones | 10 Å + analytic Long-range Corrections |
Electrostatics | 10 Å + Ewald Summation |
Prob. of Disp. Move | 0.4 |
Prob. of Rot. Move | 0.4 |
Prob. of Ins/Del Move | 0.2 |
Biasing Function Update Frequency | 1.0E6 trial moves |
Simulation Length | > 20 sweeps (see below) |
T (K) |
N_{max} |
N_{trials} |
ρ_{vap} (kg/m^{3}) |
+/- |
ρ_{liq} (kg/m^{3}) |
+/- |
p_{sat} (bar) |
+/- |
lnz_{sat} |
+/- |
300 | 296 | 1.458E+11 | 7.373E-03 | 3.253E-05 | 9.981E+02 | 2.928E+00 | 1.017E-02 | 4.516E-05 | -1.524E+01 | 5.724E-02 |
325 | 305 | 1.326E+11 | 3.061E-02 | 8.897E-05 | 9.830E+02 | 1.082E+00 | 4.551E-02 | 1.320E-04 | -1.384E+01 | 2.240E-03 |
350 | 300 | 9.577E+10 | 1.007E-01 | 1.508E-04 | 9.648E+02 | 1.281E+00 | 1.599E-01 | 2.177E-04 | -1.267E+01 | 1.875E-03 |
375 | 295 | 4.792E+10 | 2.761E-01 | 1.718E-04 | 9.438E+02 | 1.177E+00 | 4.635E-01 | 2.867E-04 | -1.168E+01 | 2.879E-03 |
400 | 290 | 2.859E+10 | 6.588E-01 | 2.240E-03 | 9.197E+02 | 4.015E-01 | 1.157E+00 | 3.548E-03 | -1.085E+01 | 1.275E-03 |
425 | 285 | 1.726E+10 | 1.409E+00 | 5.614E-03 | 8.932E+02 | 1.919E+00 | 2.556E+00 | 8.567E-03 | -1.013E+01 | 2.166E-03 |
450 | 280 | 1.143E+10 | 2.778E+00 | 2.839E-02 | 8.632E+02 | 2.276E+00 | 5.136E+00 | 4.432E-02 | -9.525E+00 | 1.885E-03 |
475 | 275 | 8.446E+09 | 5.138E+00 | 2.065E-02 | 8.307E+02 | 1.346E+00 | 9.524E+00 | 4.108E-02 | -8.998E+00 | 2.500E-03 |
500 | 270 | 6.125E+09 | 9.025E+00 | 9.615E-02 | 7.930E+02 | 2.395E+00 | 1.651E+01 | 3.900E-02 | -8.543E+00 | 5.692E-04 |
525 | 265 | 4.406E+09 | 1.547E+01 | 1.812E-02 | 7.494E+02 | 2.921E+00 | 2.719E+01 | 9.094E-02 | -8.148E+00 | 1.752E-03 |
Remarks:
Uncertainties were obtained from three independent simulations and represent 95% confidence limits based on a standard t statistic. Liquid-vapor coexistence was determined by adjusting the activity such that the pressures of the liquid and vapor phases were equal. Here, the pressure is not the conventional virial pressure [5,6] but is the actual thermodynamic pressure, based on the fact that the absolute free energies can be obtained from the distributions determined from simulation [7]. Alternative methods, for example Gibbs-ensemble Monte Carlo and combination grand-canonical Monte Carlo and histogram re-weighting, can be used to determine liquid-vapor coexistence. A review of standard methods of phase equilibria simulations can be found in Ref. 8.
As introduced in Refs. 5 and 6, the activity, z, is defined as
$$ \Large z = \Lambda^{-3} \exp \left( \beta \mu \right) $$
where Λ is the de Broglie wavelength, β = 1/(k_{B}T) (where k_{B} is Boltzmann's constant), and μ is the chemical potential. It is sometimes more convenient to work with ln z in the simulations and in post-processing. The reported activity has units of σ^{-3}. The density is obtained from the simulation in units of average number of molecules per volume. To convert to units of kilograms per cubic meter, the molar mass of a water molecule was taken to be 18.015 28 grams per mole [9].
The Ewald "k" vectors as described in the reference calculation were truncated at k^{2} < 38. Two optimizations to the computation of the potential were utilized to speed up the simulations. First, the truncation of the real space part of the potential for all atoms in a pair of molecules were based on the intermolecular oxygen-oxygen distance, rather than each individual atomic distance. Second, the Coulombic real space term involving the complimentary error function was tabulated into 1E4 bins/Å and interpolated by the Newton-Gregory forward difference method [5]. Monte Carlo trial moves included translation, rotation, deletions and insertions of rigid-body molecules. A translation or rotation move was four times more likely to be attempted than an insertion or deletion move. The parameter associated with the maximum possible translation or rotation was optimized, via a 5% change every 1E6 trials, to yield approximately 25% acceptance of the trial move. Configurational-bias was employed to facilitate insertion and deletion of SPC/E water molecules by performing 10 multiple-first-bead insertions of the Lennard-Jones interactions of an oxygen atom. The Wang-Landau update factor was initially set to unity, and was multiplied by 0.5 whenever the flatness criteria of 0.8 was met. After the update factor was smaller than 5E-5, the collection matrix was updated. After the update factor was smaller than 1E-6, transition-matrix Monte Carlo was performed with an update to the biasing function every 1E6 trials [10]. To parallelize the simulation, each isotherm was divided into 12 overlapping windows to span the entire density range. Because lower density simulations are faster, the parallelization was load balanced by decreasing window size with increasing density by a power scaling with exponent of 1.75. Neighboring windows overlapped by three molecules. The free energy of the entire density range was recovered by setting the free energy of neighboring windows equal at the middle overlapping molecule number, and discarding the largest and smallest number of molecules (when neighbor present) in each window. Simulations were terminated after all windows swept at least 20 times. A simulation is said to have conducted one sweep when each macrostate has been visited from a different macrostate at least 100 times.