Grand and Semigrand Canonical Basin-Hopping

We introduce grand and semigrand canonical global optimization approaches using basin-hopping with an acceptance criterion based on the local contribution of each potential energy minimum to the (semi)grand potential. The method is tested using local harmonic vibrational densities of states for atomic clusters as a function of temperature and chemical potential. The predicted global minima switch from dissociated states to clusters for larger values of the chemical potential and lower temperatures, in agreement with the predictions of a model fitted to heat capacity data for selected clusters. Semigrand canonical optimization allows us to identify particularly stable compositions in multicomponent nanoalloys as a function of increasing temperature, whereas the grand canonical potential can produce a useful survey of favorable structures as a byproduct of the global optimization search.


INTRODUCTION
Structure prediction is essential in many areas of computational science, ranging from molecular physics and biochemistry to soft and condensed matter. For a given system with definite size, the global optimization problem is usually nontrivial owing to highdimensional potential energy landscapes, and many methods have been proposed to locate low-energy configurations. 1−4 The global minimum is fundamentally important and often carries essential insight into the interactions responsible for the emergence of specific morphologies, and it plays an important role in explaining self-assembling motifs and symmetries. 5,6 However, in many applications, temperature can play a significant role and the entropic contribution to the free energy of individual configurations becomes important. Examples of entropy-driven structural transitions have been reported in atomic clusters, 7,8 proteins, 9 colloids, 10 glasses, 11 and pressurized materials. 12 The determination of configurations that are low in free energy can proceed by the a posteriori analysis of molecular simulations, often employing biases in order to sample the energy landscape more efficiently and based on system-dependent order parameters. 13−17 Such free energies are global and can encompass many potential energy minima, as expected in the context of phase transitions. Local free energies can also be defined for individual isomers, which are related to global quantities through suitable grouping procedures that require additional knowledge about the connectivity of the minima. 18,19 It is possible to locate the global free energy minimum among an existing database of structures by evaluating the entropy using the harmonic approximation 20,21 and, when affordable, incorporating anharmonic corrections. 21−23 More recently, calculating the free energy directly on-the-fly during global optimization was proposed, 24 producing a promising procedure for exploring energy landscapes, since the free energy minimum was encountered faster this way than by postprocessing a sample based on optimizing the potential energy alone. In the present contribution, we further extend this approach by addressing systems with variable size or composition, which should be treated in the grand canonical or semigrand canonical ensembles, respectively. Fluctuations in the number of particles occur in the case of nucleation of fluids and their absorption into porous materials, as well as in the increasingly important problem of reversible gas storage for energy production. These ensembles correspond to situations in which the system exchanges particles with a (possibly fictitious) reservoir, thereby controlling size or composition at fixed temperature.
Grand canonical ensembles are characterized by the chemical potential (or chemical potential difference) and a finite temperature, with the Gibbs free energy being the potential of interest that, in turn, controls the size or composition around equilibrium. As the chemical potential varies, changes in the Gibbs free energy are indicative of different regimes in which the system grows or shrinks or reaches equilibrium values in absorption isotherms. In semigrand canonical ensembles, particularly stable compositions should be manifested by plateaux in the segregation isotherms. Compared to free energy global optimization, the need to calculate local Gibbs free energies for systems with varying size or composition requires sampling these additional variables as well. The extra degrees of freedom further justify the use of the harmonic approximation to approximate the entropy component in a computationally efficient manner. In practice, the harmonic approximation requires calculating the vibrational frequencies at local minima, which involves constructing and diagonalizing the dynamical matrix (the mass-weighted Hessian).
Our results for atomic clusters in the grand canonical ensemble indicate that increasingly large clusters are obtained as the chemical potential is increased or the temperature is decreased, in agreement with conventional nucleation theories. The results of grand canonical basin-hopping simulations in the harmonic approximation are also found to agree with a model for the grand canonical partition function fitted to reproduce the size-dependent heat capacity for specific clusters. Our other example application deals with model nanoalloys treated in the semigrand canonical ensemble, which is at fixed total size but varying alloy composition. Here, we show that the semigrand canonical basin-hopping method efficiently locates the stability plateaux in the composition isotherms that were previously reported based on alternative simulation methods. 25 The article is organized as follows. The next section describes the method in its general formulation and details the harmonic expression employed for the local Gibbs free energy associated with individual potential energy minima. Practical details regarding the implementation of the basin-hopping method are also given in relation to Monte Carlo moves that change the system size. The results for the grand canonical ensemble are presented in Section 3 and discussed in the context of the homogeneous nucleation problem. The semigrand canonical application to model nanoalloys is described in Section 4, followed by concluding remarks in Section 5.

METHODS
2.1. Grand Canonical Formulation. In the grand canonical ensemble, the volume , temperature T, and chemical potential μ are fixed, whereas the pressure, energy, and number of particles, N, can fluctuate. The grand partition function describing an equilibrium distribution is then where β = 1/k B T, with k B the Boltzmann constant, and Z N T ( , , ) is the canonical partition function of the system with fixed number of particles, N. In previous work, we have constructed Z N T ( , , ) as a sum over contributions from local minima, i.e., the superposition approach. 20,21,26−29 Here, the classical vibrational density of states can be written as where κ(N) = 3N − 6 is the number of nonzero eigenvalues for the is the geometric mean vibrational frequency of minimum α, ν α N (j) is the normalmode frequency of the jth mode in this minimum, and E α N is the corresponding potential energy. The approximation in eq 2 corresponds to using harmonic vibrational frequencies. The superposition approach can also incorporate quantum effects 23 and anharmonicity, 21−23,28,30−32 but it is usually employed in the harmonic approximation to obtain a rapid survey of thermodynamic properties, which is guaranteed to be ergodic by construction.
In recent work, we have demonstrated how the superposition framework can usefully be applied within grand and semigrand canonical formulations to examine equilibrium thermodynamics. 25 In the present contribution, we show how a potential function based on the grand (and semigrand) canonical ensembles can be used in the context of basin-hopping global optimization. This approach is a natural extension of the free energy basin-hopping method, 24 generalized so that the size (or composition) is permitted to change. To find the largest contribution to the grand potential from local minima of any size N, we adapt the acceptance criterion to use the potential where we have included the rotational partition function π β = | | ℏ q I 8 / A Metropolis acceptance criterion was applied at each step, u s i n g a n a c c e p t a n c e p r o b a b i l i t y b a s e d o n min{1,exp[−(ξ α(new) N(new) − ξ α(old) N(old) )/k B T GCBH ]}, where T GCBH is a fictitious temperature parameter that determines how often uphill moves in ξ are accepted. Along with the temperature parameter, another key choice in basin-hopping that affects efficiency is the coordinate perturbation scheme applied before each local minimization. Here, we employed perhaps the simplest scheme based on Cartesian coordinate displacements, drawn from a uniform distribution with a fixed maximum value. Many other possibilities have been considered in the literature, along with variations in the acceptance condition, and the present approach could be combined with any of these methods in future work. So long as the key local minimization is included, 2,36,37 efficiency gains might be possible. Here, we adopted one other modification, since moves that involve changes in N are likely to be much more disruptive than geometrical perturbations. We therefore considered moves changing the number of atoms only at intervals of Δ basinhopping steps. Before changing the number of atoms, the structure of the current minimum in a Markov chain over blocks was saved, along with the value of ξ α N . After adding or removing an atom, the structure obtained after minimization was used as the initial seed for a local Markov chain of Δ coordinate perturbations, accepting moves according to the condition At the last step of each block, the minimum with the lowest value of ξ α N was used in a block accept/reject test and compared with the current minimum in the block Markov chain. If the block move was rejected, then the current structure was reset to the one saved in the block Markov chain before proposing another move that changed the number of atoms. The same BH temperature

Journal of Chemical Theory and Computation
Article parameter was employed for both types of acceptance check, although different values, or indeed different criteria, could certainly be considered.
Two different schemes were compared for grand canonical basin-hopping (GCBH) in the present work to analyze the role of different contributions to the potential. The first scheme used ξ α N , as defined above, and the second scheme based the sampling on ξ α omitting the last term of the right-hand side of eq 3 involving the canonical partition function of the minimum. In effect, this version corresponds to neglecting entropic contributions. The results are compared in Section 3.
The grand canonical potential in eq 3 can be reinterpreted as a semigrand canonical potential for a binary system with variable numbers of A-and B-type particles, N A and N B , respectively, but N = N A + N B fixed. The potential employed in the accept/reject criterion is then In Section 3, we describe a semigrand canonical basin-hopping (SGCBH) scheme with and without the rotational contribution ( metric perturbations disabled, and particle insertion/deletion moves replaced by exchange moves that transmute particles.

Steps That Change the Cluster Size.
Steps within the blocks of constant size were proposed in a manner similar to that in previous work, including procedures to move weakly bound or surface atoms. 37 Steps to change size were proposed by adding or deleting single atoms, with probability p + and p − = 1 − p + , respectively. All of the results reported below simply used p + = p − = 0.5 throughout. Atomic clusters bound by a pair potential were chosen for the first application of grand canonical basin-hopping. For pair potentials (but also for embedded-atom potentials, as used in our semigrand canonical basin-hopping examples), we can easily identify the most weakly bound atom for each minimum, and this was the atom removed in steps that reduced the cluster size. To add an atom, the center of coordinates was first located, along with the largest atomic radial distance, r max . An atom was then added at a random point on the sphere with radius r max + δr (where δr = σ for Lennard-Jones systems) by generating a three-component vector with each entry drawn from the normal distribution with zero mean and unit variance 38 and then normalizing appropriately. Gaussian random variables were generated using the Box− Muller algorithm. 39 Following each atom addition or deletion, the resulting configuration was immediately minimized; if this quench failed, then the attempted size change was simply rejected. Local minima were also rejected if they did not correspond to connected single clusters. Here, we used a depth first search 40 to check for a percolating network of atoms, which proved to be particularly useful in previous studies of clusters bound by shortrange potentials. 41 The present results therefore exclude fragmented systems, focusing on single clusters.

APPLICATION OF BASIN-HOPPING TO A GRAND CANONICAL POTENTIAL
We considered clusters bound by the Lennard-Jones (LJ) potential, 42 where the potential energy is ϵ is the equilibrium pair well depth, 2 1/6 σ is the equilibrium pair separation, and r ij is the distance between particles i and j.
Depending on the values of T and μ, the GCBH runs either exhibit sizes that shrink to a single atom or grow quite steadily to contain hundreds of atoms for the run lengths considered here. To prevent indefinite growth, an upper limit of 1000 atoms was applied.
3.1. Model for the Size-Dependent Density of States. In order to interpret the observed behavior, a model for the sizedependent partition function was constructed for comparison with the GCBH runs using data available for selected LJ clusters. Here, the input included results for predicted global minimum energies, point groups, normal modes, and moments of inertia up to LJ 1610 . 37 where the subscript 0 refers to the putative global minimum for N atoms, and Stirling's approximation was used for ln N!. The first term is the usual Boltzmann factor, associated with the global potential energy minimum for N. It serves as the lower bound to the potential energy distribution of local minima at each size. The next two terms derive from the harmonic approximation to the vibrational density of states, again for the global minimum, which serves as a baseline. The fourth term is the N-dependent part of the rotational partition function for the global minimum, and the N! term comes from the number of distinct permutation-inversion isomers, with the corresponding point group order in the final term. This last term is a representation of the quadrature over the potential energy

Article
In selecting the above model, we have been guided by the densities of minima calculated in previous work 31 for LJ 31 and LJ 75 . The energy range and growth rate, ΔE(N) and γ(N), were fitted to best reproduce the heat capacity curves from parallel t e m p e r i n g M o n t e C a r l o d a t a , u s i n g ΔE(N) = ΔE 0 + NΔE 1 + N 2/3 ΔE 2 and γ(N) = γ 0 + Nγ 1 + N 2/3 γ 2 . Various other representations were also considered; the choice in eq 7 produces an acceptable fit for the melting peak in the heat capacity curves ( Table 1). The fitting is only qualitative, since we make no attempt to reproduce additional peaks due to phase-like transformations below the melting point. 26,29,37,49−57 We require only a qualitative description of the growth in the number of local minima, convoluted with the change in rotational and vibrational densities of states relative to the global minimum reference values. The results presented below correspond to H = 30 quadrature slices, which is sufficient to converge the properties of interest. The resulting heat capacities for some selected cluster sizes are shown in Figure 1.
It is instructive to consider the individual summation terms in eq 1 associated with a particular N at selected fixed values of the chemical potential μ and temperature T. In Figure 2, we plot  8) with Z N T ( , , ) evaluated from the model defined in eq 7 and μ W N T ( ; , , ) representing a relative (i.e., not normalized) probability of cluster size N for a given μ, , and T. Note that test runs where the rotational contribution was omitted gave essentially the same results (data not shown).
The nonmonotonic trends in Figure 2 exhibit a minimum, which corresponds to the least probable cluster size N lp , i.e., the N that minimizes μ W N T ln ( ; , , ) . (The relationship between N lp and the definition of a critical nucleus will be investigated in more detail in future work.) Figure 3 shows that, for a range of μ and T values, N lp (μ, T) is consistently intermediate between N min = 3 and N max = 1610. In contrast, Figure 4 shows that the value of N that maximizes μ W N T ln ( ; , , ) , denoted by N mp as short-hand for "most probable", always appears either at the maximum (N max ) or the minimum (N min ) size considered. Since N mp accounts for the largest contribution to μ Ξ T ( , , ) and always occurs at one end of the size range considered, one result of a GCBH run is a prediction of cluster growth or dissociation for a given μ and T. If a survey of low-lying minima is not required, then each run can be terminated once the growth or dissociation is clear.
3.2. GCBH Results. We find that cluster growth is favored by low temperature and larger μ, and the values of μ above which GCBH runs produced growth, μ*(T), are summarized as a function of temperature in Table 2. The same values were obtained in runs of 10 7 GCBH steps with size moves attempted

Journal of Chemical Theory and Computation
Article every 50 or 1000 steps. The corresponding block sizes for accept/reject testing after a size move were set to 40 and 900. Runs that exhibited cluster growth did not reach the limit of 1000 atoms. The maximum Cartesian coordinate perturbation was fixed at 0.4, and the temperature parameter employed in both the standard BH steps (constant size) and in the block accept/reject tests was fixed at T GCBH = 1.5 reduced units throughout in this initial survey. The temperature parameter must be carefully distinguished from the physical temperature that appears in the grand potential. As in standard basinhopping, effective values for the temperature parameter should be large enough for the algorithm to escape from low-lying local minima and small enough for low-lying regions to be properly explored. The GMIN code includes a replica exchange option where exchanges are attempted between parallel runs with different values of the basin-hopping temperature parameter, 58 but we have not explored this extension in the present work.
The observation of a minimum in the relative probability size distribution for ranges of T and μ seems to be consistent with the appearance of a critical nucleus in classical nucleation theory. 59−62 In the latter framework, clusters containing fewer particles than the critical value are predicted to shrink and clusters containing more particles should exhibit spontaneous growth. However, we note that various subtleties exist in the treatment of cluster nucleation by simulation, 63,64 and connections with the procedure adopted in the present work require further investigation. Our results simply represent a search to optimize the potential ξ α N , defined in eq 3, with no attempt to describe a vapor/droplet equilibrium. For reference, recent calculations of the chemical potential for saturation (equilibrium), μ s , with alternative potentials for argon produce values around −8.7 and −10.4 for k B T/ϵ ≈ 0.61 and 0.78, respectively. 65 Runs using ξ α N (0) = E α N − Nμ, which do not require normalmode analysis, gave results consistent with the low-temperature limit of ξ α N , with μ* = −5. The lowest minima located as a function of cluster size are an interesting and potentially useful byproduct of the GCBH runs. Since systematic global optimization is not being performed explicitly for each cluster size, we would not expect to locate the true global minima for all N in the range encountered. Nevertheless, searches that involve changes in size, and the block move structure, may provide a useful survey of the low-lying morphologies. To illustrate this possibility, the lowest potential energies encountered in two selected GCBH runs compared to the lowest known minima are plotted as a function of N in Figure 5. These plots include a range of μ values for temperatures of 0.001 and 0.1 and illustrate the trends that also appear for higher temperatures. At T = 0.001, runs with μ values large enough for cluster formation often encounter the lowest known minimum up to N around 200 atoms. For larger sizes, the potential energy difference generally increases with N.
For T = 0.1, the lowest minima encountered are close to the global minimum for smaller sizes, and the energy difference Note that N mp always occurs at (or very close to) one limit of the cluster size range considered.
The μ values considered were in integer steps.

Journal of Chemical Theory and Computation
Article again generally rises with N. However, for μ = −2, we see the difference decrease above about N = 370. These trends basically reflect the proportion of steps that the GCBH run spends around each cluster size. The better a given N value is sampled, the smaller the deviation from the energy of the putative global minimum. The results for higher temperatures (omitted for brevity) are consistent with this pattern. To obtain a survey of low-energy structures as a function of size, we would therefore choose a low temperature, together with a large enough value of μ to produce clustering.

APPLICATION OF BASIN-HOPPING TO A SEMIGRAND CANONICAL POTENTIAL 4.1. Model and Calculation Details.
To demonstrate semigrand canonical basin-hopping (SGCBH) and compare with previous work, 25 we consider icosahedral Ag n Pd 55−n nanoalloys modeled by the same many-body (Gupta) potential (with the same parameter values). Using independent SGCBH runs, each at a different value of Δμ, we seek the composition and mixing pattern that minimize the value of ξ α N Ag defined in eq 5, with the rotational contribution omitted, as in ref 25 Each SGCBH run consisted of 10 5 basin-hopping steps at k B T SGCBH = 0.01 eV, with the stoichiometry reset to N Ag = 28 and species distribution randomized every 500 steps to enhance exploration.
4.2. Equilibrium Composition. The Ag content of the lowest encountered configuration obtained at the three physical temperatures of 0, 100, and 300 K is plotted versus Δμ in Figure  6. At low temperature, the segregation isotherms display particularly striking plateaux for N Ag = 12 and 42, as well as secondary steps at intervening sizes. These plateaux indicate particularly stable compositions and are associated here with highly symmetric motifs in which the nanoalloy is enriched in silver at the surface, either at the vertices (N Ag = 12) or entirely (N Ag = 42). Those motifs were also identified in ref 25, and after scaling the chemical potential difference by a factor of 2, the results agree well with those obtained here.
In ref 25, the harmonic superposition method was applied to the same system, and the equilibrium fraction of silver atoms contains all of the important thermal fluctuations. Among the contributions to the grand partition function included in the superposition method, the semigrand canonical basin-hopping method seeks only the term corresponding to the potential energy minimum that makes the largest contribution. Only at T = 0, when thermal occupation of higher energy minima is completely suppressed, do the two calculations coincide precisely. We show in Figure 7 the importance of such fluctuations by comparing the SGCBH and HSA results at 300 K, correcting the previous results by scaling the chemical potential difference, as explained above. In this graph, we have also superimposed the basin-hopping results obtained when the rotational contribution to the partition function is included.
The neglect of higher energy minima in the SGCBH method explains the staircase character of the silver fraction in Figures 6  and 7, whereas the superposition method gives smoother variations even at 100 K. As temperature increases, thermal fluctuations also increase and further smoothen the equilibrium silver fraction. 25 A similar effect is found for the SGCBH approach, which results from the increasingly similar Gibbs free energies of the various isomers. While the fluctuations predicted by the HSA appear relatively significant at 300 K, they are actually well represented by the global minimum of the semigrand partition function. It is also worth noting that some of the differences between the SCGBH and HSA calculations are related to the inclusion of nonicosahedral structures in the latter results, which contribute significantly to the stabilization of the composition 43:12 in Ag/Pd near Δμ = 1.2 eV. 25 This discrepancy observed between the superposition and basinhopping calculations would likely be reduced if the SGCBH simulations were no longer restricted to sampling the icosahedral homotop. Figure 6. Ag content in the lowest lying Ag n Pd 55−n icosahedron at a given (relative) chemical potential Δμ = μ Ag − μ Pd , found using SGCBH based on the potential defined in eq 5 for three temperatures. In the ball-and-stick representation of the most persistent stoichiometries, Ag atoms are gray (lighter) and Pd atoms are magenta (darker). Figure 7. Ag content in the lowest lying Ag n Pd 55−n icosahedron at a given (relative) chemical potential Δμ = μ Ag − μ Pd , found using semigrand canonical basin-hopping (SGCBH) based on the potential defined in eq 5 for T = 300 K, and compared to the results of harmonic superposition approximation (HSA) that include the contribution of a database of minima (solid line). The results of SGCBH accounting for the rotational partition function are also shown.

Journal of Chemical Theory and Computation
Article Finally, Figure 7 confirms that the rotational partition functions plays essentially no role on the segregation isotherms for this cluster. This result was expected because all minima are homotops of the two-layer Mackay icosahedron; hence, the moments of inertia differ only in the ways the atomic masses are allocated within this framework.
4.3. Some Remarks on the High Symmetry of Stable Compositions. The high symmetry of the most persistent compositions is consistent with the principle of maximum symmetry, 6,26,66 which suggests that structures with a higher symmetry content are most likely to lie in the low-and highenergy tails of the distribution. Since exceptions to this principle are not difficult to find, we do not expect it to be universally applicable, but there seems to be no doubt that many global minima observed for a wide variety of systems have high symmetry content. An explanation can be found by writing the total energy as a sum over contributions from a many-body expansion, involving single atom, pairwise and three-body terms, etc. If these terms are drawn from the same distribution, then geometrical symmetry (degeneracies) would be manifested as correlation. The variance is larger when correlation is present. Symmetrical structures are therefore more likely to have particularly high or particularly low energy. Low-lying structures are therefore likely to exhibit symmetry.
More formally, denote the mean and variance of a variable, X, drawn from probability distribution p(X) as μ and σ 2 . Contributions to the total energy are correlated when symmetry is present, and for the corresponding probability distributions, p(X,X′) ≠ p(X) p(X′). If we sum M terms with mean μ and variance σ 2 from the same distribution, then the variance of the mean is Mσ 2 + M(M − 1)ρσ 2 where the correlation ρ is defined by For ρ = 0, the variance is Mσ 2 , but for ρ = 1, it rises to M 2 σ 2 . On average, we therefore expect to find structures with a higher symmetry content in the tails of the distribution. 6 We have previously exploited symmetry-biased moves in the context of basin-hopping global optimization, and the core orbits scheme produced efficiency gains of 1 or 2 orders of magnitude for some benchmark atomic clusters. 67 The present results suggest that the maximum symmetry principle may also apply to nanoalloys.

CONCLUSIONS
Grand canonical ensembles are commonly encountered in problems relevant to nucleation or absorption. In this article, a methodology is introduced to search directly for the configuration that minimizes the Gibbs free energy and therefore makes the largest individual contribution to the grand potential. The method is based on basin-hopping global optimization and extends a recent effort 24 where the free energy landscape for a given system is minimized using the harmonic approximation for the entropy at finite temperature. Here, we have further developed the method to treat a variable number of particles (grand canonical ensemble) or, at fixed size, the composition in heterogeneous systems (semigrand canonical ensemble). Two applications illustrating both situations were presented, dealing with the nucleation of LJ clusters at fixed temperature or with the progressive segregation in nanoalloys as the chemical potential is varied.
In the grand canonical case, the method predicts that the global minimum in the local Gibbs free energy switches from dissociated states to bound clusters for larger chemical potentials and lower temperatures. Those results are in agreement with a specific model for the size-and temperature-dependent canonical partition function fitted to reproduce the heat capacities of reference clusters. They are also consistent with one of the primary conclusions of nucleation theory, namely, the existence of a critical nucleus that minimizes the Gibbs free energy. 59−62 Our semigrand canonical application to silver−palladium nanoalloys can be compared with the predictions of the harmonic superposition method, 25 which is similar to the present approach but incorporates the contributions of different minima to the grand partition function, not just the largest individual terms. The two approaches generally agree well, with discrepancies appearing at higher temperature, where new families of minima are stabilized, which are omitted in the present work.
Compared with the superposition method for the same statistical ensembles, grand canonical basin-hopping is less demanding because only the global minimum of the corresponding ensemble is sought instead of entire (ergodic) samples covering all relevant regions of the free energy landscape for the different system sizes or compositions. For the semigrand canonical ensemble, our approach provides a powerful way of solving the combinatorial problem of chemical ordering for a nonrigid lattice at finite temperature.
Future applications of the present methodology include the important case of absorption into porous materials, especially at low temperature and high densities, for which conventional grand canonical simulations are not practical. In particular, the contribution of zero-point motion can be straightforwardly incorporated in the expression of the partition functions from knowledge of the individual vibrational frequencies. Moreover, the notorious difficulty of sampling high-density phases would be alleviated owing to the systematic local minimization step of the basin-hopping procedure, thus increasing the chance of accepting the corresponding moves.
Funding F.C. acknowledges generous computational resources granted by the regional Pole Scientifique de Modeĺisation Numeŕique in Lyon. D.J.W. and D.S. acknowledge financial support from the EPSRC and the ERC.