Nucleation parameters of SPC/E and TIP4P/2005 water vapor measured in NPT molecular dynamics simulations

Abstract

Nucleation rates for droplet formation in water vapor are measured in molecular dynamics (MD) simulations of SPC/E and TIP4P/2005 water by monitoring individual nucleation events. The nucleation process is simulated in the NPT ensemble to evaluate the steady-state nucleation rate in accordance with the assumptions of classical nucleation theory (CNT). Nucleation rates measured between 300 and 425 K for the SPC/E model, and between 325 and 475 K for the TIP4P/2005 model, agree with the CNT predictions roughly within the standard deviation of the MD measurements of the nucleation rates.

Appendices

Appendix 1. Classical nucleation theory overview

The following Eqs. (A1)–(A7) summarize the formulas of the classical nucleation theory (CNT) for the case of unary homogeneous nucleation of droplets. A thorough derivation of the formulas governing the nucleation parameters can be found in nucleation textbooks, e.g., Kalikmanov [18]. The nucleation parameters reflect the thermophysical state of the nucleating parent phase characterized by pressure p and temperature T. Within the CNT model, these two parameters are assumed constant during the nucleation process, and the required thermophysical properties of water are being evaluated at the nucleation pressure and temperature.

The nucleation work $\Delta G$ (J) within CNT is a function of the cluster size given by the number of molecules n

$$\begin{aligned} \Delta G(n) = \sigma \root 3 \of {36 \pi (v_l n)^{2}} - n k T \ln \frac{p}{p_{sat}} \end{aligned}$$

(A1)

where $v_l$ (m$^{3}$) is molecular volume of water molecules in the liquid phase, $p_{sat}$ (Pa) is the saturation pressure of water, and $\sigma$ (N/m) is the surface tension.

The energy barrier to nucleation denoting the maximum of the nucleation work (A1), i.e., the critical nucleation work $\Delta G ^{\star }$, takes the form

$$\begin{aligned} \Delta G ^{\star } = \frac{16 \pi }{3} \frac{\sigma ^3}{\left( k T \ln \frac{p}{p_{sat}} \right) ^2} \end{aligned}$$

(A2)

The critical cluster size, conveniently expressed as the number of molecules in the critical cluster $n^{\star }$, is calculated as

$$\begin{aligned} n^{\star } = 36 \pi v_l^2 \left( \frac{2}{3} \frac{ \sigma }{k T \ln \frac{p}{p_{sat}}} \right) ^3 \end{aligned}$$

(A3)

The nucleation rate J (m$^{-3}$s$^{-1}$), representing the number of droplets crossing the critical size and growing to macroscopic proportions, is given by

$$\begin{aligned} J = J_0 \exp \left( - \frac{\Delta G ^{\star }}{kT} \right) \end{aligned}$$

(A4)

where the pre-exponential factor $J_0$ (m$^{-3}$s$^{-1}$) is evaluated as

$$\begin{aligned} J_0 = Z \beta A^{\star } \rho _v \end{aligned}$$

(A5)

which is a product of the non-equilibrium (Zeldovich) factor Z, impingement rate of monomers onto the cluster surface $\beta$ (m$^{-2}$s$^{-1}$), surface area of the critical cluster $A^{\star }$ (m$^2$), and water vapor number density $\rho _v$ (m$^{-3}$), which is calculated as a product of supersaturation ratio $p/p_{sat}$ and saturation vapor density $\rho _{sat,v}(T)$ (m$^{-3}$).

The Zeldovich factor in Eq. (A5) takes the form

$$\begin{aligned} Z = \frac{\root 6 \of {36}}{3} \root 3 \of { \frac{v_l }{\pi }} \sqrt{ \frac{\sigma }{k T}} (n^{\star })^{-2/3} \end{aligned}$$

(A6)

and the impingement rate $\beta$ (m$^{-2}$s$^{-1}$) is calculated as

$$\begin{aligned} \beta = \frac{p}{\sqrt{2 \pi m_1 k T}} \end{aligned}$$

(A7)

To sum up, CNT allows us to predict a steady-state nucleation rate in metastable water vapor, which is a kinetic property of the macroscopic parent system at a pressure p and temperature T. The CNT relations utilize the following four thermophysical properties of macroscopic water substance:

• Surface tension $\sigma (T)$

• Saturation pressure $p_{sat}(T)$

• Saturated vapor number density $\rho _{sat,v}(T)$

• Liquid number density $\rho _{l}(T)$

The above thermophysical properties of the SPC/E and TIP4P/2005 water models are summarized in Appendices 4 and 5, respectively.

Appendix 2. Impingement rate evaluation from cluster growth data

The procedure to deduce the average impingement rate of vapor molecules attaching onto the cluster surface from molecular simulation data is based on the calculation of instantaneous growth and decay of the maximum cluster during the cluster growth stage. Assuming a constant impingement rate $\beta _{MD}$ (m$^{\mathrm {-2}}$ s$^{\mathrm {-1}}$) and a single supercritical growing cluster, the cluster growth can be described as a change of the number of cluster molecules n over time t in the following manner

$$\begin{aligned} \frac{d n}{d t} = \root 3 \of {36 \pi (v_l n)^{2}} \beta _{MD} \end{aligned}$$

(B8)

where the right-hand side of Eq. (B8) represents the cluster surface area multiplied by the impingement rate. Differential equation (B8) can be readily integrated and its solution takes a simple form

$$\begin{aligned} n = C_0 (t - t_{0})^3 \end{aligned}$$

(B9)

Finally, the parameters $C_0$ and $t_0$ can be evaluated by fitting the cubic function (B9) to the cluster growth data from MD simulation, i.e., number of cluster molecules vs. simulation time (as shown in Figure 5, blue line), during the growth stage of the simulation run. The fitted value of $C_0$ is then related to the impingement rate $\beta _{MD}$ as

$$\begin{aligned} \beta _{MD} = \frac{C_0}{\root 3 \of {36 \pi v_l^{2}} } \end{aligned}$$

(B10)

Appendix 3. Nucleation work evaluation from cluster distribution data

The procedure to evaluate the nucleation work $\Delta G_{MD}$ as a function of cluster size n from cluster size distribution has been described in previous MD nucleation studies of water nucleation [3, 9, 12]. In this work, the cluster size distribution is calculated from saved configurations of all simulation runs in the MD simulations set, which in our case consisted of 20 simulation runs. For cluster size n, the nucleation work is a function constrained-equilibrium distribution of clusters

$$\begin{aligned} \frac{\Delta G_{MD} (n)}{k T} = - \ln \frac{\rho _{eq}(n)}{\rho _{eq}(1)} \end{aligned}$$

(C11)

where the constrained-equilibrium distribution $\rho _{eq}$ is related to the steady-state distribution $\rho$ through a recursive relation [18]

$$\begin{aligned} J_{MD} \sum \limits _{m=1}^{n-1} \frac{1}{\beta _{MD} A(m) \rho _{eq}(m)} = \frac{\rho (1)}{\rho _{eq}(1)} - \frac{\rho (n)}{\rho _{eq}(n)} \end{aligned}$$

(C12)

The steady-state distribution $\rho (n)$ is the actual cluster distribution observed in MD simulation and is therefore calculated

$$\begin{aligned} \rho (n) = \frac{1}{20} \sum \limits _ {i=1}^{20} \frac{1}{C_i} \sum \limits _ {j=1}^{C_i} \frac{c_{ij}(n)}{V_{MD,ij}} \end{aligned}$$

(C13)

where i denotes the simulation run and j denotes the j-th configuration within run i out of its total $C_i$ saved configurations. Furthermore, $c_{ij}(n)$ is the number of clusters of size n found in configuration j of run i, and $V_{MD,ij}$ is the volume of the simulation cell in configuration j of run i.

Note that to solve Eq. (C12) for $\rho _{eq}(n)$, two of the nucleation parameters, $J_{MD}$ and $\beta _{MD}$, need to be known in advance. Indeed this is the case, since both parameters are measured in our MD simulations independently. However, any errors in measuring $J_{MD}$ or $\beta _{MD}$ will inevitably lower the accuracy of evaluating $\Delta G_{MD}$.

Appendix 4. Thermophysical properties of SPC/E water

Parametrizations of thermophysical properties of SPC/E water valid for the temperature range 300–425 K are summarized below.

Surface tension $\sigma$ (N/m) was parametrized by Vega and Miguel [19] as

$$\begin{aligned} \sigma (T) = s_1 \left( 1 - \frac{T}{T_c} \right) ^{11/9} \left( 1 - s_2 \left( 1-\frac{T}{T_c}\right) \right) \end{aligned}$$

(D14)

where $s_1$ = 0.20532, $s_2$ = 0.6132, and $T_c$ = 625.7.

Saturation pressure $p_{sat}$ (Pa) data reported by NIST [16] have been fitted to an Antoine-type of function

$$\begin{aligned} p_{sat} (T) = 10^5 \exp \left( p_1 + \frac{p_2}{T + p_3} \right) \end{aligned}$$

(D15)

where $p_1$ = 12.74080, $p_2$ = −4712.723, and $p_3$ = −25.67480.

Saturated vapor density $\rho _{sat,v}$ (m$^{-3}$) data reported by NIST [16] have been fitted to a polynomial function in exponent as

$$\begin{aligned} \rho _{sat,v} (T) = \exp \left( \sum \limits _{i=0}^4 \frac{a_i}{T^i} \right) \end{aligned}$$

(D16)

where $a_0$ = 2.855$\times$10$^{1}$, $a_1$ = −2.897$\times$10$^{4}$, $a_2$ = 1.304$\times$10$^{7}$, $a_3$ = −3.116$\times$10$^{9}$, and $a_4$ = 2.722$\times$10$^{11}$.

Liquid density $\rho _{l}$ (m$^{-3}$) data reported by NIST [16] have been fitted to a polynomial function

$$\begin{aligned} \rho _{l} (T) = \sum \limits _{i=0}^4 \frac{b_i}{T^i} \end{aligned}$$

(D17)

where $b_0$ = −4.686$\times$10$^{3}$, $b_1$ = 7.169$\times$10$^{6}$, $b_2$ = −3.533$\times$10$^{9}$, $b_3$ = 7.954$\times$10$^{11}$, and $b_4$ = −6.815$\times$10$^{13}$.

Appendix 5. Thermophysical properties of TIP4P/2005 water

Parametrizations of thermophysical properties of TIP4P/2005 water valid for the temperature range 325–475 K follow below.

Surface tension $\sigma$ (N/m) was parametrized by Vega and Miguel [19] as Eq. (D14), where $s_1$ = 0.22786, $s_2$ = 0.6413, and $T_c$ = 641.4.

Saturation pressure $p_{sat}$ (Pa) data reported by Vega et al. [20] have been fitted to an Antoine-type of function (D15), where $p_1$ = 12.46120, $p_2$ = −4476.552, and $p_3$ = −41.49840.

Saturated vapor density $\rho _{sat,v}$ (m$^{-3}$) data reported by Vega et al. [20] have been fitted to function Eq. (D16), where $a_0$ = 3.578$\times$10$^{1}$, $a_1$ = −3.935$\times$10$^{4}$, $a_2$ = 1.825$\times$10$^{7}$, $a_3$ = −4.216$\times$10$^{9}$, and $a_4$ = 3.526$\times$10$^{11}$.

Liquid density $\rho _{l}$ (m$^{-3}$) data reported by Vega et al. [20] have been fitted to function Eq. (D17), where $b_0$ = −5.385$\times$10$^{3}$, $b_1$ = 8.130$\times$10$^{6}$, $b_2$ = −3.958$\times$10$^{9}$, $b_3$ = 8.6345$\times$10$^{11}$, and $b_4$ = −7.067$\times$10$^{13}$.