Water bilayers on ZnO(101̄0) surfaces: data-driven structural search

Abstract

Data science methods hold enormous promise for enhancing the efficiency of multiple aspects of theoretical materials science. Here we demonstrate an approach for the use of data science methods for a structural search for high-stability atomic structures in ab initio simulations, via the analysis of a large set of candidate structures. By partitioning a large data set of structures over an appropriate set of variables, we are able to identify a small fraction of structural space into which all low-energy structures are concentrated. Structural search methods may then be applied in the identified area. The method is demonstrated on the problem of finding stable structures of the water bilayer on the ZnO(101̄0) surface.

---

Identification of thermodynamically stable atomic structures of a chemical system under specific conditions is one of the most important problems in theoretical materials science. Although substantial progress has been made in recent years on structural search algorithms, particularly for crystalline solids, via methods such as Ab Initio Random Structure Search¹ (AIRSS), particle swarm optimisation² and evolutionary algorithms,³ the problem remains in general very difficult for systems with more than a small number of atoms per unit cell, since the number of local minima expands exponentially with the number of atoms in the system.⁴ The difficulties are particularly pronounced for the problem of finding optimal structures of surface reconstructions and terminations, since the space of possible structures at surfaces, compared to the space of crystalline structures, tends to contain a large number of local minima which are fairly similar in energy. The extension of existing structural search methods to broader classes of systems is highly desirable for many applications in chemistry and materials science.

Finding the global optimum structure of an atomistic system is equivalent to finding the global minimum on the potential energy surface (PES) that describes the energy of the system as a function of atomic coordinates. Searching for the global minimum of a general arbitrary PES is an NP-hard problem for which no general shortcuts exist.⁵ However, since we are interested exclusively in PESes that represent atomistic configurations, we may take advantage of certain properties of chemical PESes to make the search easier. On an atomistic PES, the maximum curvature is limited by chemical considerations. As a particular consequence, it is known that low-energy minima tend to occur in the vicinity of other low-energy minima. It is also known that the basins which terminate at lower-energy minima tend to have a larger hyper-area than those with higher minima⁶ and that the hyper-areas of the basins of attraction appear to follow a power law distribution,⁷ meaning that low-energy minima are both few in number and surrounded by a relatively large hypervolume. The global minimum is thus almost certain to be concentrated within an area of the PES with many low minima.

Methods such as AIRSS and genetic algorithms, combined with the rapidly growing availability of massively parallel but low-communication bandwidth computing resources, allow the generation of large numbers of local minima on the potential energy surface. With the growing sophistication of data science methods to allow the extraction of insight from large datasets, it is reasonable to wonder if the use of such techniques on datasets of possible structures can be used to accelerate the process of structural search by identifying, in an appropriate basis set, targeted regions of the PES which contain low-lying minima and identifying structural characteristic of this region. It is not necessarily obvious what form such an algorithm should take, nor which variables in which the PES should be expressed in order to best localise the interesting regions of the PES in as few regions as possible. Statistical analysis of the potential energy surfaces of chemical systems is likely to produce new insights that can be used for understanding the general properties of such PESes, allowing efficient data-driven structural search algorithms to be developed.

In this work, we explore the use of data science methods in conjunction with both classical and ab initio simulation to extract insight from atomistic PESes and accelerate the discovery of ab initio optimal structures, by studying as an initial example the potential energy surface of the water bilayer on the ZnO(101̄0) surface, a technologically important system offering a reasonably high degree of complexity. We take advantage of the existence of a good classical potential for the system of interest, which allows a large number of candidate structures to be generated and analysed quickly. We identify structural features of the PES of this system by a partitioning of the space of local minima on the surface over an appropriate set of variables. The analysis of the classical PES can then be used as a map of the ab initio PES, allowing the identification of small portions of the PES (in an appropriate basis set) in which all low-energy structures are concentrated, allowing a structural search concentrated in the regions of interest.

1 ZnO–water background

The ZnO/water interface is of great importance in a number of important processes including photocatalytic water splitting.⁸ Adsorbed surface water also plays a significant role in controlling the morphology⁹ and electronic properties^10,11 of ZnO nanostructures. Due to the complexity of water's hydrogen bonding networks¹² the atomic structure of water at surfaces is often complex and difficult to characterise, however an understanding of these structures is essential for increased control over many important reaction at ZnO surfaces.

The (101̄0) surface of ZnO is the most stable¹³ and best characterised surface facet. Water on this surface is known to form as stable monolayers, bilayers, and higher coverages depending on conditions, however only the monolayer's structure has been well characterised so far.^14–19 Water coverage on ZnO(101̄0) of one monolayer has attracted attention with a variety of experimental^14,19,20 and ab initio theoretical methods.^15–19 Isolated water molecules on ZnO adsorb non-dissociatively,¹⁶ while the monolayer has been shown to form a “half-dissociated” phase in which undissociated H₂O and singly dissociated OH + H motifs exist on alternate Zn–O surface dimers. The monolayer is known to dynamically change between undissociated and fully dissociated structures at room temperature.¹⁹ Quantum chemical methods have also been used to study the interaction of water with small ZnO clusters.^15,21

2 Generation of the structural data set

The first step in the data-driven structural search process is the construction of a large dataset of possible structures representing local minima on the (classical) PES for the water bilayer system. This is undertaken via an unbiased random structure search using the classical ReaxFF potential.

We generated a structural data set consisting of 9073 unoptimised random starting structures, which we label Dataset A. Each surface geometry studied consists of a ZnO slab with boundary conditions periodic in the x and y directions combined with four randomly positioned O and eight randomly positioned H atoms. We set up a wurtzite ZnO slab, oriented in the (101̄0) direction, consisting of sixteen Zn and sixteen O atoms, with two surface dimers along the dimer row (112̄0) direction and four bilayers (approx 10 Å) thick. The size of the periodic supercell in the x and y directions was 5.30 by 6.57 Angstroms, which was chosen to be optimal in the DFT rather than the classical calculations. The four oxygen and eight hydrogen atoms were added to the system at random x and y coordinates and at a z coordinate randomly chosen from a uniform distribution between 0.9 Å and 3.4 Å above the top layer. The cells used had periodic boundary conditions in all three dimensions, with a vacuum spacing of 14.7 Å between the bare slab and its periodic image, ensuring at least 10 Å of vacuum after termination; this is of no importance for the ReaxFF calculations but is important in the later ab initio calculations. The unit cell is shown in Fig. 1.

Fig. 1 Side and end views of ZnO slab used in all calculations. The green shaded region depicts the volume in which the four O and eight H atoms are randomly placed.

3 Variables for the description of the potential energy surface

The distribution of energies of the random starting structures in Dataset A is shown in Fig. 1. Unsurprisingly, most random configurations represent structures which are deeply unphysical resulting in a distribution which is both broad and high in potential energy. The mean energy was −97.25 eV (nearly 40 eV higher than the mean energy of the optimised structures) and the standard deviation is 13.21 eV (Fig. 2).

Fig. 2 Distribution of classical ReaxFF energies of the 9073 randomly generated structures in the data set upon generation (Dataset A), and following optimisation with the HFTN method and ReaxFF potential (Dataset B).

Each starting structure in Dataset A was subjected to a geometry optimisation using the Hessian-Free Truncated Newton method and the ReaxFF potential as implemented in the LAMMPS code²² to generate a new data set of points representing local minima on the ReaxFF energy surface, which we call Dataset B†. The distribution of the energies of the structures in Dataset B is shown in Fig. 1. The minimum energy was −140.5 eV with a mean of −135.5 eV and a standard deviation of 2.4 eV. The distribution of energies in Dataset B exhibits multiple peaks, in contrast to the single-peaked distribution of energies in Dataset A.

4 Variables for the description of the potential energy surface

The 9073 local minima in Dataset B collectively contain an enormous amount of information about the shape of the classical potential energy surface. The challenge is to characterise the relationship between the potential energy and the chemical structure in a way which offers predictive power. In order to do this, we need a convenient set of variables, or descriptors, in which to express each structure which are capable of characterising a given structure in a way that offers predictive power.

Several criteria for such a set of variables are obvious: geometrically similar structures should have similar variables, each unique structure should be characterised by a unique set of variables, and that these variables should not vary when obvious structure-preserving symmetry operations are applied. The PES should also be reasonably smooth and have as few minima as possible when expressed in the new basis. One possible set of variables for describing structures and which exhibits many of these properties which has received some recent attention is the bispectrum representation of the atomic neighbourhoods, as used in the GAP²⁵ and SNAP²⁶ general interatomic potential methods. The bispectrum, however presents a high degree of implementation and interpretational complexity.

We propose instead a set of variables which combines predictive power with computational simplicity. We describe each structure in terms of the sorted shortest distances between atoms of each pair of species. We label the nth-shortest distance between two atoms of species X and Y (taking periodic boundary conditions into account) as L(X,Y)_n−1, using a zero-origin numbering convention. We will refer to this basis set as a basis of Sorted Neighbour Distances, or SNDs. In this particular implementation we ignore all lattice atoms except for those in the top layer, since distances between bulk-like atoms do not differ meaningfully between structures, and compute distances only for the four dimer and twelve adsorbate atoms. Since the surface-layer and adsorbate-layer oxygen species play different roles in the structure, we label them as separate species, O_s and O_a respectively. This gives a total of 120 SND values which characterise each structure. The 120 variables contain some redundancy of information, which is desirable, since we are seeking the most efficient way to subdivide the data.

By representing the set of structures in the SND basis set we may immediately extract insight into the chemistry of the system. Fig. 3 shows the ReaxFF energy of the 9073 structures in the dataset plotted with respect to four selected SND descriptors. Plotting energy against L(H,H)₂ shows a clear distinction between several clumps of clusters – one with a small value of L(H,H)₂ and a high classical energy representing structures with H₂ molecules, and several other clumps with larger L(H,H)₂ values; the lowest-energy structures invariably have a value around 1.7 Å. The plots of L(H,O_a)₄, L(O_a,O_a)₀ and L(O_a,H)₁ likewise show clustering of structures into groups with clear relationships between SND values and energy (Fig. 4).

Fig. 3 Classical energies of the structures in Dataset B plotted with respect to four descriptors in the SND basis set: L(H,H)₂, L(H,O_a)₄, L(O_a,O_a)₀, and L(O_s,H)₁.

Fig. 4 (a) First few splits of the regression tree as applied to the 9073 classical minimum structures in Dataset B to classify their potential energies. (b) Structure of the full regression tree. At each node, the split is to the left if the condition is true.

Of the 9073 local optima in Dataset B, it is sensible to wonder how many are duplicates. Owing to the large number of possible structures it is not possible to distinguish duplicates by the similarity of their energies, however comparison of the SND values can identify duplicates. If we define two structures to be identical if the supremum norm of the difference between the SNDs of the two structures is less than 0.02 Å, we find just 88 structures, exclusively concentrated in the lowest 1 eV, are identical to some other structure, including one structure found five times, two structures found four times, and many structures found two or three times. The higher duplication rate among lower-lying minima is consistent with the well-known tendency of low-lying structures to lie in larger basins.¹ Our 9073 starting structures have thus reached in excess of 9000 individual local minima. The low level of duplication suggests that not every possible minimum of the PES has been reached, however the level of duplication among the lowest-energy structures suggests that a reasonable level of coverage has been achieved among these structures.

5 Partitioning of the potential energy surface

The next task is to divide the potential energy surface into regions which are characterised by similar structures. Many clustering algorithms exist which classify structures into groups according to their distance apart in variable space, however we wish instead to classify structures according to their similarity in the variables which best determine the structure's potential energy, while giving low importance to variables which have no significant impact upon the energy. Regression trees²⁷ are a class of algorithms which recursively split a variable space according to the variable which best predicts a given output variable. The regression tree splits the potential energy surface into groups of points whose potential energy is as similar as possible in the most effective possible way, by recursively splitting the data points into groups according to an algorithm which maximises the extent to which the remaining data points in each group are similar in energy. By using a regression tree we partition the space in a way that at each step of the algorithm ensures that the energy values points in each partition (or bin) are as close to constant as possible.

Fig. 3(a) shows the first few nodes in a regression tree fitted to Dataset B as expressed in the SND basis set. The regression tree was generated using the rpart package²⁸ in R.²⁹ At the first node of the tree, it is found that the best possible split in the data is L(H,H)₁ ≥ 1.33 Å, i.e. the third-shortest hydrogen–hydrogen distance in the sample is greater than this threshold. In physical terms, this distinguishes those structures which have two or more H₂ molecules (higher in energy) from those which do not (lower in energy). The next split for structures on the low-energy side is L(H,H)₀ ≥ 1.324 Å. The third split for those in the more stable set is on the criterion L(O_a,O_a)₀ > 2.39 Å. The tree is constructed using a complexity cost parameter²⁸ of 10⁻⁸. Fig. 3(b) shows the full regression tree, consisting of 794 classification bins containing 11.4 structures on average. Each structure is judged by up to 28 criteria.

Fig. 5(a) shows the distribution of classical energies for the structures in Dataset B as a function of partition bin. Bin indices run from 1 to 1587 with some omitted due to the numbering of non-terminal nodes in the tree. The bin identity predicts the energy to within 0.29 eV (one standard deviation) on average. For lower-energy bins the energies are better constrained, with those bins whose average energy lies below −138 eV exhibiting a standard deviation in energy of 0.17 eV.

Fig. 5 (a) Distribution of classical energies of 9073 ReaxFF-optimised structures into the classification bins. (b) Distribution of DFT energies of 2382 structures, sampled as the three lowest-energy members of each of the classically determined bins, following optimisation using DFT methods. Most of the minima found using classical methods may be seen to have optimised into a relatively small number of structural bins.

6 Characterisation of the ab initio PES

So far we have used random structural search combined with statistical methods to gain insight into the shape of the classical PES and the distribution of minima thereon. The goal of this work is to utilise the insight derived from the classical PES to gain insight into the ab initio PES. In attempting to transfer the classical PES binning to the ab initio PES we make an assumption of similarity between the two PESes – that minima on the ab initio PES are located close to minima on the classical PES. The relative depths of structures on the classical and ab initio PES cannot, however, be assumed to match closely.

From each of the 794 classification bins we chose the three lowest-lying structures as bin representatives, and conducted a DFT geometry optimisation on each to the nearest minimum on the ab initio PES. For each DFT geometry optimisation, we used the VASP code with pseudopotentials of the PAW type³⁰ and the exchange-correlation functional of Perdew, Burke and Ernzerhof.³¹ Wavefunctions were expanded in a basis set with a cutoff of 300 eV and only the gamma point was used; these parameters are somewhat coarse but are suitable for screening the large data set. Geometry optimisation used a conjugate gradient algorithm, and proceeded until the change in total energy at two subsequent steps was less than 10⁻⁷ eV. The set of 2382 DFT-optimised bin representative geometries is labelled as Dataset C.

Fig. 6 shows the DFT energies for the classically-optimal and DFT-optimal structures as a function of the classical energy of the starting structure. For the starting geometries (blue) the classical and DFT energies are similar and close to normally distributed, with a standard deviation of 1.24 eV. Notably we see a clustering of post-optimisation DFT energies into five distinct clumps.

Fig. 6 Blue: DFT energies for the 2832 bin representative structures picked from Dataset B, prior to DFT optimisation, compared to the classical energy of the same structure. Green: DFT energies following optimisation of each of the structures, vs. the classical energy of the classical structure from which it was generated.

The post-DFT-optimisation structures in Dataset C may be sorted into bins using the same decision tree used for the classical structures. Fig. 5(b) shows the U_DFT values and bin locations as a function for all structures in Dataset C. Notably, many structures from Dataset B have moved to a different bin during the optimisation, and 243 of 794 bins are occupied. Also notable in Fig. 5(b) is the exclusive segregation of all structures into a small number of bins in which all energies lie below −200 eV, and the remainder of bins in which all structures have energies above −200 eV. These two sets of structures are separated by a wide gulf in excess of −200 eV.

There are a total of 396 optimised structures whose final energy is less than −200 eV. These come from 187 origin bins, but have final structures in only 27 destination bins. It thus appears that only 27 out of the original 794 bins contain structures which are likely to be competitive. This suggests that concentration of the structural search into these particular regions of the PES is likely to produce energetically competitive structures in an efficient way.

Fig. 7 depicts a decision tree by which we may determine whether a given structure exists in one of the 27 low-energy destination bins.

Fig. 7 Decision tree for determining whether a structure lies within one of the 27 bins of interest. At each node, we split to the left if the condition is false and to the right if it is true.

7 Targeted optimisation of structures in region of interest

Having determined the regions of the PES in which all low-lying DFT minima exist, we may now target our energy search at structures within this region. Dataset B contains 328 structures which lie in these bins. It would be possible with relatively little computational time to generate a larger pool of structures, however for this initial proof of concept these 328 structures will suffice. Each was subjected to a DFT geometry optimisation using the same parameters used in the previous Section. The set of stable structures was labelled Dataset D.

As anticipated, structures optimised from the targeted bins exhibit a high degree of stability, with a mean energy of −201.7 eV and a standard deviation of 0.24 eV. All but two of the optimised structures have an energy below −201 eV, and the most stable structure has an energy of −202.1 eV. Fig. 8 plots the final DFT energy against the starting classical energy for each of the structures in Dataset D.

Fig. 8 Post-DFT-optimisation DFT energy plotted against pre-DFT-optimisation classical energy for the 328 targeted-region structures in Dataset D.

We find several geometrically distinct structures with extremely similar energies to be the most stable structures in the dataset. Fig. 9 shows the geometries of the three lowest-energy structures. The second- and third-lowest energy structures have DFT energies which are just 0.009 eV and 0.01 eV higher than those of the first. Four additional structures were found to be within 0.025 eV (k_BT at room temperature) of the lowest-energy structure. Distinguishing which of these structures as the true ground state is well beyond the capabilities of density functional theory. Moreover, the question is not of importance, since the energetic differences between structures are less than k_BT at room temperature, implying that under realistic conditions the water bilayer will exist in a mix of structures.

Fig. 9 The (a) the lowest, (b) second-lowest and (c) third-lowest DFT-energy structures in Dataset D, combined with (d) the lowest-classical-energy structure from Dataset B. Hydrogen, oxygen and zinc atoms are shown in white, red and grey respectively. Hydrogen–oxygen distances less than 1.2 Å (corresponding to covalent bonds) are shown in green while those between 1.2 Å and 2.4 Å (corresponding to hydrogen bonds) are shown in blue.

Several common characteristics of these stable structures may be identified, however; in particular, the lower-lying monolayer of all structures have the same “half-dissociated” character as is present in the water monolayer, with one out of every two surface Zn atoms being bonded to a hydrogen atom. It is notable that the three lowest-energy structures show different structural characteristics to the lowest classical energy structure in Dataset B (Fig. 9(d)) in which the water molecules lie substantially closer to the surface.

Due to near-degeneracy of the distinct lowest-energy structures in the data set, we do not concern ourselves with identifying the true ground state structure of the water bilayer, and satisfy ourselves instead that we have a procedure for reliably generating geometrically distinct structures out of the set of low-energy water bilayer structures. Further work could generate a larger set of these structures and study them with more accurate methods.

8 Discussion

These results function as a proof of concept for the use of data science methods upon large datasets of candidate structures in an appropriate basis set in order to extract insight about the potential energy surface and identify classes of structures which are close to the global minimum. These methods hold promise to enhance structural search by providing a scheme by which to identify regions of interest on the PES which can be searched more thoroughly. The procedure used in this exploratory work is by no means optimal, and many avenues exist for improving the efficiency of this process for this or other systems.

The SND basis set has been shown to be hold a great deal of predictive power for the water bilayer system. It is to be anticipated that other surface and interfacial systems such as crystal surface reconstructions, self-assembled monolayers and small clusters may benefit from expression and partitioning in this basis set. The utility of the SND basis for the description of solid crystalline, defective or amorphous systems is a possible object of future study.

The use of the regression tree algorithm represents a relatively simple yet powerful way of partitioning the set of structures. We have not explored other possibilities which may involve making splits on multiple variables simultaneously. The combination of regression trees with clustering algorithms is another potential area of interest. The incorporation of principal component analysis into the description of the structural data³² may also prove useful by providing an alternative basis set in which the partitioning may be more efficient.

The implementation in this case requires the existence of an efficient empirical potential. With the growing sophistication of general-form empirical potentials previous work on high-accuracy low-transferability empirical potentials to be used for ab initio thermodynamic calculations provides additional motivation for improved methods for on-the-fly fitting of empirical potentials,³³ which could be achieved with new methods such as SNAP.²⁶

For illustrative purposes we have in this case computed potential energies of a much larger set of bin-representative candidate structures than would be required in practice. Once the degree of mismatch between the classical and ab initio PESes has been established it is not necessary to compute representative energies of bins we may be confident will not be competitive. For an exhaustive structural search, however, it would be desirable to continue generating and DFT-optimising candidate structures within the optimal structural bins until all low-lying structures have been found multiple times.

Data-driven structural search is best viewed as a complement to existing methods rather than a replacement. In this case, for instance, we have used data-driven searching to enhance ab initio random structure search. We may make some remarks about the utility, efficiency and range of applicability of the present method in comparison to the straightforward implementations of existing algorithms including ab initio random structure search and genetic algorithms.

The method presented in this paper involves the production of a large dataset of local minima on a classical potential energy surface, the classification of the data set into structural bins and the ab initio optimisation of bin representative structures. The time taken by the analysis is insignificant, but the time taken to generate the classical minima and evaluate the ab initio energies are significant. Compared to ab initio random structure search, this method will be worthwhile if the time spent on these steps is less than the time which would otherwise be spent generating and optimising points in the non-preferred region, which in the present example is around 99% of structures. The present implementation is ideal for structural search on chemical PESes which contain an extremely large number of minima of roughly similar energy. For this reason, systems such as surface terminations and reconstructions are likely to be most suited for this method.

An important question is the relationship between the size of the random datasets and the probability of finding the correct ground state structure. As with other structural search methods, however, this question defies an easy answer. In this particular case, the probability depends on details of the shapes of both the classical and ab initio PESes and the relationship between the two which are not easily accessible. Future work may assist in developing heuristics for determining the appropriate choice of these parameters.

All existing structural search methods produce a large amount of data about the structure of the potential energy surface as they run, which in the future could be put to better use. In ab initio random structure search, the locations and energies of the minima evaluated so far are not used to inform the generation of future structures, resulting in a vast amount of data left unused. Particle swarm optimisation methods take the global structure of the PES into account in a way, by letting points during a single optimisation communicate information with each other. Genetic algorithms take the global structure of the PES into account in a different way, by cross-breeding existing structures in ways likely to produce new low-energy structures. Data-driven partitioning of the PES in an appropriate basis set to identify regions of greater interest may be able to enhance implementations of each of these methods.

9 Conclusions

We have demonstrated a method for structural search which makes use of data science methods for partitioning the potential energy surface to identify regions which contain low-lying minima. We have demonstrated the value of the SND variables basis set in expressing the geometry of a structure in a way that is computationally and conceptually simple and has predictive power. Using these methods we have identified a number of structures for the water bilayer on the ZnO(101̄0) surface with near-identical energies. All structures maintain the alternation between one terminated and one unterminated Zn surface dimer atom seen in the water monolayer structure.

The use of data science methods for structural search provides a promising new avenue of research which may allow greatly enhanced efficiency and automation of structural search in the future. Further research may reveal better descriptors and better algorithms. Ultimately, ab initio simulation gives us an enormous amount of data about the shape of potential energy surfaces which has not previously been taken advantage of, and by mining this data we may be able to produce predictors which understand chemistry better than we do.

References

1. C. J. Pickard and R. J. Needs, J. Phys.: Condens. Matter, 2011, 23, 053201.
2. Y. Wang, J. Lv, L. Zhu and Y. Ma, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 094116.
3. A. R. Oganov and C. W. Glass, J. Chem. Phys., 2006, 124, 244704.
4. F. H. Stillinger, Phys. Rev. E, 1999, 59, 48.
5. D. Wales, Energy Landscapes: Applications to Clusters, Biomolecules and Glasses, Cambridge University Press, 2003.
6. J. P. K. Doye, D. J. Wales and M. A. Miller, J. Chem. Phys., 1998, 109, 8143.
7. C. P. Massen and J. P. K. Doye, Phys. Rev. E, 2007, 75, 037101.
8. X. Yang, A. Wolcott, G. Yang, A. Sobo, R. C. Fitzmorris, F. Qian, J. Z. Zhang and Y. Li, Nano Lett., 2009, 9, 2331.
9. Z. Fu, Z. Wang, B. Yang, Y. Yang, H. Yan and L. Xia, Mater. Lett., 2009, 61, 4832–4835.
10. A. J. Nozik and R. Memming, J. Phys. Chem., 1996, 100, 13061–13078.
11. Y. Li, F. D. Valle, M. Simonnet, I. Yamada and J.-J. Delaunay, Appl. Phys. Lett., 2009, 94, 023110.
12. D. Marx, Science, 2004, 303, 634–636.
13. C. Tang, H. F. Wilson, M. J. S. Spencer and A. S. Barnard, Phys. Chem. Chem. Phys., 2015, 17, 27683–27689.
14. B. Meyer, D. Marx, O. Dulub, U. Diebold, M. Kunat, D. Langenberg and C. Wöll, Angew. Chem., 2004, 43, 6641–6645.
15. J. B. L. Martins, J. Andres, E. Longo and C. A. Taft, Int. J. Quantum Chem., 1996, 57, 861–870.
16. B. Meyer, H. Rabaa and D. Marx, Phys. Chem. Chem. Phys., 2005, 8, 1513–1520.
17. D. J. Cooke, A. Marmier and S. C. Parker, J. Phys. Chem. B, 2006, 110, 7985–7991.
18. A. Calzolari and A. Catellani, J. Phys. Chem. C, 2009, 113, 2896–2902.
19. O. Dulub, B. Meyer and U. Diebold, Phys. Rev. Lett., 2005, 95, 136101.
20. G. Zwicker and K. Jacobi, Surf. Sci., 1983, 131, 179.
21. J. B. L. Martins, C. A. Taft, E. Longo, E. A. S. de Castro, W. F. da Cunha, J. R. S. Politi and R. Gargano, Int. J. Quantum Chem., 2012, 112, 3223–3227.
22. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
23. D. Sheppard, R. Terrell and G. Henkelman, J. Chem. Phys., 2008, 128, 134106.
24. E. Bitzek, P. Koskinen, F. Gaehler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
25. A. P. Bartok, M. C. Payne, R. Kondor and G. Csanyi, Phys. Rev. Lett., 2010, 104, 136403.
26. A. P. Thompson, L. P. Swiler, C. R. Trott, S. M. Foiles and G. J. Tucker, J. Comput. Phys., 2015, 285, 316–330.
27. L. Breiman, J. Friedman, R. Olshen and C. Stone, Classification and Regression Trees, Wadsworth and Brooks, Monterey, CA, 1984.
28. T. Therneau, B. Atkinson and B. Ripley, rpart: Recursive Partitioning and Regression Trees, 2015.
29. R Core Team, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, 2015.
30. P. E. Blochl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
31. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
32. M. Fernandez, M. Trefiak and T. K. Woo, J. Phys. Chem. C, 2013, 117, 14095–14105.
33. H. F. Wilson, Comput. Phys. Commun., 2015, 197, 1–6.

---

Footnote

† The HFTN algorithm implemented in LAMMPS was chosen as the method for structural optimisation after tests showed it to produce a lower-energy distribution of structures than three alternative methods: steepest descent, conjugate gradient and two different damped dynamics methods (quickmin²³ and fire²⁴).

---