SAT-TMMC: Liquid-Vapor coexistence properties - TIP4P Water (LRC)

Liquid-vapor coexistence properties of Water, modeled by the TIP4P Force Field [1], obtained by grand-canonical Wang-Landau/Transition-matrix Monte Carlo and histogram re-weighting. The implementation of the TIP4P model is as described in the SPC/E reference configuration calculation described elsewhere on this website, but with the appropriate TIP4P geometric and interaction parameters [1,2]. Mean values of the saturation pressure, density, and activity (chemical potential- see below) for each phase are reported.

Remarks:

Mean values and standard deviations of the saturation pressure, density, and activity for each phase are reported. These are based upon three independent simulations, initiated from a different random number generator seed. The activity, z, is defined as

$$ \Large z = \Lambda^{-3} \exp \left( \beta \mu \right) $$

where Λ is the de Broglie wavelength, β = 1/(k_BT) (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. 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 [4].

The simulation domain was a cube of side length 20 Å with periodic boundary conditions. The details of the implementation of the potential are as described for the SPC/E reference configuration calculation on this website, but with the appropriate TIP4P geometric and interaction parameters [1,2]. Because the original reference for the TIP4P model used a slightly different definition of the interaction parameters, and different units, it is important to point out that the Lennard-Jones interaction parameters for this data were taken as ε = 0.648694333 kJ/mol and σ = 3.153577942 Å [1]. The Ewald "k" vectors as described in the reference calculation were truncated at k² < 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 [3]. 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 steps, to yield approximately 25% acceptance of the trial move. Configurational-bias was employed to facilitate insertion and deletion of 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 preformed with an update to the biasing function ever 1E6 steps [5]. 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. Phase coexistence was obtained by reweighting the macrostate distribution to the chemical potential where the probability of observing vapor or liquid phases is equivalent.

References

1. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein , J. Chem. Phys., 79, 926-935, 1983.
2. http://www.nist.gov/mml/csd/informatics_research/spce_refcalcs.cfm
3. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, New York, 1989).
4. J. Meija, et al., Pure and Applied Chemistry, 88, 265-291, 2016.
5. V. K. Shen and D. W. Siderius, J. Chem. Phys., 140, 244106, 2014.
6. V. K. Shen and J. R. Errington, JPC B 108, 19595, 2004.
7. V. K. Shen and J. R. Errington, JCP 122, 064508, 2005.
8. V. K. Shen, R. D. Mountain, and J. R. Errington, JPC B 111, 6198, 2007.
9. D. W. Siderius and V. K. Shen, JPC C 117, 5681, 2013.