Thermodynamics and the potential energy landscape: case study of small water clusters

Abstract

We investigated the structure and the thermodynamic properties of small water clusters with the nested sampling computational technique, using two different water models, the coarse-grained mW (up to 25 molecules) and the flexible version of the TIP3P (up to 16 molecules). By mapping the entire potential energy landscape of the clusters, we calculated the heat capacity curves, located the structural transitions and identified those local minimum basins which contribute the most to the total partition function. We found that in the case of the mW model, trends in first-order-like and continuous-like transitions can be very well matched to the characteristics of the landscape: cluster sizes with fewer and narrower local minimum basins show a sharper ‘melting’ peak on the heat capacity curve. Trends in the case of the TIP3P model were not easily assigned to the changing occupation of basins, and the contribution of local minima was negligible, except for n = 7, 15 and 16.

---

1 Introduction

Understanding the structural and thermodynamic properties of small water clusters has been in the focus of both experimental and theoretical research for several decades. Water clusters play a crucial role in atmospheric processes,¹ affect the geometry and stability of macromolecules,² and as small building blocks they have importance in helping us to understand the anomalous behaviour of bulk water.^3,4 Moreover, in the field of computational chemistry, water clusters serve as a testing case for global optimisation and sampling techniques, as they pose a significantly harder challenge than atomic systems with similar numbers of particles, because of the translational and orientational degrees of freedom of the molecules. Due to the increased computational cost, the majority of theoretical studies concentrate on water clusters consisting of fewer than 25 molecules, and on determining the global minimum structure only. Such studies have been performed with a variety of methods (e.g. with simulated annealing,⁵ basin hopping algorithm⁶ and genetic algorithm⁷) and for a large number of different potential models of water, such as the coarse-grained model of mW,⁶ TIP3P,^7–9 different parametrisations of TIP4P,^7,9,10 TIP5P,¹¹ SPC/E⁷ and ab initio methods.¹² These works have shown that different potential models often predict significantly different configurations as the global minimum structure.^9,11 The water hexamer is particularly controversial and considerable effort has been invested into studying the ground state of (H₂O)₆ both computationally and experimentally.^13–19

However, to be able to gain a better understanding of the physical processes and the behaviour of a system, a more extensive search of the potential energy landscape is necessary, determining thermodynamic properties and local minimum configurations at finite temperatures. A number of local minima have been mapped for a few smaller cluster sizes,²⁰ and the disconnectivity graph of (H₂O)₂₀ has been constructed with the TIP4P potential.^21,22 Although a smaller number of minima were located than in the case of atomic Lennard-Jones clusters with the same number of degrees of freedom, the landscape was found to be frustrated, with competing low-energy minima of different morphologies, separated by high barriers. Heat capacity curves of several water models have been calculated for small water clusters, using parallel tempering calculations,^20,23 Wang–Landau sampling^24,25 and multicanonical-ensemble molecular dynamics.²⁶ These studies concluded that the thermodynamic properties of the clusters strongly depend on the potential model, and while certain cluster sizes show a sharp peak on the heat capacity curve, others have only a shoulder representing a change. This suggests that certain cluster sizes (e.g. (H₂O)₈ and (H₂O)₁₂) undergo a first-order-like melting transition,⁵ while in other cases (e.g. (H₂O)₇ and (H₂O)₁₁) continuous-like melting occurs.²⁶ It has to be noted that, although phase transitions are defined for macroscopic systems, the terms melting and evaporation are routinely used in finite systems as well, to describe the temperature region where the cluster starts to explore a much larger phase space, similarly to a phase transition.^27,28

In the present work we used the nested sampling technique to sample the potential energy landscape of small water clusters, using two different potential models: the coarse-grained model mW²⁹ and the flexible version of the TIP3P³⁰ potential model. Nested sampling not only allows us to calculate the partition function, and thus the thermodynamic properties such as the heat capacity, but also gives us a unique insight into the structure of the landscape, allowing us to identify the basins with significant thermodynamic contribution, and determine their relative phase space volume.

2 Computational details

2.1 Nested sampling

Nested sampling (NS) is a Bayesian statistical approach^31,32 adapted to explore the atomic phase space.³³ It is a top-down sampling technique, starting from the gas phase and progressing towards the ground state structure through a series of nested energy levels, estimating the corresponding phase space volume of each. This way the method samples the different basins proportional to their volume, and instead of providing an exhaustive list of the local minima, it identifies the thermodynamically most relevant states without any prior knowledge of the structures or phase transitions. The NS algorithm allows the direct computation of the partition function, and thus gives access to thermodynamic properties, such as the heat capacity or the free-energy. The NS method has been applied to study Lennard-Jones clusters^33,34 and the hard-sphere model,³⁵ and it has been shown that it enables the automated calculation of the complete pressure–temperature–composition phase diagram.^36–38

The nested sampling calculations were performed as presented in ref. 33, using the pymatnest program package.^36,39 All simulations were run at constant volume, chosen such that the density is 0.00032 water molecules per ų (0.0095 g cm⁻³) for every cluster size in the case of TIP3P, and 0.000625 water molecules per ų for mW. The initial configurations were generated randomly in the case of the mW model, and by 6 × 10⁵ Hamiltonian Monte Carlo³⁷ (all-atom) moves for every sample in the case of the TIP3P model. New samples were generated by performing Hamiltonian Monte Carlo moves for both potential models, using the LAMMPS package⁴⁰ for the dynamics. The number of walkers was chosen such that the difference in the melting temperature predicted by independent parallel runs is less than 10 K, and the solid–solid transition is found reliably (Table 1).

Table 1 Typical nested sampling parameters used in the case of different cluster sizes: the number of walkers, K, which determines the resolution of the sampling, and the length of the walk used for generating new configurations, L. The table also includes the number of total energy evaluations in a single run, N_e

(H2O)n	K	L	N e
mW
5–10	1000	600	3 × 10^8
11–15	2500	800	6 × 10^8
16–19	4000	800	1 × 10^9
20–25	5000	800	2 × 10^9
TIP3P
5–10	2000	4000	5 × 10^9
11–15	6000	6000	3 × 10^10
16	11520	6000	6 × 10^10

---

2.2 Potential models

The first water model we chose is mW,²⁹ to explore the behaviour of particles preferring tetrahedral arrangement without electrostatic interactions, thus being computationally more affordable. mW is a simple coarse-grained model of water, reparametrised from the Stillinger–Weber potential⁴¹ to reproduce the density and melting temperature of water at atmospheric pressure. Despite having only short-range interactions and describing the hydrogen bonded network through the angular dependent term without electrostatic interactions, it reproduces the energetics, anomalies, and the liquid and hexagonal ice structure of water remarkably well. Although the mW model is practically a single atom, for the ease of readability we will refer to a single mW site as a molecule throughout the paper.

Our second choice for the water model was TIP3P, as it is expected to show a larger dependency in thermodynamic behaviour on cluster size than other TIPnP models,^24,25 thus making it an ideal case for testing. The TIP3P water model is a three-site description of the water molecule,³⁰ where both hydrogen and oxygen atoms carry a charge, and a single Lennard-Jones interaction site is placed on the oxygen. Details of the flexible model are shown in Table 2. We used a cutoff distance of 12.0 Å to truncate the Lennard-Jones interactions, and the Coulombic interaction was calculated explicitly within a cutoff of 0.6 × (box length), outside which the Ewald summation was used.

Table 2 Potential parameters for the flexible version of the TIP3P water model

Atom	q/e	ε/eV	ε/Å
Hydrogen	0.415	0.00000	0.00
Oxygen	−0.830	0.00660	3.15

---

K/eV rad^−2	θ 0/°	K/eV rad^−2	r 0/Å
2.39	104.52	19.51	0.9572

---

3 Results

3.1 mW clusters

Heat capacity curves calculated for clusters with the mW potential model are shown in Fig. 1. The most prominent peaks represent the gas–liquid transition, which depends on the density. There is a slight trend in the evaporation transition shifting to higher temperatures as the cluster size increases, as expected. At temperatures between 200 and 300 K a smaller peak can be observed, corresponding to the solidification of the cluster, shifting towards lower temperatures as the size of the cluster increases. For sizes n = 8, 12, 14, 16, 19, 21, 23 and 25 this peak is sharper, indicating a first-order-like transition, while for the rest of the clusters it only appears as a shoulder. In two of the clusters, n = 7 and 20, an additional small peak can be found (see the inset of Fig. 1), suggesting an additional solid–solid transition at very low temperatures.

Fig. 1 Heat capacity of (mW)_n clusters as a function of temperature for 5 ≤ n ≤ 25. For every cluster size the results of at least two independent runs are shown. The curves are vertically shifted for ease of readability. The inset shows the low temperature results for clusters n = 7 and 20.

In order to distinguish between different structures found during the sampling, we calculated the Steinhardt bond-order parameters Q₄ and Q₆,⁴² and the average size of the rings formed within the cluster (using 3.45 Å as the cutoff distance), for all configurations generated during the nested sampling calculation.

With (mW)₁₃ as an example, Fig. 2 demonstrates this; showing the order parameters as a function of temperature and the two transitions, which are clearly visible. Below T ≈ 550 K, rings are being formed and the variance in the order parameter decreases as the molecules condense and create a cluster. Below T ≈ 280 K, where the shoulder on the heat capacity function can be observed, two groups of points become clearly distinguishable, corresponding to two different structures, the global minimum and one of the local minima. This indicates that although a large number of other molecular arrangements must be possible, these two configurations have far the largest phase space volume and thus thermodynamic relevance.

Fig. 2 Average Q₄ bond order parameter as a function of temperature of configurations generated during nested sampling for the (mW)₁₃ cluster. Each dot corresponds to a configuration and is coloured according to the average ring size within the cluster. Temperatures corresponding to the peaks on the heat capacity curve are shown by vertical dotted lines.

Using this approach, we are able to identify different basins and assign which configurations belong to these, and then use this information to calculate the partition function, Z, of the different basins separately, as

where β is the inverse thermodynamic temperature, K is the number of walkers, and U is the energy of the i-th configuration that belongs to basin A. This allows us to determine the partition function ratios and thus the structures most populated under certain conditions. Partition function ratio results are shown in Fig. 3, together with the structures of the global and local minima identified. The plots only show results of individual runs, however parallel calculations found the same minimum basins in each case, with their relative contribution to the partition function not differing more than approx. 10% of the ones shown in the figure.

Fig. 3 Partition function ratio of different basins (Q_basin/Q_total) as a function of temperature for water clusters with the mW potential. Dark blue corresponds to the global minimum structure; orange, green, light blue and purple are different local minima, while grey corresponds to basins of smaller volume or structures difficult to assign to well defined basins (e.g. melt-like configurations). Snapshots of the structures and corresponding energies are also shown.

Global minima found by nested sampling are all identical to those published previously by Farrell et al.⁶ At very low temperatures these are the most relevant structures, as expected, but as the temperature increases the contribution of certain local minima can increase significantly, depending on the size of the cluster. For n < 18, clusters with an even number of molecules, the contribution of local minima to the partition function is negligible, while for odd numbers of molecules at least one of the local minima becomes the most dominant structure at a relatively low temperature. This trend turns to its opposite once one of the molecules appears to be in the middle of the cluster, thus solvated, the first such ground-state being n = 19. It is interesting to note that the first solvated molecule appears in the case of the highest energy local minima of n = 17, see the purple band in Fig. 3. Most importantly this trend aligns perfectly well with observations on the heat capacity curves: those clusters that show well-defined peaks – suggesting first-order-like transitions – are the ones where local minima do not have a significant contribution, while clusters having a potential energy landscape with multiple basins contributing substantially to the partition function have a continuous-like transition, thus only a shoulder on the heat capacity function. This means that we can directly connect the thermodynamic behaviour with the inherent structure of the potential energy landscape.

In the studied range of clusters, sizes n = 7, 13, 15, 17 and 20 are the cases where the global minimum basin is not the most populated basin below the melting transition. The (mW)₂₀ cluster is the most extreme case, where the fraction of the partition function corresponding to the global minimum basin is almost zero at 200 K, indicating that at that point the highly symmetrical cage-like structure has a very low proportion among the configurations, thus the funnel leading to the ground state structure is extremely narrow, making the sampling more challenging. As a consequence, we needed a higher number of walkers to converge n = 20 sufficiently during the nested sampling calculations. Although the location of the low temperature heat capacity peak observed in the case of n = 20 can be aligned well with the temperature where the proportion of the global basin contribution is around 0.5, this is not obvious in the case of n = 7. Moreover, no additional solid–solid transition heat capacity was observed in the case of n = 13, 15 or 17, however it could have been expected according to the partition function ratios.

Fig. 4 shows the energy per molecule for the different global and local minima identified during the sampling and coloured by the average number of neighbours (molecules were treated as nearest neighbours if closer than 3.3 Å). Clusters of even number of molecules have lower energy per molecule, up to n < 18, a trend which aligns well with our observations of first-order-like transitions on the heat-capacity curves. In the case of these global minimum configurations, all the molecules have exactly three neighbours. The trend changes when configurations with a solvated molecule start to appear, and above n = 20 clusters with an odd number of molecules become lower in energy. In the studied range the largest number of neighbours has been seen in the second lowest energy local minima of (mW)₁₇, (structure corresponding to the green band in Fig. 3), with 3.88 average neighbours and the average angle being 93.89.

Fig. 4 Energy per molecule (mW) as a function of cluster size, for the structures identified (and shown in Fig. 3). Global minima are connected by a black line and points are coloured according to the average number of nearest neighbours within the cluster.

It could also be noted that the global minima of clusters n = 17, 19, 21 and 22 are chiral structures and both enantiomers have been found among the sampled configurations.

3.2 TIP3P clusters

Heat capacity curves calculated for clusters modelled by the TIP3P potential are shown in Fig. 5. The most prominent peaks correspond to the evaporation of the cluster, and similarly to our results with the mW model, some cluster sizes exhibit a clear second peak suggesting a first-order-like melting transition (n = 8, 9, 10, 11, 14, 15), while others have only a shoulder or no sign of transition at all. Although the water octamer shows a pronounced melting peak for both the mW and the TIP3P model, it is not surprising that this similarity does not hold for larger clusters, and the same sizes show different behaviour and trends. This is in line with the work of Yin and Landau,²⁵ who found that for certain cluster sizes the characteristics of the heat capacity curves can strongly depend on the potential model employed. A comparison of the nested sampling results to heat capacity curves calculated by parallel tempering and Wang–Landau sampling is shown in Fig. 6, showing an overall good agreement. Similarly to mW, n = 7 shows a small peak at very low temperature, suggesting a solid–solid transition in the structure.

Fig. 5 Heat capacity of (TIP3P)_n clusters as a function of temperature for 5 ≤ n ≤ 16. For every cluster size the results of at least two independent runs are shown. The curves are vertically shifted for ease of readability.

Fig. 6 Heat capacity of TIP3P clusters of 8 and 12 molecules, calculated by nested sampling, Wang–Landau sampling (results from ref. 25) and parallel tempering simulations (results from ref. 23).

In order to describe the structures, we used the Q₄, Q₆ and the average ring size as order parameters to characterise the relative position of the oxygen atoms. For clusters n = 14, 15, 16 we observed different hydrogen bonding networks of the same oxygen backbone in some cases, and if so, the relative energies of these were less than 0.001 eV. We also treated these as belonging to the same basin when calculating the partition function ratios. These results are shown in Fig. 7 along with the structure and energies of the global and local minima found.

Fig. 7 Partition function ratio of different basins (Q_basin/Q_total) as a function of temperature for water clusters with the TIP3P potential. Dark blue corresponds to the global minimum structure; orange, green, light blue and purple are different local minima, while grey corresponds to basins of smaller volume or structures difficult to assign to well defined basins. Snapshots of the structures and corresponding energies are also shown.

The global minimum structures we identified are all identical to the ones published by Kabrede et al.,⁷ except for n = 11, 13. For n = 11, we found a highly symmetrical structure lower in energy, while the previously known global minimum is the structure we identified as the single other explored local minimum, its hydrogen bonding network identical to the one reported in ref. 7 and 9 (and not the one found in ref. 8). For n = 13, the lowest energy structure identified by nested sampling is novel as well. The previously known global minimum (as in ref. 7) was found to be the second lowest energy minimum, and the structure reported in ref. 8 was the third lowest energy minimum in our calculations, although the phase space volume of both structures were found to be very small. We have to emphasize that in the current work the flexible version of the TIP3P model was used with a long-range Coulombic solver, thus we can expect some variations in the energy. We found that the interaction cutoff distance has far the strongest effect in this case: with shortening of the Lennard-Jones cutoff considerably (from 12 Å to 7.5 Å), the order of minima in energy changes, and one of the local minima becomes the lowest in energy for n = 11 and 13. However, the fact that the rest of the cluster sizes agree with previous reports, suggests that the global minimum is robust with respect to variations of the TIP3P model, but sizes n = 11 and 13, which indeed attracted more attention due to inconsistencies of their global minima before,^8,9,43 are in contrast much more sensitive.

From Fig. 7 the transition seen in n = 7 is obvious, the local minima become dominant above 20 K, consistent with the small peak observed on the heat capacity curve. For larger clusters, both the number of local minima and the ratio of their contribution to the partition function are generally small. The most dominant heat capacity peaks were observed for n = 8, 9, 10, and although in these cases the transition from global minima to other structures is relatively sharp, the temperature dependence of the partition function ratios is not significantly different from that of n = 11 and 12, where the heat capacity shows only a shoulder. Overall, the connection between the shape of the melting peak and the partition function ratio picture is not as clear as in the case of the mW model. As the size of the cluster increases the potential energy landscape becomes more complex, and we found converging n = 15 and 16 significantly more challenging. Apart from n = 7, these were the only TIP3P clusters where a dominant local minimum was identified, thus, below the melting transition, the global minimum structure is not the most populated. However, no corresponding peak on the heat capacity curve was identified.

4 Conclusions

We have sampled the potential energy landscape of small water clusters described by the mW (5 ≤ n ≤ 25) and TIP3P (5 ≤ n ≤ 16) potential models using the nested sampling method. Using different order parameters we identified the thermodynamically most relevant minimum basins, and calculated their relative contribution to the total partition function. This also allowed us to connect the thermodynamic properties to the features of the energy landscape: in systems where local minima gave an overall small contribution to the total partition function, a sharp ‘melting’ peak was observed on the heat capacity curve, suggesting a first-order-like transition. This was true for clusters with global minimum of particularly low energy. For clusters, where the global minimum basin is not the most populated at finite temperature and at least two structures compete below the ‘melting’ temperature, the heat capacity shows only a shoulder, consistent with observations of a more continuous-like transition. These findings were particularly clear for the mW model, however the behaviour of the clusters described by the TIP3P potential was more complex.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

Via the authors’ membership of the UK's HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk). JD acknowledges support of an EPSRC summer studentship at the University of Reading. LBP acknowledges support from the Royal Society through a Dorothy Hodgkin Research Fellowship. The authors also thank Albert P. Bartók for useful discussions.

References

1. P. G. Sennikov, S. K. Ignatov and O. Schrems, Complexes and clusters of water relevant to atmospheric chemistry: H₂O complexes with oxidants, Chem. Phys. Chem., 2005, 6, 392–412.
2. P. L. Privalov and C. Crane-Robinson, Role of water in the formation of macromolecular structures, Eur. Biophys. J., 2017, 46, 203–224.
3. M. F. Chaplin, A proposal for the structuring of water, Biophys. Chem., 1999, 83, 211.
4. D. J. Wales and I. Ohmine, Structure, dynamics, and thermodynamics of model (H₂O)₈ and (H₂O)₂₀ clusters, J. Chem. Phys., 1993, 98, 7245.
5. C. J. Tsai and K. D. Jordan, Theoretical study of small water clusters: Low-energy fused cubic structures for (H₂O)_n, n = 8, 12, 16, and 20, J. Phys. Chem., 1993, 97, 5208–5210.
6. J. D. Farrell and D. J. Wales, Clusters of coarse-grained water molecules, J. Phys. Chem. A, 2014, 118, 7338–7348.
7. H. Kabrede and R. Hentschke, Global minima of water clusters (H₂O)_n, n < 25, described by three empirical potentials, J. Phys. Chem. B, 2003, 107, 3914.
8. J. A. Niesse and H. R. Mayne, Global optimization of atomic and molecular clusters using the space-fixed modified genetic algorithm method, J. Comput. Chem., 1997, 18, 1233.
9. D. J. Wales and M. P. Hodges, Global minima of water clusters (H₂O)_n, n ≤ 21, described by an empirical potential, Chem. Phys. Lett., 1998, 286, 65–72.
10. S. Kazachenko and A. J. Thakkar, Improved minima-hopping. TIP4P water clusters, (H₂O)(n) with n ≤ 37, Chem. Phys. Lett., 2009, 476, 120–124.
11. T. James and D. J. Wales, Global minima of water clusters (H₂O)_n, n ≤ 21, described by a five-site empirical potential, Chem. Phys. Lett., 2005, 415, 302–307.
12. S. Maheshwary, N. Patel, N. Sathyamurthy, A. D. Kulkarni and S. R. Gadre, Structure and stability of water clusters (H₂O)_n, n = 8–20: An ab initio investigation, J. Phys. Chem. A, 2001, 105, 10525–10537.
13. K. Liu, M. G. Brown, C. Carter, R. J. Saykally, J. K. Gregory and D. C. Clary, Characterization of a cage form the water hexamer, Nature, 1996, 381, 501–503.
14. J. M. B. Y. Wang, V. Babin and F. Paesani, The water hexamer: Cage, prism, or both. full dimensional quantum simulations say both, J. Am. Chem. Soc., 2012, 134, 11116–11119.
15. R. J. Saykally and D. J. Wales, Pinning down the water hexamer, Science, 2012, 336, 814–815.
16. U. Buck, C. C. Pradzynski, T. Zeuch, J. M. Dieterich and B. Hartke, A size resolved investigation of large water clusters, Phys. Chem. Chem. Phys., 2014, 16, 6859–6871.
17. W. T. S. Cole, O. Yönder, A. A. Sheikh, R. S. Fellers, M. R. Viant, R. J. Saykally, J. D. Farrell and D. J. Wales, Terahertz vrt spectroscopy of the water hexamer-h12 cage: Dramatic libration-induced enhancement of hydrogen bond tunneling dynamics, J. Phys. Chem. A, 2018, 122, 7421–7426.
18. W. T. S. Cole, J. D. Farrell, A. A. Sheikh, O. Yönder, R. S. Fellers, M. R. Viant, D. J. Wales and R. J. Saykally, Terahertz vrt spectroscopy of the water hexamer-d12 prism: Dramatic enhancement of bifurcation tunneling upon librational excitation, J. Chem. Phys., 2018, 148, 094301.
19. J. O. Richardson, C. Pérez, S. Lobsiger, A. A. Reid, B. Temelso, G. C. Shields, Z. Kisiel, D. J. Wales, B. H. Pate and S. C. Althorpe, Concerted hydrogen-bond breaking by quantum tunneling in the water hexamer prism, Science, 2016, 351, 1310.
20. A. N. Tharrington and K. D. Jordan, Parallel-tempering monte carlo study of (H₂O)_n = 6–9, J. Phys. Chem. A, 2003, 107, 7380–7389.
21. D. J. Wales, M. A. Miller and T. R. Walsh, Archetypal energy landscapes, Nature, 1998, 394, 758–760.
22. A. Baba, Y. Hirata, S. Saito, I. Ohmine and D. J. Wales, Fluctuation, relaxation and rearrangement dynamics of a model (H₂O)₂₀ cluster: Non-statistical dynamical behavior, J. Chem. Phys., 1997, 106, 3329.
23. J. Gelman-Constantin, M. A. Carignano, I. Szleifer, E. J. Marceca and H. R. Corti, Structural transitions and dipole moment of water clusters, J. Chem. Phys., 2010, 133, 024506.
24. J. Yin and P. D. Landau, Structural properties and thermodynamics of water clusters: A wang-landau study, J. Chem. Phys., 2011, 134, 074501.
25. J. Yin and P. D. Landau, Wang–Landau approach to the simulation of water clusters, Mol. Simul., 2018, 45, 241–248.
26. T. Kaneko, T. Akimoto, K. Yasuoka, A. Mitsutake and X. C. Zeng, Size-dependent phase changes in water clusters, J. Chem. Theory Comput., 2011, 7, 3083–3087.
27. M. Schmidt and H. Haberland, Phase transitions in clusters, C. R. Phys., 2002, 3, 327–340.
28. D. Wales, Energy Landscapes, Cambridge University Press, 2003.
29. V. Molinero and E. B. Moore, Water modeled as an intermediate element between carbon and silicon, J. Phys. Chem. B, 2009, 113, 4008–4016.
30. W. Jorgensen, J. Chandrasekhar, J. D. Madura, R. Impey and M. L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys., 1983, 79, 926.
31. J. Skilling, Bayesian inference and maximum entropy methods in science and engineering, in AIP Conference Proceedings, vol. 735, 2004, p. 395.
32. J. Skilling, Nested sampling for general bayesian computation, Bayesian Anal., 2006, 1, 833–859.
33. L. B. Pártay, A. P. Bartók and G. Csányi, Efficient sampling of atomic configurational spaces, J. Phys. Chem. B, 2010, 114, 10502–10512.
34. S. Martiniani, J. D. Stevenson, D. J. Wales and D. Frenkel, Superposition enhanced nested sampling, Phys. Rev. X, 2014, 4, 031034.
35. L. B. Pártay, A. P. Bartók and G. Csányi, Nested sampling for materials: The case of hard spheres, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 022302.
36. R. J. N. Baldock, L. B. Pártay, A. P. Bartók, M. C. Payne and G. Csányi, Determining pressure–temperature phase diagrams of materials, Phys. Rev. B, 2016, 93, 174108.
37. R. Baldock, N. Bernstein, K. M. Salerno, L. B. Pártay and G. Csányi, Constant pressure nested sampling with atomistic dynamics, Phys. Rev. E, 2017, 96, 043311.
38. L. B. Pártay, On the performance of interatomic potential models of iron: Comparison of the phase diagrams, Comput. Mater. Sci., 2018, 149, 153–157.
39. N. Bernstein, R. J. N. Baldock, L. B. Pártay, J. R. Kermode, T. D. Daff, A. P. Bartók and G. Csányi, pymatnest, https://github.com/libAtoms/pymatnest, 2016.
40. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys., 1995, 117(1), 1–19.
41. F. H. Stillinger and T. A. Weber, Computer simulation of local order in condensed phases of silicon, Phys. Rev. B: Condens. Matter Mater. Phys., 1985, 31, 5262–5721.
42. P. J. Steinhardt, D. R. Nelson and M. Ronchetti, Bond-orientational order in liquids and glasses, Phys. Rev. B: Condens. Matter Mater. Phys., 1983, 28, 784.
43. F. F. Guimarães, J. C. Belchior, R. L. Johnston and C. Roberts, Global optimization analysis of water clusters (H₂O)_n (11 ≤ n ≤ 13) through a genetic evolutionary approach, J. Chem. Phys., 2002, 116, 8327.

---

Footnote

† Electronic supplementary information (ESI) available: All the minimised structures shown in Fig. 3 and 7 are available in the xyz format. See DOI: 10.1039/c9cp00474b

---