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

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.


Introduction
The aim of this work is to observe steady-state nucleation events in numerical molecular dynamics (MD) experiments and to evaluate nucleation parameters for droplet nucleation in water vapor. Two water models, SPC/E [1] and TIP4P/2005 [2], are used in this work. These two models have been subjected to several nucleation studies in the past and reference nucleation rate data can be found in the literature for SPC/E [3][4][5][6][7] and also for TIP4P/2005 [8][9][10][11]. Before the introduction of the TIP4P/2005 model, an important study on nucleation using the TIP4P model [12] appeared, in which the Yasuoka-Matsumoto method for evaluation of nucleation rate from MD simulations was introduced. This method was used in the majority of subsequent MD nucleation studies.
Existing MD water nucleation studies differ in the amount of nucleation rate data reported. Simple studies include only one simulation run [9,11], while comprehensive numerical work [7] has collected nucleation rate data over a wide range of temperatures and supersaturations. Accordingly, the computational resources used in the various studies vary by many orders of magnitude. For example, Dumitrescu et al. [10] report 12 simulation runs, each of which consumed about 2.6 thousand CPU hours. On the other hand, Angelil et al. [7] report 18 large-scale simulation runs, the largest of which consumed 8.2 million CPU hours.
An often-discussed feature of the experimental nucleation data is their difference from the theoretical predictions of classical nucleation theory (CNT). In previous works, CNT predictions of the nucleation rate typically differ by many orders of magnitude from the experimental results, which is the case for both MD nucleation simulations of various water models and physical nucleation experiments in laboratory [13]. For nucleation data listed in Table 1 from MD, variations of one to three orders of magnitude are typical for the nucleation rate. The critical cluster size is underestimated by a factor between 2 and 45× . This is similar to the situation with physical nucleation experiments in water vapor, where discrepancies in the nucleation rate have already been quantified in terms of the so-called Wölk correction function [14].
In order to perform a thorough comparison between the MD nucleation rates and the CNT predictions, attention must be paid to ensuring that the nucleation processes are compatible with the assumptions of the CNT statisticalmechanical model. Since the CNT model describes a steadystate nucleation process, i.e., a process of cluster formation in which the thermodynamic state of the parent vapor phase is not affected by the cluster formation, the same conditions must be ensured in the corresponding MD simulations.
In the terminology of MD simulations, such steady-state nucleation processes are realized as NPT processes, where the temperature and pressure are kept constant. Therefore, NPT simulations will be considered in this work, which is the main difference from the previous MD studies of droplet nucleation in water vapor, mostly performed in NVT, but also in NVE ensembles.
The most often used methods for the evaluation of nucleation rate in MD simulations can be characterized as indirect. Typically, evolution of cluster size distribution is being observed during the course of the simulation and the distribution is linked to the nucleation rate by relations originating from nucleation theory equations. In the established Yasuoka-Matsumoto method [12], the nucleation rate is evaluated according to the slopes of number of clusters of a given size as a function of simulation time. Another indirect method for nucleation rate evaluation, the mean-first passage time (MFPT) method [8], measures first appearances in time of certain-sized clusters in the simulation and links the first passages to an error-function-like dependency based on nucleation theory relations. A high number of simulation runs (several hundreds) are being performed to gather reliable statistics for MFPT evaluation.
The approach presented in this work for nucleation rate evaluation lies in a direct observation of nucleation events, i.e., no additional theoretical considerations need to be invoked. Each simulation run is set up to represent a single nucleation event, which is observed after a well-defined period of subcritical cluster size fluctuations in the metastable parent phase, followed by the growth of a supercritical droplet. These two simulation stages need to be clearly distinguishable within each simulation run as the transition between these two stages corresponds to the event of nucleation. The setup of our numerical experiments is described in the "Numerical experiments" section together with the numerical evaluation of MD nucleation parameters, i.e., nucleation rate, impingement rate, critical cluster size, and critical nucleation work. The agreement between nucleation parameters from MD simulations and CNT predictions is discussed in the "Results and discussion" section.

Numerical experiments
Numerical experiments reported in this work have been performed by using the DL_POLY 4.08 molecular dynamics software [15]. For both water models, SPC/E and TIP4P/2005, a cut-off radius of 13 Å was used. Ewald precision parameter for electrostatics was 10 6 . The cut-off radius was chosen to comply with MD simulations performed to evaluate thermophysical properties of SPC/E and TIP4P/2005 water models (Appendices 4 and 5).
Each numerical experiment consists of the following steps: A. Prediction of system pressure at target temperature B. NVT system equilibration at elevated temperature C. NPT evolution of the metastable system at target temperature and pressure The last step C is further divided into two stages: 1. Incubation stage of the metastable vapor 2. Growth of the supercritical liquid droplet The transition between stages 1 and 2 will be considered as the nucleation event, and the time instant when this transition appears will be referred to as the nucleation time.
The first step is to estimate the initial parameters of the simulated system. For a target temperature, we need to estimate the pressure of the system at which a nucleation event can be observed in a reasonable time frame. On the one hand, the nucleation time must fall within a realistic time interval that can be achieved by the numerical simulation, and on the other hand, the nucleation time must be long enough to allow a clear distinction between an incubation stage and a growth stage of the nucleation process. Nucleation times of about 10 ns meet these requirements and were therefore targeted in our simulations.
Each simulation run represents a single nucleation event. For each thermodynamic state, characterized by its temperature and pressure, a total of 20 simulation runs were performed, each starting from different initial positions of water molecules, in order to collect reliable statistics for the evaluation of the average nucleation rate and its standard deviation. Molecular systems of three sizes were simulated, i.e., 108, 500, and 2916 molecules, which resulted in nucleation rates in the range between 10 30 and 10 33 m −3 s −1 . The details of the simulation procedure are described in the following paragraphs. To illustrate our thought process in evaluating the nucleation rate using molecular simulations, an example experiment MD was selected and described in detail. The example system considered in the following paragraphs features 500 TIP4P/2005 molecules at 425 K and 590 kPa.

Classical nucleation theory predictions
We use CNT as a reference for the nucleation rate measurements in molecular simulations. For a given temperature T, CNT is used to find pressure p at which a desired nucleation rate should be observed. The CNT nucleation rate J CNT then represents a reference value for the nucleation rate J MD measured by molecular simulations. We use the standard form of CNT for unary homogeneous droplet nucleation as summarized in Appendix 1. CNT is employed in its basic form without empirical corrections introduced in the past for specific purposes in the nucleation literature. The thermophysical properties of the water models required by CNT are summarized in Appendices 4 and 5, for the SPC/E model and the TIP4P/2005 model, respectively.
The nucleation parameters of our example system as predicted by CNT are listed in Table 2. This tells us that our system with 500 TIP4P/2005 molecules, if allowed to evolve in an NPT simulation at 425 K and 590 kPa, is expected to nucleate one liquid drop after 9.5 ns on average. This corresponds to a steady-state nucleation rate of 2.2×10 31 m −3 s −1 . The technical details of the molecular simulation run are described below.

System equilibration at elevated temperature
An equilibration NVT run at 1000 K is performed at first in order to generate a pool of initial configurations for subsequent NPT nucleation simulation runs. The reason for performing the equilibration run at supercritical temperature is to prevent the formation of larger clusters during this initial simulation stage. An Andersen thermostat was used with 0.5 ps relaxation time and 0.5 softness parameter to control the temperature in this step.
The volume of the NVT system at elevated temperature is adjusted so that the density of the system corresponds to the intended water vapor density at the target simulation temperature and pressure, as given in Table 2 for our example case. This is accomplished by multiplying the saturation density of water vapor at the target simulation temperature by the target supersaturation, which gives us the system density.
The NVT equilibration is performed for 5 ns with a time step of 1 fs, storing the configuration history every 500 time steps. The generated configurations, i.e., atomic positions and velocities, in the period between 1 and 5 ns of the equilibration run are used as a pool for randomly selecting 20 initial configurations to be used in the next simulation step.

Metastable vapor
Each configuration generated at 1000 K is used as the initial configuration in a new NPT simulation run by scaling the velocities to the target simulation temperature (using the restart scale directive of DL_POLY). Since the volume of the system was adjusted in the previous equilibration run to match the density of the water vapor at the target temperature and pressure, a relatively smooth transition from the NVT ensemble to the NPT ensemble is ensured without any significant change in the system volume due to the action of the barostat.
Pressure and temperature are controlled by the Nosé-Hoover thermostat and barostat with relaxation times 0.5 ps and 0.5 ps, respectively. A certain shift in nucleation rates evaluated in this work can be expected depending on the relaxation time parameters. Two pairs of relaxation times were tested to evaluate how the coupling parameters for the thermostat and barostat affect the resulting nucleation rate. For 0.1 ps and 0.1 ps relaxation times, the nucleation time was by 11% shorter than with standard 0.5 ps relaxation time settings for a numerical experiment with 108 TIP4P/2005 molecules at 420 K and 580 kPa. For 2.5 ps and 2.5 ps relaxation times, the nucleation time increased by 30% relative to standard 0.5 ps settings. Such shifts in nucleation times represent changes in nucleation rate well within the standard deviation (one-sigma) of the nucleation rate estimation by our method, and therefore the effect of relaxation time parameters can be considered as negligible. The time step in our NPT simulations is 2 fs. The evolution of the metastable system is monitored by storing the positions of the atoms every 1000 time steps. As the metastable system is free to develop larger clusters over time, a cluster monitoring procedure was developed to calculate the number of clusters and their sizes in each of the stored configurations. The cluster search algorithm uses the Stirling criterion with 2.5 Å O-H distance [9].

Growth of the liquid cluster
The transition of the system to the growth stage is observed, when the size of the maximum cluster exceeds the critical size and progresses in its growth until the simulated system consists entirely of the liquid phase. Figure 1 shows the evolution of the maximum cluster size in a single NPT run of our example simulation.
In Figure 1, a clear distinction can be observed between the initial incubation stage, which lasts for about 10 ns, and the subsequent growth stage. The transition to continuous growth of the maximum cluster marks the nucleation event. During the incubation stage, the system resides in a metastable state for a period of time and is characterized by fluctuations in the maximum cluster size, which can even exceed the critical cluster size slightly, but they break up. Also during the incubation stage, the system fluctuates around a constant volume as can be seen in Figure 2. The growth stage is characterized by a continuous growth of the maximum cluster in the system. During the growth stage, the maximum cluster is scavenging the gas molecules and small clusters from the system. Since the system is maintained at a constant temperature and pressure by the action of the thermostat and barostat, the transition of molecules from the vapor phase to the liquid phase manifests itself by a considerable decrease in the volume of the system, as can be seen in Figure 2. Similarly, the density of both monomers and small clusters decreases as the system nucleates. This can be observed from the number of monomers in Figure 3 and the number of clusters with sizes in the range of two to nine molecules in Figure 4.
To calculate the time of nucleation t nuc in our simulation run, we iteratively search for the transition between the incubation stage and the growth stage as follows. First, we make a rough estimate of the nucleation time t nuc,0 by considering the nucleation time as the time when the maximum cluster becomes larger than 20% of the total system molecules. In the first iteration, an average of the maximum cluster sizes n inc,0 for the incubation period from the simulation start to t nuc,0 is evaluated, and a cubic fit c gr,0 of the maximum cluster sizes for the growth period between t nuc,0 and the simulation end is evaluated as well. A second estimate for the nucleation time t nuc,1 is calculated as the time, at which the cubic fit of cluster growth equals the average cluster size during the incubation period c gr,0 (t nuc,1 ) = n inc,0 . The new estimate of the nucleation time t nuc,1 is used in the second iteration to evaluate the average cluster size n inc,1 during the incubation stage and the cubic growth function c gr,1 during the growth stage, etc. After several iterations, the search for the nucleation time is finished under the condition |t nuc,i+1 − t nuc,i | < 0.01 ns. The nucleation time estimated by this numerical method is shown in Figures 1, 2, 3, and 4 by the dotted vertical line at 9.8 ns, and the final cubic fit of the maximum cluster growth is shown by the orange line in Figure 1.

Cluster growth rate
The growth of the cluster beyond the critical cluster size is characterized by a power-law dependence of the size on time. Assuming that the impingement rate of the vapor molecules on the cluster surface is constant, the size evolution of the cluster is represented by a cubic function as described in Appendix 2. This behavior can be observed in Figure 1 (orange line).
The attachment of the vapor molecules to the growing cluster can be evaluated in the MD experiment by comparing the maximum cluster size between successive configurations stored during the numerical simulation. The instantaneous growth rate calculated in this way is plotted in Figure 5 (blue line). After fitting this stochastic growth process to the cubic model according to Appendix 2, a good agreement is observed between the resulting average impingement rate MD (green dotted line) and the impingement rate model used by CNT (orange line) in Figure 5.

Nucleation rate
Since the nucleation event is a stochastic phenomenon, variations in the nucleation time of parallel NPT simulation runs starting from different initial configurations are to be expected, even though all simulation runs have the same target pressure and temperature. For the 20 parallel simulation runs of our example system, the nucleation times shown in Figure 6 span from a few nanoseconds to over 100 ns. Note that such a wide range of nucleation times shows that The lognormal average nucleation time t nuc for this set of simulation runs is 23.8 ns a nucleation rate measured in a single-run MD simulation [9,11] is a poor estimate of the steady-state nucleation rate.
To calculate the average nucleation time t nuc , a log-normal distribution of nucleation times is assumed. A log-normal distribution of nucleation times corresponds to normal distribution of the critical nucleation work due to relation (A4). In our example case, the log-normal average nucleation time is 23.8 ns, where the standard deviation (one-sigma) corresponds to the interval (11.1, 51.0) ns.
After estimating the nucleation time, the nucleation rate is evaluated as where V avg is the average volume of the system during the incubation stages (Figure 2) of the set of 20 simulation runs under evaluation. A well-established incubation phase during the NPT simulation runs is essential for a reliable measurement of V avg . The standard deviation of the average volume in the 20 simulation runs of our example case was 0.5%.
The resulting nucleation rate according to Eq. (1) in our example system is 1.1×10 31 m −3 s −1 with standard deviation between 5.1×10 30 and 2.4×10 31 . We see that the original nucleation rate estimate by CNT of 2.2×10 31 m −3 s −1 ( Table 2) falls within the standard deviation of our MD nucleation rate measurement in this example case.

Nucleation work
A standard approach to extracting the nucleation work ΔG MD from the cluster size distribution is used in this work, as described in Appendix 3. Figure 7 shows the nucleation work calculated from the cluster distributions for our (1) J MD = 1 V avg t nuc example system according to Eq. (C11) compared to the CNT nucleation work calculated according to Eq. (A1).
The vertical lines in Figure 7 represent the critical cluster sizes. For CNT (dotted orange line), the critical size is calculated according to Eq. (A3). For the set of 20 NPT simulations, the nucleation work data (blue points) are fitted to a CNT-like nucleation work function ΔG MD (n) = p 1 n 2∕3 − p 2 n , and the maximum of the fitted function is calculated analytically (dashed blue line).

Results and discussion
According to the above described procedure, nucleation parameters of two water models, SPC/E and TIP4P/2005, were evaluated from our MD simulations and compared to CNT predictions.
For SPC/E, nucleation simulations for temperatures in the range 300-425 K were performed for systems of 108 molecules in order to monitor the temperature dependence of the nucleation rate. At 375 K, additional system sizes of 500 and 2916 molecules were simulated in order to observe the pressure dependence of the nucleation rate, because larger systems at lower supersaturations allow for the simulation of lower nucleation rates.
For TIP4P/2005, simulations of 108 molecules at temperatures in the range 325-475 K were run. At 425 K, system sizes of 500 and 2916 molecules were simulated as well.

Nucleation rate
The nucleation rates measured in SPC/E simulations are plotted in Figure 8 and the corresponding data points are summarized in Table 3 together with uncertainty intervals (one-sigma in the logarithmic scale). We can see that the At 375 K, SPC/E simulations of larger systems allowed us to evaluate the pressure dependence of the MD nucleation rate. The slope of the 375 K CNT isotherm in Figure 8 is again reproducing the slope of MD nucleation rate data quite accurately.
The comparison of our SPC/E results and nucleation rates reported by Angelil et al. [7] is shown in Figure 9. Both data sets are reproduced by CNT accurately. For Angelil et al.'s data, we actually found a closer agreement with CNT than the authors reported. This is caused by the slightly different parametrizations of SPC/E thermophysical properties employed in the respective CNT calculations, as the accurate NIST SPC/E thermophysical properties [16] were not available at the time of Angelil et al.'s work.
The other SPC/E NVT nucleation rates reported in the literature (Table 1) deviate from CNT by several orders of magnitude. Therefore, the agreement between the CNT prediction and Angelil et al. NVT data might seem coincidental. However, since Angelil et al. performed large-scale molecular simulations (up to 4.1 million SPC/E molecules), the pressure during their NVT simulation runs was not decreasing substantially due to droplet formation and as such the simulated process was resembling the steady-state nucleation process also quite closely.
The nucleation rates measured in TIP4P/2005 simulations are plotted in Figure 10, and the data points are given in Table 4. The agreement between CNT predictions and MD nucleation rates is again very good. For 5 out of the total 9 cases, the CNT nucleation rates lie within the uncertainty interval of MD nucleation rate, and in the other cases the   Figure 9 Nucleation rate in SPC/E water performed in this work (dots) compared with nucleation rates measured in NVT simulations by Angelil et al. [7] (crosses) and with CNT predictions (lines) Figure 10 Nucleation rate in TIP4P/2005 water by numerical experiments (dots) vs. nucleation rate predictions by CNT (lines) CNT nucleation rate differs from the MD rates by a factor of less than 3.1× . The worst case shows a difference of 8.56× 10 31 m −3 s −1 CNT nucleation rate vs. 2.78×10 32 m −3 s −1 MD nucleation rate. The level of agreement between CNT predictions and our MD measurements is basically the same for both water models. In both cases, the best agreement can be observed for high temperatures, i.e., for simulations with larger critical clusters.
To justify our choice of 20 simulations as sufficient statistics for averaging the nucleation rate, or, in other words, to illustrate the precision of evaluating the nucleation rate from sets of 20 parallel nucleation simulations, nucleation rates were evaluated for 6 sets of 20 simulations at the same thermodynamic state of the parent vapor, i.e., 450 K and 1030 kPa, as shown in Figure 11. The data points differ by a factor less than 2 × in the absolute value of the nucleation rate, while the respective uncertainty intervals overlap to a large extent. The differences in the temperature during the incubation stages of the numerical experiment measured as an average over the set of 20 simulation runs are of the order of 0.01 K and the differences in pressure measured in the same way show an order of 1 kPa.

Impingement rate
The impingement rate evaluated from our numerical experiments according to the procedure described in the "Cluster growth rate" section is summarized in the last columns of Tables 3 and 4, for SPC/E and TIP4P/2005 models, respectively. A very good agreement with CNT impingement rates calculated according to the basic ideal-gas model (A7) is observed. For SPC/E, the difference in the worst case is by a factor of 2.6×.   For TIP4P/2005, the worst difference factor is 7.5× at 475 K, which is the highest simulated temperature. It should be noted, however, that the molecular system at 475 K is on the edge of its ability to spawn the steady-state nucleation process properly. The critical cluster size for 475 K simulation runs is almost 30, while the total number of molecules in the simulation is 108 only. In such case, the size fluctuations of the largest cluster around the critical size may reach the total size of the system. In other words, the growth/decay behavior of the nucleating cluster around its critical size may feature situations, when the cluster interacts with itself, because there is not enough vapor molecules surrounding it, effectively altering the nucleation process.

Nucleation work and critical cluster size
Critical nucleation work and the corresponding critical cluster size have been evaluated from molecular simulations according to the procedure described in the "Nucleation work" section.
For SPC/E, the results are summarized in Table 5. The critical clusters range in size from 5.8 to 15.5 molecules in our MD experiments. Smaller critical sizes have been achieved for lower system temperatures, which reflects the way of selecting the system temperatures and pressures in order to arrive at sufficiently long nucleation times in the respective simulations. Slight underprediction (of up to 40 %) of the critical cluster size by CNT can be observed for temperatures below 375 K and for smaller critical clusters. For higher temperatures, the critical cluster size is being slightly overpredicted (less than 14 %) by CNT. The critical nucleation work ΔG ⋆ MD is reproduced by CNT quite accurately over the whole range of SPC/E simulations.
For TIP4P/2005, critical cluster parameters are summarized in Table 6, showing a range of critical clusters sizes from 6.0 to 33.2 molecules. We can observe a similar level of agreement between numerical results and CNT predictions as in the case of SPC/E simulations; the differences between theoretical critical cluster size and MD experiment is less than 2 molecules. The critical nucleation work is also reproduced by CNT quite accurately over the range of simulated TIP4P/2005 systems. A relatively larger deviation in both nucleation work and cluster size can be observed only in the 475 K case, which may represent a system in which the nucleation process deviates from the desired steady-state scenario as already mentioned above.

Conclusions
Droplet nucleation processes in water vapor were simulated using molecular dynamics (MD). Two molecular models for water were selected, SPC/E and TIP4P/2005, which have Table 5 SPC/E critical cluster size n ⋆ and critical nucleation work ΔG ⋆ predicted by CNT vs. measured in MD simulation Nucleation rates were measured directly in the molecular simulations by observing the occurrence of the supercritical cluster. The agreement between the MD simulations and the predictions of classical nucleation theory (CNT) was found to be approximately in the range of standard deviation of the MD nucleation rate. A similar level of agreement was found between the impingement rate measured for the growing, just-nucleated cluster and the impingement rate calculated using the simple, ideal-gas-based model adopted by CNT. Two other nucleation parameters, i.e., the critical nucleation work and the critical cluster size, were derived from the MD cluster distributions according to classical theoretical relations between the steady-state cluster distribution and the constrained-equilibrium cluster distribution. These two nucleation parameters agree with the CNT predictions with deviations in the low double-digit percentage range, even for critical cluster sizes in the single-digit range.
Compared to the MD simulations of water nucleation previously presented in the literature, our MD simulations are distinguished by their adherence to the constant temperature and constant pressure assumptions of the CNT model. This cannot be emphasized enough; the CNT model predicts a steady-state nucleation rate and as such it has been derived assuming constant temperature and constant pressure of the parent phase during the nucleation process. Therefore, attempts to reproduce nucleation rates using CNT for processes other than NPT are bound to be inaccurate. The degree of inaccuracy of such attempts will be proportional to the degree of deviation of the simulated nucleation process from the steady-state nucleation case.
There has been a considerable progress on developing more accurate water models in the last decades [17]. Several models, e.g., OPC, OPC3, TIP3P-FB, and TIP4P-FB, feature a closer agreement with thermophysical properties of water. Future nucleation simulations with such more advanced molecular models may result in nucleation rates corresponding to those from laboratory nucleation experiments in water vapor. where the pre-exponential factor J 0 (m −3 s −1 ) is evaluated as which is a product of the non-equilibrium (Zeldovich) factor Z, impingement rate of monomers onto the cluster surface (m −2 s −1 ), surface area of the critical cluster A ⋆ (m 2 ), and water vapor number density v (m −3 ), which is calculated as a product of supersaturation ratio p∕p sat and saturation vapor density sat,v (T) (m −3 ).
The Zeldovich factor in Eq. (A5) takes the form and the impingement rate (m −2 s −1 ) is calculated as 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 (T) • Saturation pressure p sat (T) • Saturated vapor number density sat,v (T) • Liquid number density l (T) The above thermophysical properties of the SPC/E and TIP4P/2005 water models are summarized in Appendices 4 and 5, respectively. from cluster distribution data The procedure to evaluate the nucleation work Δ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 where the constrained-equilibrium distribution eq is related to the steady-state distribution through a recursive relation [18] The steady-state distribution (n) is the actual cluster distribution observed in MD simulation and is therefore calculated 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.