ARTICLE pubs.acs.org/JPCA
Geometry Optimization of Atomic Clusters Using a Heuristic Method with Dynamic Lattice Searching Xiangjing Lai,* Wenqi Huang, and Ruchu Xu School of Computer Science and Technology, Huazhong University of Science and Technology, Wuhan, 430074, China
bS Supporting Information ABSTRACT: In this paper, a global optimization method is presented to determine the global-minimum structures of atomic clusters, where several already existing techniques are combined, such as the dynamic lattice searching method and two-phase local minimization method. The present method is applied to some selected large-sized Lennard-Jones (LJ) clusters and silver clusters described by the Gupta potential in the size range N = 13140 and 300. Comparison with the results reported in the literature shows that the method is highly efficient and a lot of new global minima missed in previous papers are found for the silver clusters. The method may be a promising tool for the theoretical determination of groundstate structure of atomic clusters. Additionally, the stabilities of silver clusters are also analyzed and it is found that in the size range N = 13140 there exist 12 particularly stable clusters.
1. INTRODUCTION The determination of the global minimum structure of molecular or atomic clusters by numerical global optimization methods is one of most challenging problems in computational chemistry. The difficulties arise mainly from the fact that the number of geometrically distinct local minima on the potential energy surface (PES) of the clusters grows exponentially with the cluster size.1 Therefore, developing highly efficient optimization methods is of great importance in tackling this kind of problem. In the past 20 years, a large number of successful global optimization methods have been developed for the problem. They can be divided into three classes according to their degree of dependence on the problem structure: (1) the totally problem-independent methods, (2) the totally problem-dependent methods, and (3) the partly problem-dependent methods. The methods in class (1) are general-purpose ones, where only objective function and its gradient need to be provided. The most representative examples include basin-hopping algorithm (BH) and its variants,25 minima-hopping algorithm (MHOP),6 and random tunneling algorithm (RTA).7 Although such methods can be extended easily to other multiple-minima problems, such as protein folding problem, their efficiency is usually lower than that of the problem-dependent methods. In the methods that belong to class (2), all operations depend strongly on the specific problem. Representative examples include the dynamical lattice searching method (DLS) and its variants811 and the heuristic algorithm with surface and interior operators (HA-SIO) and its variants.1214 Such methods are highly efficient, but they are also difficult to directly extend to other optimization problems. The methods in class (3) include conformation space annealing (CSA) algorithm,15 adaptive immune optimization algorithm r 2011 American Chemical Society
(AIOA),16 population basin hopping algorithm (PBH),17 and genetic algorithms (GAs),1823 etc., and some foundational ingredients in these methods are closely related to the specific problem, such as the dissimilarity measure in CSA, AIOA and PBH, and the cut-splice crossover operator in GAs.2123 The efficiency of these methods is generally between that of methods in class (1) and that of methods in class (2). Studies indicated that the methods in classes (1) and (3) are usually very successful in locating the global minima of small and medium-sized atomic clusters.26,1517 However, for relatively large atomic clusters it seems to be difficult to locate the true global minima within reasonable computational time by using such methods. In this case, the totally problem-dependent methods may be more practical. In previous studies,4,24,25 it has been shown that the PES of an atomic cluster consists of a lot of funnels and the difficulties of the corresponding global optimization arise mainly from two cases. In the first case, the PES of the cluster is very rugged, that is, there exist many large funnels on the PES and each large funnel contains a lot of small funnels. In the second case, the global minimum of the cluster lies in a very narrow funnel on the PES, and hence the system often enters the wider funnels which do not contain global minimum. As a result, the probability that the system enter the correct funnel is small. Based on the above considerations, we believe that a highly efficient optimization algorithm should be able to reach rapidly the bottom of the current funnel and to visit the narrow low-energy funnels with a relatively high probability. Received: November 6, 2010 Revised: April 17, 2011 Published: April 28, 2011 5021
dx.doi.org/10.1021/jp110620x | J. Phys. Chem. A 2011, 115, 5021–5026
The Journal of Physical Chemistry A
ARTICLE
In this paper, we present an efficient method for the geometry optimization of atomic clusters. In the method, several already existing techniques are combined to overcome the two difficulties mentioned above. At first, in order to overcome the first difficulty, the DLS method is adopted as a high-level local search routine, which can cross high energy barriers on the PES and locate rapidly the bottom of the current funnel without jumping out the funnel.26 To deal with the second difficulty, the seeding technique and two-phase local minimization method are employed to generate the starting structure of a cluster for DLS. Furthermore, the dynamic lattice searching method with interior operation (DLS-IO),10 which is an unbiased method, is also incorporated into the present algorithm to explore efficiently the whole configuration space. To examine the efficiency of the proposed method, we have applied it to the optimization of several selected large-size LJ clusters and silver clusters modeled by the Gupta potential in the size range N = 13140 and 300. Computational results show that the method is highly efficient, and a large number of new global minima missed in previous papers were found for the studied silver clusters.
2. THE POTENTIALS 2.1. Lennard-Jones Potential. Lennard-Jones (LJ) potential is one of most intensively studied pair potentials, which has been considered as a popular test system for evaluating the performance of the newly proposed cluster optimization algorithms. It can be written as 2 ! !6 3 12 N1 N σ σ 4 5 ð1Þ ELJ ¼ 4ε rij rij i¼1 j¼i þ 1
∑ ∑
where N is the number of atoms in the cluster, rij represents the distance between atoms i and j, and the reduced units with ε = 1.0, σ = 1.0 are used. The potential energy of the atom i in a cluster can be calculated according to the following eq 2: 2 ! !6 3 12 σ σ 5 ð2Þ ELJ ðiÞ ¼ 4ε 4 rij rij j6¼ i
∑
2.2. Gupta Potential. Gupta potential27 is a relatively reliable
and widely used potential in the optimization of the metal clusters, which can be written as EGupta ¼
UN N Vi 2 i¼1
∑
ð3Þ
0 1 1=2 rij r ij Vi ¼ A exp p 1 @ exp 2q 1 A r0 r0 j6¼ i j6¼ i
∑
∑
ð4Þ In the above equations, N is the cluster size, UN is a function of the cluster size, rij represents the distance between atoms i and j, and r0 is the equilibrium nearest-neighbor distance in the bulk metal. In the eq 4, Vi is composed of a pairwise repulsion energy of the BornMayer type and a many-body attractive contribution, the parameters p and q represent the repulsive interaction range and the attractive interaction range, respectively, and the
parameter A is fitted to experimented values of the cohesive energy.28 The following values of parameters were used for silver clusters in this work: A = 0.09944, p = 10.12, q = 3.37.29 The reduced units with UN = 1.0, r0 = 1.0 were adopted. The potential energy of the atom i in a cluster C can be calculated by the following equation, EGupta ðiÞ ¼ EGupta ðCÞ EGupta ðC \figÞ
ð5Þ
where EGupta(C) represents the potential energy of the cluster C and EGupta(C\{i}) represents the potential energy of the cluster obtained by removing the atom i from C. Note that the definition of atomic potential energy is important for the DLS method, and its roles in DLS will be stated in section 3.3.
3. COMPUTATIONAL DETAILS 3.1. Generation of the Initial Configuration. Since decahedral (D) structures are typical motifs for a wide variety of atomic clusters and are usually difficult to obtain in the optimization, an n 1 shell decahedron is adopted as the inner core of the initial structure of a cluster in the present method, where n represents the number of shells needed for locating all atoms of the cluster. The remaining atoms are randomly placed in the outside of the inner core, and the energy of the resulting cluster is then locally minimized by a local minimization procedure (the limited-memory quasi-Newton method (L-BFGS)30 was used throughout this study). Thus, the global-minimum structures with decahedral motif can be easily located. 3.2. Interior Operator. Generally, the energies of surface atoms are higher than those of interior atoms in a locally minimized cluster. Based on this consideration, Takeuchi proposed an interior operator in a previous study,12 where several surface atoms with high energy are moved near the center of mass of a cluster. In this study, a simplified version of the interior operator is presented, where the first m atoms with highest energy are randomly moved into a sphere having a radius of re/2. Here re is equilibrium distance between two atoms and the number m is randomly selected from 1 to 3. The center of the sphere coincides with the center of mass of the cluster. Note that the combined use of interior operator and local minimization is particularly excellent for lowering the energy of cluster, as mentioned in ref 12. 3.3. Dynamic Lattice Searching (DLS) Method. The DLS8 method starts from a randomly generated and locally minimized structure of a cluster. Then, a cycle of dynamic lattice construction, dynamic lattice searching and local minimization is performed until no new structure with lower potential energy can be found. The dynamic lattice construction procedure will find out all the stable locations on the surface of the cluster for an additional probe atom. These locations and Nmov (where Nmov is a predetermined integer) locations occupied by the atoms with highest energy are called dynamic lattice sites, and a set of these dynamic lattice sites is called a dynamic lattice. The dynamic lattice searching procedure will search for several low-energy candidates by iteratively moving the atom with highest energy to an unoccupied dynamic lattice site with lowest energy. These candidates are then locally minimized by means of a local minimization procedure, and the obtained best local minimum will be selected as the starting structure of the next iteration if it is lower in energy than the starting structure. In addition, it should be noted that the dynamic lattices are generally constructed by a series of minimizations of the energy of an added 5022
dx.doi.org/10.1021/jp110620x |J. Phys. Chem. A 2011, 115, 5021–5026
The Journal of Physical Chemistry A
ARTICLE
atom. A more detailed description of the DLS method can be found in refs 8,31,32. 3.4. Two-Phase Local Minimization. One two-phase local minimization involves two local minimizations.33,34 First, the first local minimization is performed with a modified potential function to guide the search to a promising region of the configuration space. Then, the second local minimization is performed with the original potential function starting from the point returned by the first local minimization to obtain corresponding local minimum. The modified potential function used in this work can be expressed as follows: Emod ¼ E þ
N1
N
½μ 3 rij þ β 3 ðmaxf0, rij2 ðη 3 RÞ2 gÞ2 ∑ ∑ i¼1 j¼i þ 1 ð6Þ
where N is the √ number of atoms, μ = 0.1, β = 3.0, η = 1.3 and R = re 3 [(3N)/(4π 2)]1/3. Here re is the equilibrium distance between two atoms, rij is distance between atoms i and j, and E is the original potential function. Equation 6 is a simplified version of the modified potential function proposed by Locatelli and Schoen for the optimization of LJ clusters.33 In fact, it can be found that the modified potential function presented in eq 6 favors the compact and spherical configurations since large distances between atoms are penalized. 3.5. Identification of the Central Atom. In this study, the central vacancy problem is also considered, since some globalminimum configurations of the LJ clusters and the metal clusters are without the central atom. Generally, it is difficult to identify the central atom of a cluster since the central atom is not necessarily the one nearest to the center of mass of the cluster. In order to efficiently identify and move the central atom, a definition called the locally atomic density is introduced here. The locally atomic density of atom i, D(i), can be written as N
∑
1 rij 3
ð7Þ j ¼ 1ð 6¼ iÞ rij e re where N is the number of atoms, rij is the distance between atoms i and j, and re is the equilibrium distance between two atoms. According to eq 7, it can be found that the atom with the maximal locally atomic density is usually the central atom of a cluster because the distances between the atoms in the inner layers are usually smaller than those between the atoms in the outer layers for a locally minimized cluster.35 3.6. Large Perturbation. To implement the structural transition of the cluster, a perturbation strategy is adopted in the present method, which can be stated as follows: (1) First, select the first m atoms with highest energy, where m is the maximal integer such that m e N/4. (2) Then, the selected atoms are randomly placed √ in a sphere with a radius of R/2, where R = re 3 [(3N)/(4π 2)]1/3 and the center of the sphere coincides with the center of mass of the cluster. (3) Finally, the resulting configuration is locally optimized by the L-BFGS method. After such a perturbation move, the structure of the cluster usually becomes disordered. 3.7. Description of the Present Algorithm (1) Generate the initial structure of a cluster by the method in section 3.1. (2) Perform the DLS procedure with the current structure as the initial structure. DðiÞ ¼
Figure 1. Flowchart of the present algorithm.
(3) Perturb the current cluster by the method in section 3.6. (4) A cycle of the interior operator and two-phase local minimization is carried out until no new structure with lower energy can be found during the last 10 attempts. (5) Perform the DLS procedure with the resulting structure as the initial structure. (6) Create a central vacancy for the current cluster, i.e., the central atom of the cluster is randomly moved on the surface of the cluster. (7) Perform the DLS procedure with the resulting structure as the initial structure. (8) Perturb the current cluster by the method in section 3.6. (9) A cycle of the interior operator and the L-BFGS method is carried out until no new structure with lower energy can be found during the last 10 attempts. (10) Perform the DLS procedure with the resulting structure as the initial structure.. The above algorithm was implemented in C language, and all the calculations were carried out on a PC with 2.4 GHz CPU and 1G RAM. Note that (1) the stop criterion for the first local minimization in the twophase local minimization method is |G| e 101, where |G| is the root mean-square (rms) gradient of potential functions; for other local minimizations used in this study, the stop criterion is |G| e 105; (2) the calculations will terminate once the known global minimum is found after DLS. The algorithm is summarized in Figure 1. The present method has two important objectives. The first objective is to combine the advantages of several existing methods to intensively exploit some typical low-energy funnels on the PES of 5023
dx.doi.org/10.1021/jp110620x |J. Phys. Chem. A 2011, 115, 5021–5026
The Journal of Physical Chemistry A
ARTICLE
Table 1. Comparison of the Performances of the Present Method, DLSc-RO and DLS-IO for Some Large-Size LJ Clusters DLSc-RO11
this work N
NLM
CPU time (h)
hit rate (%)
NLM
CPU timea (h)
DLS-IO10 hit rate (%)
b
500
1774
0.52
17
561b
451
0.18
54.5
660c
273
0.15
66
665d
13208
8.16
3
666b
9344
5.99
4
17608
667b
6404
4.46
6
6582
6.72
670b
2306
1.5
16
5188
6.76
2.8
685c 923b
504 2480
0.3 3.41
47.5 16
138 18976
0.16 75.2
70.2 1.0
NLM
CPU timea (h)
hit rate (%)
4163
2.63
2.5
9943
10.76
0.6
840
0.72
20.5
9302
15.09
0.8
144
0.14
65.2
124028
181.40
0.1
29582
38.03
0.5
152205
209.17
0.1
22.5
0.8 21208
43.02
0.6
2.4
a DLSc-RO and DLS-IO were performed on the computers with 2.8 GHz CPU and 1 G RAM. b Icosahedral motif. c Decahedral motif. d Central vacant icosahedral motif.
the clusters, respectively. In fact, the present method involves three typical methods, i.e., DLS with constructed core (DLSc),9 DLS-IO10 and two-phase local minimization. The DLSc method with decahedral inner core is able to efficiently exploit the decahedral funnels on the PES of the clusters; the combined use of interior operator, two-phase local minimization and DLS can guide the search into icosahedral (IC) funnels; and the DLSIO method is able to efficiently exploit close-packed funnels according to a previous study10 and our experience. On the other hand, the decahedral, icosahedral and close-packed structures are typical motifs for a wide range of atomic clusters, so it is useful to combine the above methods. The second objective of the present method is to achieve a trade-off between intensification and diversification. Intensification focuses on optimizing the objective function as far as possible within a limited search region while diversification should be able to drive the search to explore new promising regions of the search space.36 Study indicates DLS-IO is an unbiased optimization method which is able to explore the whole configuration space. Thus, the diversification of the search can be implemented by integrating DLS-IO into the present method. Moreover, for those clusters with a multifunnel PES, it is usually difficult to obtain the optimal configuration within a reasonable computational time using one simple method only, so the combination of the some existing simple methods is a significant way to construct a highly efficient optimization method.
4. RESULTS AND DISCUSSION 4.1. Optimization of Large LJ Clusters. To assess the performances of the present method, it was applied to the optimization of several large LJ clusters, including LJ500, LJ561, LJ660, LJ665667, LJ670, LJ685 and LJ923. The optimized results over 200 independent runs are listed in Table 1, together with the results of DLSc-RO11 and DLS-IO.10 The average number NLM of local minimizations and CPU time needed for a successful hit of the global minimum are usually adopted as the criteria to evaluate the speed of cluster optimization methods. From Table 1, it can be seen that NLM of the present method is obviously less than that of DLSc-RO for LJ500, LJ666, LJ670 and LJ923. As for LJ561, LJ660, LJ667 and LJ685, these values are comparable with those of DLSc-RO. Also, it can be easily found that NLM of the present method is much less than that of DLS-IO for all the clusters considered here. On the other
hand, the average CPU time per hit of global minimum for the present method is less than 9 h for all clusters under investigation. According to results of DLSc-RO, however, for LJ665, LJ666 and LJ923 the corresponding values are, respectively, 38.03 h, 22.5 and 75.2 h on the computers with 2.8 GHz CPU and 1 G RAM. This means that our method is in some sense faster and more robust compared with DLSc-RO. The success rate is another important parameter to examine the efficiency of optimization algorithms. From Table 1, it can be seen that the success rates for the present method are above 3% for all clusters considered here, and they are also higher than those of DLSc-RO and DLS-IO except for LJ685. For LJ500, LJ561, LJ670 and LJ923, the optimal configurations are based on icosahedral packing and the success rates of the present method are relatively high. The major reason is that two-phase local minimization used in this study favors compact and spherical configurations, and thus the clusters tend to form icosahedral structures in the optimization process. For LJ660 and LJ685, the success rates of the present method and DLSc-RO are in the same level, and the major reason is that the same strategy to generate initial structure of a cluster is used in both methods. On the other hand, from Table 1, it can easily be found that the global minima with icosahedral motif are not difficult to locate by the present method. This means that the perturbation strategy in section 3.6 is useful to implement the structural transition of the cluster, due to the fact that initial structures of the clusters have decahedral motifs. Moreover, the success of the present method on LJ665 which has an icosahedral motif without central atom indicates that the approach of creating the central vacancy for the cluster is obviously efficient. 4.2. Optimization of Silver Clusters with the Gupta Potential. To illustrate the applicability of the present method to the many-body potentials, we applied it to the silver clusters described by the Gupta potential containing up to 140 atoms. The present method located a large number of new global minima missed in previous papers and reproduced the previously putative global minima for the remaining clusters. The structures and energies of these new global minima are summarized in Table 2. For comparison, the previous results are also given in Table 2. It can be seen from Table 2 that the 16 new global minima, which are lower in energy than those reported in previous papers, were found in the present study, and their structures are different from those reported in the literature. For example, the structures 5024
dx.doi.org/10.1021/jp110620x |J. Phys. Chem. A 2011, 115, 5021–5026
The Journal of Physical Chemistry A
ARTICLE
Table 2. Comparison of the New Global Minima with Previously Putative Ones for Silver Clusters previous resultsa
new results N
a
energy
motifb
energy
motifb
previous resultsa
new results N
energy
60
62.5295
D
62.4991
IC
94
99.6872
66
69.0840
D
69.0746
D
95
100.8219
motifb FCC D
energy
motifb
99.6225
Dþ
100.7840
Dþ D
86
90.9531
FCC
90.8851
Dþ
96
101.8791
D
101.8538
88
93.0999
FCC
93.0812
FCC
98
104.1045
FCC
104.0593
D
89
94.1852
Dþ
94.1628
Dþ
100
106.3429
D
106.3212
D
90
95.2876
Dþ
95.2509
FCC
123
131.7289
D
131.7207
D
91
96.4109
FCC
96.3656
Dþ
132
141.6854
D
141.6770
D
93
98.5832
FCC
98.5700
Dþ
138
148.3256
D
148.3044
D
The previously putative global minima are from refs 37 and 38, respectively. b Dþ represents decahedral motif with an antilayer.
Figure 2. The finite difference ΔE of the energy of silver clusters.
of the new and previous global minima for Ag60 have decahedral and icosahedral motifs, respectively. For Ag86, Ag91, Ag93, Ag94 and Ag98, the optimal structures obtained in this study have facecentered-cubic (FCC) motifs, whereas the previous optimal structures are based on decahedral packing. For Ag66, Ag89, Ag96, Ag100, Ag123, Ag132 and Ag138, although the new and previous global minima have the same structural type, the arrangements of the outmost-shell atoms are different. Moreover, it can be easily found that in the size range N = 13140 the ratio of face-centered-cubic structures clearly increases compared with the previous results. Prior to this work, there are several global optimization methods that have been systematically applied to Ag clusters described by the Gupta potential, including RTA,28 DLS,37 DLS with constructed core (DLSc)38 and Aufbau/Abbau.39 Therefore, the fact that in this work the 16 new global minima were found clearly shows the effectiveness of the present method. On the other hand, the average CPU time and the success rate are two important parameters to measure the efficiency of optimization methods. In this work, the average CPU time needed for a hit of global minimum is less than one hour for most clusters in this size range, and the corresponding success rate is ca. 10%. However, it should be noted that there are still a few cases are difficult to optimize, such as Ag41 and Ag96. For these clusters, a considerable computational effort is required for searching their global minima. In order to further illustrate the efficiency of the present method for large silver clusters, we have optimized the structure of Ag300. The global minimum obtained in this study has the
Figure 3. The second finite difference Δ2E of the energy of silver clusters.
energy of 329.0460, which is 0.1471 lower than the previous global minimum obtained by DLSc.40 The average CPU time needed for a hit of global minimum over 200 independent runs is approximately 12 h, and the success rate is ca. 2.5%. Moreover, it took about 0.3 h for an independent run on average. In DLSc,40 however, it took 38 h for an independent run for those clusters from Ag170 to Ag310; furthermore, for some clusters, it took more than 1000 h for a successful hit of global minimum.40 In contrast, the present method is higher in terms of both the solution quality and the computing speed. There are two possible reasons that can be used to explain why the present method outperforms the DLSc method used in refs 38 and 40. At first, the definition for the energy of a single atom in a cluster, EGupta(i), is particularly suitable for DLS in the optimization of metal cluster modeled by the many-body potentials. With the definition in eq 5, the DLS procedure performs well in our implementations. Second, the present method is essentially unbiased and can explore the whole configuration space, while the DLSc method will become biased if the size of cluster is large. 4.3. Stability Analysis of Silver Clusters. To illustrate which clusters are particularly stable compared with their neighbors, the finite difference ΔE of the energy and the second finite difference Δ2E of the energy are plotted as a function of the cluster size in 5025
dx.doi.org/10.1021/jp110620x |J. Phys. Chem. A 2011, 115, 5021–5026
The Journal of Physical Chemistry A
ARTICLE
Figures 2 and 3, respectively. The ΔE and Δ2E are defined as follows: ΔEðNÞ ¼ EðNÞ EJ ðNÞ
ð8Þ
Δ2 EðNÞ ¼ EðN þ 1Þ þ EðN 1Þ 2EðNÞ
ð9Þ
where EJ(N) = a þ bN1/3 þ cN2/3 þ dN is a four-parameter fit of the energy of global minimum. The valleys in ΔE and the peaks in Δ2E correspond to the particularly stable clusters compared with their neighboring clusters. For silver clusters in the size range N = 13140, EJ(N) = 0.4806 þ 0.414N1/3 þ 0.4136N2/3 1.167N, according to the results obtained in this work. At first, it can be easily found from Figure 2 that there exist 8 apparent valleys at N = 13, 38, 55, 64, 71, 75, 86, 101, showing the stability of the corresponding clusters. Then, from Figure 3, we can find 12 apparent peaks at N = 13, 38, 55, 64, 71, 75, 86, 101, 108, 117, 126, 135. Thus, we can conclude that in the size range of N = 13140, there are 12 particularly stable clusters that appear at N = 13, 38, 55, 64, 71, 75, 86, 101, 108, 117, 126, 135. These results are basically consistent with those reported in a previous paper38 except for Ag86.
5. CONCLUSIONS Based on several already existing optimization techniques, we have designed a heuristic algorithm for geometry optimization of atomic clusters in this work. By application of the algorithm to the optimization of large LJ clusters and silver clusters described by the Gupta potential, it is shown that the algorithm is highly efficient. Using the present method, we have found a large number of new global minima missed in previous papers for silver clusters under investigation. Moreover, it is found that in the size range of N = 13140 the ratio of face-centered-cubic structures clearly increases compared with the previous results reported in the literature.38 ’ ASSOCIATED CONTENT
bS
Supporting Information. Cartesian coordinates and energies of the silver clusters obtained in this work. This material is available free of charge via the Internet at http://pubs.acs.org.
(7) Jiang, H.; Cai, W.; Shao, X. Phys. Chem. Chem. Phys. 2002, 4, 4782. (8) Shao, X.; Cheng, L.; Cai, W. J. Comput. Chem. 2004, 25, 1693. (9) Yang, X.; Cai, W.; Shao, X. J. Comput. Chem. 2007, 28, 1427. (10) Shao, X.; Yang, X.; Cai, W. J. Comput. Chem. 2008, 29, 1772. (11) Wu, X.; Cai, W.; Shao, X. Chem. Phys. 2009, 363, 72. (12) Takeuchi, H. J. Chem. Inf. Model. 2006, 46, 2066. (13) Takeuchi, H. J. Phys. Chem. A 2008, 112, 7492. (14) Takeuchi, H. J. Chem. Inf. Model. 2007, 47, 104. (15) Lee, J.; Lee, I.; Lee, J. Phys. Rev. Lett. 2003, 91, 080201. (16) Cheng, L.; Cai, W.; Shao, X. Chem. Phys. Lett. 2004, 389, 309. (17) Cassioli, A.; Locatelli, M.; Schoen, F. Optim. Methods Software 2009, 24, 819. (18) Deaven, D. M.; Tit, N.; Morris, J. R.; Ho, K. M. Chem. Phys. Lett. 1996, 256, 195. (19) Wolf, M. D.; Landman, U. J. Phys. Chem. A 1998, 102, 6129. (20) Hartke, B. J. Comput. Chem. 1999, 20, 1752. (21) Chen, F. Y.; Curley, B. C.; Rossi, G.; Johnston, R. L. J. Phys. Chem. C 2007, 111, 9157. (22) Deaven, D. M.; Ho, K. M. Phys. Rev. Lett. 1995, 75, 288. (23) Johnston, R. L. Dalton. Trans. 2003, 22, 4193. (24) Doye, J. P. K.; Miller, M. A.; Wales, D. J. J. Chem. Phys. 1999, 110, 6896. (25) Doye, J. P. K.; Wales, D. J.; Miller, M. A. J. Chem. Phys. 1998, 109, 8143. (26) Cheng, L.; Feng, Y.; Yang, J.; Yang, J. J. Chem. Phys. 2009, 130, 214112. (27) Gupta, R. P. Phys. Rev. B 1981, 23, 6265. (28) Shao, X.; Liu, X.; Cai, W. J. Chem. Theory Comput. 2005, 1, 762. (29) Michaelian, K.; Rendon, N.; Garzon, I. L. Phys. Rev. B 1999, 60, 2000. (30) Liu, D. C.; Nocedal, J. Math. Prog. 1989, 45, 503. (31) Cheng, L.; Cai, W.; Shao, X. Chem. Phys. Chem. 2007, 8, 569. (32) Cheng, L.; Cai, W.; Shao, X. Chem. Phys. Chem. 2005, 6, 261. (33) Locatelli, M.; Schoen, F. Comput. Optim. Appl. 2003, 26, 173. (34) Doye, J. P. K.; Leary, R. H.; Locatelli, M.; Schoen, F. INFORMS J. Comput. 2004, 16, 371. (35) Liu, H. H.; Jiang, E. Y.; Bai, H. L.; Wu, P.; Li, Z. Q. Chem. Phys. Lett. 2005, 412, 195. (36) L€u, Z.; Huang, W. Phys. Rev. E 2009, 80, 026130. (37) Zhan, H.; Cheng, L.; Cai, W.; Shao, X. Chem. Phys. Lett. 2006, 422, 358. (38) Yang, X.; Cai, W.; Shao, X. J. Phys. Chem. A 2007, 111, 5048. (39) Alamanova, D.; Grigoryan, V. G.; Springborg, M. J. Phys. Chem. C 2007, 111, 12577. (40) Shao, X.; Yang, X.; Cai, W. Chem. Phys. Lett. 2008, 460, 315.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT The authors thank referees for their comments and questions which help us to improve the paper. This work was partially supported by the National Natural Science Foundation of China (Grant No. 60773194 and 61070235). ’ REFERENCES (1) Doye, J. P. K. Phys. Rev. E 2000, 62, 8753. (2) Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A 1997, 101, 5111. (3) Wales, D. J.; Scheraga, H. A. Science 1999, 285, 1368. (4) Leary, R. H. J. Global Optim. 2000, 18, 367. (5) Zhan, L.; Chen, J. Z. Y.; Liu, W. K.; Lai, S. K. J. Chem. Phys. 2005, 122, 244707. (6) Goedecker, S. J. Chem. Phys. 2004, 120, 9911. 5026
dx.doi.org/10.1021/jp110620x |J. Phys. Chem. A 2011, 115, 5021–5026
