The performance of Minima Hopping and Evolutionary Algorithms for cluster structure prediction
We compare Evolutionary Algorithms with Minima Hopping for global optimization in the field of cluster structure prediction. We introduce a new average offspring recombination operator and compare it with previously used operators. Minima Hopping is improved with a softening method and a stronger feedback mechanism. Test systems are atomic clusters with Lennard-Jones interaction as well as silicon and gold clusters described by force fields. The improved Minima Hopping is found to be well-suited to all these homoatomic problems. The evolutionary algorithm is more efficient for systems with compact and symmetric ground states, including LJ, but it fails for systems with very complex energy landscapes and asymmetric ground states, such as LJ and silicon clusters with more than atoms. Both successes and failures of the evolutionary algorithm suggest ways for its improvement.
To find the structural ground state of a cluster is a non-trivial global optimization task. One has to find the global minimum of the potential energy surface which is a function of all the atomic coordinates. Even for a relatively small cluster of 30 atoms the configuration space has already 90 dimensions. Because knowing the structure is a prerequisite for the study of all other physical and chemical properties the problem is of great importance and many algorithms have been developed to solve this global optimization problem. We compare Evolutionary Algorithms which have successfully been used in many diverse fields with the Minima Hopping Goedecker (2004) method.
For the prediction of the ground state structure of crystals, the USPEX (Universal Structure Predictor: Evolutionary Xtallography) Oganov and Glass (2006); Glass et al. (2006) method has turned out to be extremely powerful and has already allowed material scientists to find interesting and unexpected new crystal structures Gao et al. (2008); Oganov et al. (2006, 2008, 2007); Ma et al. (2007). Recently Evolutionary Algorithms have also been successfully used to predict surface phenomena such as steps on silicon crystals Chuang et al. (2005); Briggs and Ciobanu (2007). A widespread application of global optimization methods is the prediction of the structure of various clusters, In this field the majority of the work has been done with genetic or evolutionary methods as well Deaven and Ho (1995); Deaven et al. (1996); Wales et al. (); Johnston (2003). We note that different evolutionary algorithms developed for various types of structure prediction problems (molecules, clusters, crystals) have significant differences. Even for the same type of problem (e.g. crystal structure prediction) the previously proposed algorithms are very different in their construction and performance.
The presented evolutionary algorithm is similar to and inspired by the structure prediction algorithm USPEX. But since we work with non-periodic systems without a crystal lattice, some modifications of the original method were required. The version of minima hopping we are using is based on an improvement of the two key features of the original minima hopping method Goedecker (2004). The feedback mechanism is enhanced and the Bell-Evans-Polanyi principle Roy et al. (2008) is exploited in a more efficient way by moving preferentially along soft directions in the molecular dynamics (MD) part of the minima hopping algorithm.
Our comparison of Minima Hopping and Evolutionary Algorithms is based on Lennard-Jones systems, especially the cluster with atoms, which is an example of an easy one-funnel structure, and the -atom system, which is known to have a complicated double-funnel structure. We also apply the algorithms to more realistic systems, namely silicon clusters described by a force field and gold clusters described by an embedded atom potential.
This paper is structured as follows: We first introduce the evolutionary method used, after a quick introduction to Minima Hopping and its modifications we present the results section containing a comparison between Minima Hopping and the Evolutionary Algorithm and finally we also test different aspects of the EA and MH.
Ii The Evolutionary Algorithm
Evolutionary Algorithms (EA) implement a very simple model of biological evolution. They work on a set of samples — a population — which is gradually improved by selection and reproduction of fit members of the population — individuals. Each individual is a solution candidate. A single iteration step leads from a population to the next and is called a generation. The algorithm optimizes the fitness function — in our context the negative energy of the configuration. The operators applied to the population to obtain the next generation are the heart of the algorithm as they determine its quality and properties.
In contrast to the original simple genetic algorithm, modern applications in the field of chemical structure prediction all use real value encoding instead of binary strings and phenotypical operators acting directly in real space instead of gene modification Pullan (1997); Johnston (2003). Also state of the art is the application of local optimization to each individual thus reducing the search space to basin bottoms. Local optimization is done using standard techniques such as steepest descent and conjugate gradient methods.
FIG. 1 presents an overview of the Evolutionary Algorithm used. It starts with a randomly initialized population, this is our generation . In each generation the algorithm goes through three steps: Selection, Application of Operators and Acceptance. The steps Selection + Operator produce new offspring and mutations and put them into an intermediate population. In Acceptance the next generation is selected out of all available offspring together with individuals from the old population.
The algorithm possesses many tunable parameters. Most important values are populationsize, the number of offspring produced, the number of individuals kept from the last generation (elitism) and the mutationrate, for a complete list see TABLE 1.
|Number of offspring produced||offspring||populationsize|
|Number of individuals taken from former population||elitism||populationsize|
|Last rank with selection probability||cutoff||populationsize|
|Relative rate of offspring produced with average method||avgoffspring|
|Only one individual allowed within this energy interval||energyslot|
|Total rate of mutation||mutationrate|
|Random walk mutations (relative to total mutations)||mutrwalk|
|Strain mutations (relative to total mutations)||mutstrain|
|Probability of random rotation before recombination||raterndrot|
|Convergence criterion for force norm in local optimizer||fnrmtol|
As it is often the case in evolutionary algorithms in this field an energy slot restriction allows only one candidate per energy interval energyslot in the population. This method of preventing multiple copies of the same configuration in the population may be dangerous for it might reject an important candidate having almost the same energy as an individual already known. Using force fields allows to calculate energy and forces with very high precision since the numerical noise is extremely low. It is thus easily possible to identify structures by their energy.
Recently, a different structure of EA, optimized for parallel machines, has been presented Bandow and Hartke (2006). Instead of a stepwise evolution this approach handles a big pool of individuals which is subject to continuous application of operators. This is closer to a biological population without sharply defined generation gaps and it solves the load balancing problem in parallel implementations of EA.
Operators are used to evolve a population to a next generation. We use two different kinds of operators, heredity operators which take two individuals as input and produce a child sharing properties of both parents. The second kind is an operator applied to a single individual altering its configuration (mutation).
The selection step determines to which individuals the operators are applied, it is dependent on the operator. For a heredity operator there are two parents selected whereas for mutation operators only one individual is chosen. Selection is done using a linear ranking scheme. Individuals are sorted with respect to their fitness values and then assigned a probability depending linearly on the rank . The probability of selecting the individual with rank is in this case
where is the rank, starting at , the parameter cutoff and the first selection probability determined by normalization constraint. The cut-off value is the last rank with a selection probability above zero, all following ranks are assigned zero selection probability. The same method is also used in USPEX method Oganov and Glass (2006); Glass et al. (2006), it has turned out to be more efficient than Boltzmann selection where the selection probability follows a Boltzmann distribution depending on relative fitness values. The selection of the same individual serving as both parents is prevented. Mutation operators are applied randomly to the whole population.
The first heredity operator used is the cut and splice (cutting-plane) method introduced by Deaven and Ho Deaven and Ho (1995); Deaven et al. (1996). Both clusters are centered at their center of mass (COM). A randomly oriented plane cuts the two clusters apart. The new offspring cluster consists of one half of the first and the other half of the second cluster. Though able to obtain two offspring by this method we only produce one. The plane usually contains COM and its cut preserves the total number of atoms in the child cluster.
A new way of producing offspring is implemented in the average offspring method. Both clusters are centered at COM and for each atom of the first cluster the closest lying atom of the other cluster is identified. The atom of the child is now placed randomly on the connecting line of the two parent atoms. The randomness of this operator is necessary to prevent producing a lot of identical offspring.
Before application of either heredity operator a random rotation can be performed on one of the parents. The frequency at which this rotation is used can be adjusted via raterndrot.
Mutations are introduced to keep the diversity of the population high and prevent premature convergence. Three different methods are used. The easiest of those is the random walk mutation where atoms are randomly displaced. The displacement is approximately normal distributed with a mean displacement in the order of the two-body potential equilibrium distance (bond length). A strain mutation applies a geometrical deformation to the whole cluster, inspired by mechanical stress. The deformations include (anisotropic) compression and shear. Such strain transformations gave increased efficiency in the USPEX method. If the compression is high this method relates to the Big-Bang-Algorithm where configurations are relaxed from very high compressions Leary (1997). We only apply a moderate compression. In a third type of mutation the cluster is cut into two pieces similar to the plane-cutting method. One of the pieces is rotated around an axis perpendicular to the cutting plane by a random angle. This method was introduced as twinning mutation Wolf and Landman (1998).
Since our operators are able to produce configurations with atoms lying very close to each other we use a pre-relaxation method which is essentially a steepest descent with a very small step size. After a few steps the real self-adjusting steepest descent is started until a nearly quadratic region around the local minimum is reached. The final minimization is then carried out by a conjugate gradient method.
If an intermediate set of offspring has been created the Accepting step is triggered. In this step it is decided whether an offspring is accepted into the new population or is discarded. The algorithm first accepts the best individuals from the former population (elitism). The number of individuals chosen that way is given by elitism. In a second step the individual with the worst fitness value is replaced by the best offspring if energy slot constraints are fulfilled. This is repeated until all offspring are processed or the population is complete.
The algorithm is left running until a given limit of generations has been reached or the (known Wales et al. ()) global minimum has been found.
Iii Minima Hopping
Minima Hopping is a recently developed global optimization algorithm Goedecker (2004) which makes use of a Bell-Evans-Polanyi (BEP) principle for Molecular Dynamics (MD) trajectories Roy et al. (2008). The BEP principle states that low energy MD trajectories are more likely to enter the basin of a lower lying adjacent minimum than high energy trajectories. The algorithm also incorporates a history to repel it from previously visited regions. Using local optimization and molecular dynamics simulation it jumps between basins of attraction. The kinetic energy is kept as low as possible to escape the local minimum but is increased if this minimum has already been visited before.
Minima Hopping works with two self-adapting parameters, the kinetic energy of a MD escape step ekin and an acceptance threshold ediff for new minima to introduce a further downward preference. Starting from a local minimum a MD escape trial is started with kinetic energy ekin. After a few steps it is stopped and the configuration locally optimized again. If escaped the new minimum is only accepted when the new energy lies at maximum ediff higher than that of the previous minimum.
Minima Hopping is adjusted by tuning five feedback parameters: decreases ediff when a new minimum is accepted whereas increases the threshold on rejection, increases the kinetic energy when a MD escape trial fails, increases ekin when the new minimum is already known and decreases the kinetic energy if the minimum is unknown. Each visited minimum is added to a history list and marked as known. The algorithm currently uses the energy value to identify different minima. For more details we refer to the original paper Goedecker (2004).
Minima Hopping was used with the standard parameter set presented in the original paper: and .
The algorithm is efficient since it inhibits revisiting the same configuration many times which is likely to be the case in thermodynamically inspired methods such as Simulated Annealing.
We present two modifications of the original algorithm: The initial velocity vector of a MD escape trial is moved towards a direction with low curvature (softening) and a stronger feedback mechanism is used.
MD escape trials in the MH algorithm need an initial velocity distribution which is then
rescaled to fit the desired kinetic energy. The velocities are randomly
directed for each atom with Gaussian distributed magnitudes.
Regardless of the actual distribution chosen it has proved very useful
to use softening, to choose velocities along low-curvature directions.
In this way one can typically find MD trajectories with a relatively small
energy that cross rapidly into another basin of attraction.
In the original MH method low kinetic energy trajectories could only be obtained by
using large values for mdmin which results in long trajectories.
A direction of low curvature is found using a modified iterative dimer method
which only uses gradients, no second derivatives need to be calculated Henkelman and Jonsson (1999).
Starting at a local minimum with an escape direction the method calculates a second point at a distance along the escape direction. The forces are evaluated at and the point is moved along a force component perpendicular to :
After a few steps the iteration is stopped before a locally optimal lowest curvature mode is found. Initial velocities for the MD escape are then chosen along the final escape direction .
If the softening procedure is executed until it converges the performance drops again. It is important not to overdo softening. Always escaping into the same soft mode direction of a given minimum reduces the possibilities of different escape directions and therefore weakens the method. A good indicator was the mean kinetic energy during a run. For a few softening iterations the value decreases whereas it starts to increase again at a certain number of softening iterations. We set the iteration count to the value where the mean kinetic energy was minimal.
Typically iterations are done with a step size and a dimer length of .
In original MH ekin is increased by a factor if the current minimum has already been visited before — regardless of the number of previous visits. An enhanced feedback method uses a value of depending on the previous visits according to
where is the original value of and the number of previous visits to this minimum. The parameter has been set to after tests on bigger Lennard-Jones clusters and gold systems. This feedback mechanism reacts slightly stronger if the minimum is visited many times. If the system has only one energy funnel this enhanced feedback can even be slightly disadvantageous since it increases the kinetic energy too much and thus weakens the BEP effect of MD. The increased feedback mechanism improves the efficiency however considerably for large systems where the system can be trapped in huge structural funnels. If a cluster has for instance both low energy icosahedral and decahedral structures it takes a very long time for the MH algorithm without enhanced feedback to switch from one structure to the other.
Parallel Minima Hopping.
Parallelizing Minima Hopping is straightforward. On each processor an own MH process is started, all sharing the same history list. Only energy values of visited minima have to be shared. Due to the feedback mechanism and the common history list overlap of search areas is penalized in this parallel setup and running on several processors can thus yield an almost linear speedup in runtime. A further effect can be exploited in parallel runs. A single run might easily get trapped in a metastable funnel with a high escape time but the probability of all processors getting trapped in this local minimum is exponentially reduced with the number of processors running in parallel. The total expected runtime until success is therefore less influenced by the long escape times of relatively stable local minima in parallel runs. On the other hand there is a minimal number of local minima that have to be visited before the global optimum can be found. This minimal number of hops renders the use of too many processors less efficient again thus leading to an optimal number of processors depending on the structure of the PES. The idea of parallel sampling is known to have a positive effect Shirts and Pande (2001).
Iv Computational Experiments on Clusters
Atomic clusters are an ensemble of bound atoms, bigger than a dimer but smaller than bulk matter. They show interesting properties in the transition region of single atoms to bulk matter. From a global optimization point of view they are complicated multi-dimensional systems usually difficult to optimize since they contain a lot of local minima. Therefore they are of interest to test the capabilities of an optimization algorithm.
A simple model of chemical interaction of two non-charged atoms is the Lennard-Jones potential. The potential well depth and the equilibrium distance are the only parameters, both are set to 1.
Such models represent rare-gas clusters reasonably well. Lennard-Jones clusters are thoroughly studied model systems with well-known global minima up to atoms. Some of those clusters have a non-trivial multi-funnel potential energy surface (PES) where the global minimum is not easily found. The use of the Lennard-Jones two-body potential is fast compared to more accurate potentials or even Density Functional Theory (DFT) calculations. For these reasons they are well-suited to test global optimization methods. Our algorithms are applied to Lennard-Jones systems consisting of , , , and atoms. All clusters chosen show a different aspect, LJ is a very easy system, LJ has a double funnel structure, LJ is a very hard double-funnel system and LJ and LJ are examples of bigger clusters. Additionally, we perform single runs on LJ, LJ and LJ to find the ground state motif dependence of the algorithms. For these additional clusters we have pairs of similar-sized clusters with very different ground state geometries.
Silicon clusters are of more practical interest as silicon is so widely used in research and technology. To obtain realistic results usable in those fields it would be necessary to calculate energies using at least DFT methods. But those methods are consuming a lot of CPU time and one should highly optimize the algorithms applying DFT. Especially in the field of global optimization it is of importance to use an algorithm which is efficient to save computing time. To investigate the behavior of the presented algorithm we use only a force field to evaluate the energy of a configuration. Silicon systems are chosen as model systems since they possess directed bonds and show frustrated behavior. They have non-trivial minimum structures which are in general neither compact nor spherically symmetric. Silicon systems containing , , and atoms are explored. The force field used is Bazant’s Environment-Dependent Interatomic Potential (EDIP) Bazant and Kaxiras (1996); Bazant et al. (1997); Justo et al. (1998); Goedecker (2002). This force field has been chosen because it tends to elongated minimal structures which are not spherical (FIG. 2). There exist other force fields preferring spherical ground states. But the performance with spherical ground states is already investigated in Lennard-Jones and gold systems. Such highly elongated structures might not be very realistic as more accurate DFT calculation suggest only slightly elongated ground state configurations.
In gold systems different configurations can have very similar energy and the global optimum is often only very slightly lower than the next-best solution. This leads to a situation where the global geometry of the gold cluster can change completely by adding only a single atom. The energies of gold clusters are calculated using an empirical many-body potential by Rosato, GuillopÃ¨ and Legrand (RGL) Baletto et al. (2002). Gold clusters with , and atoms are investigated. Additionally, Au and Au are added to the test set in order to have different ground state geometries at similar sizes.
To compare different algorithms we measure two quantities, the total number of calls to the function calculating energy and forces and the total number of local geometry optimizations. The local optimization is the most expensive step in both algorithms and the speed is mainly determined by its number. To obtain meaningful numbers we perform between and runs on each problem.
The numerical convergence criterion for local optimizations is for the total euclidean force norm which leads to an energy precision of almost .
We gathered results for total performance on the investigated systems. Additionally, we also did comparative runs concerning the new modifications mentioned above. The numerical results of the different performance runs can be found in TABLE 3 whereas special test cases are addressed later. The table shows the results obtained using the best parameter set known to us. Those values are usually a combination of both heredity methods and all mutations. Minima Hopping found the global minimum in all problems investigated whereas the evolutionary algorithm could not find the global minima for LJ, LJ, Si, Si, Au and Au.
v.1 General Performance
The comparison shows a good applicability of Minima Hopping in all problems. The Evolutionary Algorithm is capable of finding most of the ground states but fails in some cases. Where it is successful it can even outperform MH. Problematic for the EA are the non-icosahedral ground states, such as LJ, LJ and the elongated silicon clusters.
The failure of the Evolutionary Algorithm to find the elongated ground state structure of silicon systems shows a clear advantage of Minima Hopping, it is not geometry dependent as the current Evolutionary Algorithms are. Moving via MD it is applicable to any system for which forces are available. In the past there have been approaches inspired by genetic algorithms which used standard crossover and in the beginning even binary strings as coordinate representation. Those ideas seem to be a bit less geometry dependent than the nowadays applied heredity methods. But the original simple genetic methods have been shown to be less effective than more elaborate geometrical operators in Lennard-Jones systems Pullan (1997).
When comparing two clusters of similar size but with different ground state structures it is obvious that the evolutionary algorithm can in general not reproduce non-icosahedral ground states with the operators presented here. The results of this tests can be found in TABLE 2. This table contains results of comparative runs without enough statistics to compare performance. The problem has already been identified and solved by niches to prevent domination of the whole population by only one geometry type Hartke (1999). Minima Hopping is able to find the non-icosahedral motifs without further modification in the standard configuration though with a decreased performance comparing to the icosahedral configuration of similar size.
The results show that a further adaption of EA to some specific features of clusters is necessary to enable the EA to work as efficiently for clusters as they do for periodic solids.
A further difference can be noted when comparing function calls. Minima Hopping tends to need in total more calls to the energy and forces function than the evolutionary algorithm in the cases where the EA succeeded, at least in smaller systems. This is clearly due to the use of MD and softening in Minima Hopping. But we remark that it is not necessary to perform MD escape with the full accuracy. A DFT calculation can be done using a reduced basis set and therefore in significantly less computer time. On the other hand Minima Hopping never produces energetically awful candidates since it moves via MD. Since heredity and mutation operators do not necessarily generate chemically reasonable configurations the local geometry optimization requires more force evaluations in the EA. Since in DFT applications the early stages of local optimization do not need high precision they could even be done using force fields thus reducing the computational demand.
|Cluster||Minima Hopping||Evolutionary Algorithm|
|Energy||#GO||calls (GO)||calls (MD)||runs||#GO||calls||configuration||runs|
v.2 Heredity Methods
While comparing both heredity methods used we observed that average offspring method performed better in systems with a compact optimal structure. The plane-cut method could only produce offspring as good as the other method by choosing a very low cutoff. When Deaven et al. Deaven et al. (1996) applied this method they produced every single combination of two parents from a very small population so having actually a very low cut-off level too. With the lower cut-off level the plane-cut method performed well overall.
The samples obtained using average offspring method are energetically closer to the fittest member of the population where the cutting plane method samples broader with more diversity (FIG. 3.a/b). It is a known fact that decreasing diversity in the population can lead to premature convergence. On the other hand, if samples in the complete energy range are allowed the method resembles too much a random search and loses efficiency. To observe a well evolving population it is necessary to have a balanced distribution of selected individuals. It seems therefore necessary to lower the cut-off parameter in plane-cut runs to select more of the fitter individuals as parents. A mix of both heredity methods delivered the best results (FIG. 3.c).
The plane-cut method is better suited to a general application of the algorithm, it can partially solve geometries less compact than the average offspring method.
In big clusters it proved useful to combine all operators available. The results, especially LJ, were best with a combination of all presented methods. In general the mixture was at a ratio or even more preferring plane-cuts in systems with known tendency towards non-spherical ground state (e.g. Si in TABLE 4).
We have also tested different plane-cut setups with a slightly modified method where a minimal distance between the two cluster halves is enforced. This method performed poorly and was always weaker than all different methods tested. Another modification where COM is not enforced to lie in the plane was also considered and dropped since there were no improvements.
The random rotation before recombination is necessary for average offspring method but only of advantage for the plane-cut method (see TABLE 5). If the system prefers non-spherical configurations raterndrot should be small.
v.3 Parameter Tuning
In systems with a double-funnel structure where the global maximum is located in the narrower funnel it might turn out advantageous to disable elitism, or (better) include in elitism only sufficiently different structures. In our tests LJ performed best when filling the population with offspring only. But we should remark that average offspring method has a rather preserving character often reproducing candidates already known.
A drawback of the evolutionary algorithm is the need to tune many parameters. A solution working with good performance on many different systems without adjustment would definitely be of interest. Minima Hopping performed better in this aspect, it never needed additional tuning, all runs were done using the same set of standard parameters. Using a standard set (TABLE 1) for all problems resulted in performance loss of the EA. The problems still converged but with performances down to half of the tuned versions shown in TABLE 3. A possible way to overcome this limitation is to fix the relative rates of elitism and mutation etc. and only adjust populationsize to the specific problem. Another possibility would be an automatically self-adapting version which tunes the parameters during runtime. In this case a stable and efficient scheme of parameter adaptation would be needed which is clearly not a trivial task.
We note that for crystals with up 30-60 atoms in the unit cell, the USPEX algorithm Glass et al. (2006); Oganov and Glass (2006) proved to perform very well with essentially a universal set of parameters without any parameter tuning. Evolutionary optimization of clusters, which are more complex systems, is more sensitive to parameter values.
v.4 Modification of Minima Hopping
Minima Hopping has been considerably improved using softening, in all studied cases.
We have tested an Evolutionary Algorithm capable of finding ground state structures of atomic clusters. In spite of the success of EA for periodic systems and on surfaces the current Evolutionary Algorithm is overall less efficient than Minima Hopping in the current implementation. It is not yet able to find global minima with geometrically difficult structures such as elongated silicon clusters and non-icosahedral ground states without the concept of niches. In contrast, Minima Hopping was able to find all ground states. Where the EA succeeds its performance is comparable to or even better than that of MH.
Further improvements of the EA could make it superior to MH for cluster optimization if the specifics of cluster systems are taken into account, as it was already done for periodic systems in the USPEX algorithm. The MH algorithm, on the other hand, shows a high stability and little need for further adaptation, at least for homoatomic systems.
Minima Hopping was improved by doing escape steps in directions with relatively low curvature of the PES and by using an enhanced feedback method.
- Goedecker (2004) S. Goedecker, J. Chem. Phys. 120, 9911 (2004).
- Oganov and Glass (2006) A. R. Oganov and C. W. Glass, J. Chem. Phys. 124, 244704 (2006).
- Glass et al. (2006) C. Glass, A. Oganov, and N. Hansen, Comput. Phys. Commun. 175, 713 (2006).
- Gao et al. (2008) G. Gao, A. R. Oganov, A. Bergara, M. Martinez-Canales, T. Cui, T. Iitaka, Y. Ma, and G. Zou, Phys. Rev. Lett. 101, 107002 (pages 4) (2008).
- Oganov et al. (2006) A. R. Oganov, C. W. Glass, and S. Ono, Earth Planet. Sci. Lett. 241, 95 (2006).
- Oganov et al. (2008) A. R. Oganov, S. Ono, Y. Ma, C. W. Glass, and A. Garcia, Earth Planet. Sci. Lett. 273, 38 (2008).
- Oganov et al. (2007) A. Oganov, Y. Ma, C. Glass, and M. Valle, Psi-k Newsletter 84, 142 (2007).
- Ma et al. (2007) Y. Ma, A. R. Oganov, and C. W. Glass, Phys. Rev. B 76, 064101 (2007).
- Chuang et al. (2005) F. Chuang, C. Ciobanu, C. Predescu, C. Wang, and K. Ho, Surf. Sci. 578, 183 (2005).
- Briggs and Ciobanu (2007) R. M. Briggs and C. V. Ciobanu, Phys. Rev. B 75, 195415 (2007).
- Deaven and Ho (1995) D. M. Deaven and K. M. Ho, Phys. Rev. Lett. 75, 288 (1995).
- Deaven et al. (1996) D. Deaven, N. Tit, J. Morris, and M. K., Chem. Phys. Lett. 256, 195 (1996).
- (13) D. J. Wales, J. P. K. Doye, A. Dullweber, M. P. Hodges, F. Y. Naumkin, F. Calvo, J. HernÃ¡ndez-Rojas, and T. F. Middleton, The cambridge cluster database, URL http://www-wales.ch.cam.ac.uk/CCD.html.
- Johnston (2003) R. L. Johnston, Dalton Trans. (2003).
- Ghasemi et al. (2008) S. A. Ghasemi, S. Goedecker, A. Baratoff, T. Lenosky, E. Meyer, and H. J. Hug, Phys. Rev. Lett. 100, 236106 (pages 4) (2008).
- Roy et al. (2008) S. Roy, S. Goedecker, and V. Hellmann, Phys. Rev. E 77, 056707 (2008).
- Pullan (1997) W. J. Pullan, Comput. Phys. Commun. 107, 137 (1997).
- Bandow and Hartke (2006) B. Bandow and B. Hartke, J. Phys. Chem. A 110, 5809 (2006), ISSN 1089-5639.
- Leary (1997) R. H. Leary, J. Glob. Optim. 11, 35 (1997).
- Wolf and Landman (1998) M. Wolf and U. Landman, J. Phys. Chem. A 102, 6129 (1998).
- Henkelman and Jonsson (1999) G. Henkelman and H. Jonsson, J. Chem. Phys. 111, 7010 (1999).
- Shirts and Pande (2001) M. R. Shirts and V. S. Pande, Phys. Rev. Lett. 86, 4983 (2001).
- Bazant and Kaxiras (1996) M. Z. Bazant and E. Kaxiras, Phys. Rev. Lett. 77, 4370 (1996).
- Bazant et al. (1997) M. Z. Bazant, E. Kaxiras, and J. F. Justo, Phys. Rev. B 56, 8542 (1997).
- Justo et al. (1998) J. F. Justo, M. Z. Bazant, E. Kaxiras, V. V. Bulatov, and S. Yip, Phys. Rev. B 58, 2539 (1998).
- Goedecker (2002) S. Goedecker, Comput. Phys. Commun. 148, 124 (2002).
- Baletto et al. (2002) F. Baletto, R. Ferrando, A. Fortunelli, F. Montalenti, and C. Mottet, J. Chem. Phys. 116, 3856 (2002).
- Hartke (1999) B. Hartke, J. Comput. Chem. 20, 1752 (1999).