J. Phys. Chem. 1992,96,4337-4341
4331
Performance of the Shift Method of Global Minimization in Searches for Optimum Structures of Clusters of Lennard-Jones Atoms Jaroslaw FWardy, Krzysztof A. Olszewski, and Lucjan piela* Quantum Chemistry Laboratory, Department of Chemistry, University of Warsaw, Pasteura 1 , 02-093 Warsaw, Poland (Received: June 6, 1991)
The optimal geometry of a finite cluster of atoms is difficult to find because the number of minima of the potential energy is huge for larger clusters. Until now, the only deterministic procedure designed to find the lowest-energy structures of such clusters was the diffusion equation method, which performed well for large clusters but sometimes failed for smaller ones. In the present paper we propose the so-called shift method (also deterministic), which deforms the original potential energy hypersurface by shifting the atomatom Lennard-Jones potential. As a result, the whole system collapses into a single point. Then, by using the reversing procedure, the original potential is gradually restored while optimizing the structure for each step. This results in a kind of big bang, because the atoms emerge from the collapse point and the cluster geometry evolves to the final one corresponding to the original undeformed potential. The procedure applied to the clusters containing up to 19 atoms turned out to produce in all cases under study the structures of the lowest energy.
1. Introduction In any theoretical approach to large and flexible molecules or to molecular clusters, one faces a very serious problem: that one is interested in a few low-energy conformational states of the system, while the number of all possible states is extremely large. The number of distinct conformations may be of the order of 10", which means that there is no hope to investigate them all and to select those of low energy. The problem is, of wurse, more general since it has to do with minimization of a function with multiple minima (the so-called global optimization problem'J). However, keeping in mind possible chemical applications, our particular purpose is to search for global minima of functions of very many variables with a large number of local minima. In such cases the task cannot be solved simply by making grid points because their number grows exponentially. In a series of papers3" we developed a deformation method designed to solve this problem. In the method one deforms the function to be minimized by increasing a certain deformation parameter in such a way as to smooth the fimction to be minimized and to induce fmt vanishing of some shallow bashs and then, when the parameter continues to grow, the deeper and deeper ones. During this procedure some other basins grow at the expense of those vanishing ones. An assumption of the procedure is that the global minimum basin survives until the very end, Le., until the moment when only one basin exists. There were indications in the literature that a smoothing of the potential energy hypersurface dramatically diminishes the number of its minima. Hoare and McInnes6.' have found that even a minor smoothing of the atom-atom potential (replacing the Lennard-Jones potential by the Morse potential) reduces the number of minima for the cluster of 13 atoms from 988 to 36. In practice, the above-described deformation is performed in a single step; Le., the deformation is suddenly turned on by setting a large value of the deformation parameter. With the deformed potential energy hypersurface, one easily gets the position of the minimum of the single basin by any minimization procedure, starting from any point of the multidimensional space. This feature is important, because here is the reason for the economy of the method. Then, this position is used as the starting point in the so-called reversing p r o c e d ~ r e .In ~ this procedure another minimization is performed which involves a slightly less deformed hypersurface, Le., the hypersurface corresponding to a slightly smaller value of the deformation parameter than that which produced the single-basin situation. This new minimization gives a new position of the minimum. The process is continued iteratively, and finally one gets the position of a minimum on the To whom correspondence should be addressed.
original (undeformed) hypersurface that hopefully is the global minimum position. Although the deformation method is quite general, in practical calculations we used only one particular def~rmation,~ in which one changes the original potential energy hypersurface Ax), where x is a set of the appropriate (e.g., conformational) variables, according to the diffusion equation (DEM, the diffusion equation method)
with the initial condition F(x,O) =Ax). The idea is especially transparent when applied to one dimension. Indeed, let us imagine a thin rod made of a homogeneous material extending along the x axis. Let the temperature distribution along the rod beAx) at time t = 0 (time t is the deformation parameter). Now, we let the time go on and the original shape of the temperature distributionflx) = F(x,O) changes to F(x,t) according to eq 1. The function F(x,t) becomes usually simpler and simpler as time goes on. In particular, its number of minima usually dramatically drops when t increases. We have shown3that the method works similarly as a high-frequency filter; Le., when applied to a Fourier representation of a function, the method makes high-frequency components vanish first (as t increases). The vanishing is faster when the frequency is higher, the damping factor being e&', where w is the frequency. This method has been used to describe some simple LennardJones systems? as well as for the so-called standard test functions for global o p t i m i z a t i ~ n . ~Recently, ,~ the method has been used to locate low-energy configurations of clusters of Lennard-Jones atoms.5 It turned out that the method finds the global minima in a majority of these applications. When applied to the standard test functions2 the method has found the global minima in all cases? and the computation time was as a rule much shorter than when other methods known in the literature were used. The exception represented the so-called Shekel where the corresponding global minima have been found, but the computation time was larger, the reason being that this was the only (1) Dixon, L.C. W., Szegb, G . P., Eds. Towards Global Optimisation I; North-Holland: Amsterdam, Holland, 1975; pp 29-54. (2) Dixon, L. C. W., Szegb, G . P., Eds. Towards Global Oprimisation II; North-Holland: Amsterdam, Holland, 1978; pp 1-1 5. (3) Piela. L.;Kostrowicki, J.; Scheraga, H. A. J . Phys. Chem. 1989, 93, 3339. (4) Kostrowicki, J.; Piela, L. J . Opt. Theory Appl. 1991, 69, 269. (5) Kostrowicki, J.; Piela, L.;Cherayil, B. J.; Scheraga, H. A. J. Phys. Chem. 1991, 95.41 14. (6) Hoare, M. R.; McInnes, J. Faraday Discuss. Chem. Soc. 1976,62, 12. (7) Hoare, M. R.; McInnes, J. Ado. Phys. 1983, 32, 791.
0022-3654/92/2096-433I$Q3.OQ/OQ 1992 American Chemical Society
4338 The Journal of Physical Chemistry, Vol. 96, No. 11, 1992 case when a numerical integration was involved. When applied to clusters of the Lennard-Jones atoms (up to 55 atoms), the method turned out to locate the global minimum configurations especially well for large clusters. For small clusters the method gave sometimes the next-lowest-energy configuration or some other low-energy one. For the number of atoms equal to 55, the number of the possible local minima was astronomical (of the order of 1045even without counting permutationally equivalent structures), but nevertheless, the method located the global minimum configuration correctly. The global minimum configuration is in this case known8-I0from some extensive local energy minimizations starting from very many configurations chosen in a systematic way (and assuming the atoms being of the same kind). Although the DEM has a clear mathematical meaning and is usually easy to use in practice, it may pose problems for some type of potentials. This may happen when their functional form is too complicated in order to get an analytical solution (recall the Shekel functions2g4)or when the potential contains poles (as e.g. the Lennard-Jones potential). In order to eliminate the unwanted (and unimportant) poles, we used in ref 5 an approximation of the Lennard-Jones potential by Gaussians. Note, however, that the main reason why the diffusion equation method works is that it deforms the hypersurface by diminishing its maxima regions and therefore forcing the basins to merge. One may think of some other deformation of that kind which may be simpler to apply in some cases. The aim of the present paper is to propose such a method that is applicable to a pairwise potential. We use the simple Lennard-Jones potential, because it is of the pairwise type and because the global treatment of the exact (e.g., ab initio) potential is prohibited due to huge computation time required to calculate the potential itself in a single point of space let alone its extremely complicated and unknown functional form. 2. Method Let us assume that the potential energy of a system composed of N atoms with positions given by a vector x may be represented as a sum of pair interactions atom-atom
V(x) =
xuij i 1, nothing would change. Then, in the reversing procedure the structure spreads out (the big bang) from the mllapse point and changes during the releasing of the deformation. When the number of steps in the reversing procedure is sufficiently large, the results are the same, Le., independent of the number of steps. In practice, we divided the whole range of the deformation parameter [0,1] into four parts, each part characterized by a different reversing procedure step. The step was smaller when the (I range was closer to 1, because the cluster sizes were smaller for a closer to 1 For each cluster a test was performed ensuring a decreasing step length keeps the final result unchanged. The I
b
I /
I
/
I / 3.8
f
01 333.144
1
I
150
0.5
>
0
t
a
0
Figure 3. (a) A typical trajectory of an argon atom in the DEM for a small cluster (Ar5);t is the deformation parameter: and x is a coordinate of an atom (A). The value o f t = 333.144 represents the collapse time.’ (b) A typical trajectory of an argon atom in the DEM for a large cluster and x is a coordinate of an atom (ArI3);t is the deformation ~arameter,~ (A). The value of t = 333.144 represents the collapse t h e e s (c) A typical trajectory of an argon atom in the shift method for small as well as large clusters (here is shown Ar12);a is the deformation parameter, and x is a coordinate of one of the atoms (A).
computation time varied from a few hours for the smallest cluster to a hundred hours for the largest one on the IBM 386 (33 MHz) personal computer. The results for the argon atoms are displayed in Table 11. As one can see in all cases studied the structures that emerged from the big bang represented those of the lowest energy. This demonstrates the power of the method, which was able to pick the global minimum out from a sometimes huge number of local minima. Only the case N = 8 appeared to be exceptionally difficult for numerical reasons; for the whole range of deformation parameters a very small step length had to be applied, which caused significant increase of computational time. Note that in fact a quasi-degeneracy occurs in that case, because two lowest-energy structures (see Table 11) differ in energy by about 0.3% in the 18-dimensional space, having the same number of locally favorable contacts between pairs of argon atoms. When the step of the reversing procedure was larger, the procedure indicated the other competing minimum (that one closest in energy). If the global minimum were not known in this case, applying our standard step length would probably result in the next-to-theglobal minimum instead of the global one.
J. Phys. Chem. 1992,96,4341-4346 TABLE III: Geometry Changes of the Argon Clusters during the Reversing Procedure no. of deformation cluster
atoms parametef 7 1.000 0.999 0.990
9
cluster geometry
0.0054 trigonal prism capped with one
half-octahedron pentagonal bipyramid pentagonal bipyramid point trigonal prism capped with two half-octahedra trigonal prism capped with two deformed half-octahedra inclined toward each other inclined trigonal prism capped with two deformed half-octahedra capped pentagonal bipyramid capped pentagonal bipyramid
0.999
0.061 6.1 0.0000 0.006
0.990
0.065
0.755
0.094
0.400 0,000 1.000 0.999
4.33
0.0000 point 0.0062 trigonal prism capped with three
0.350
4.84
0.000
7.45
0.000
8
size, A O.oo00 point
1.000
7.22
half-octahedra'' pentagonal bipyramid with two atoms added adjacently pentagonal bipyramid with two atoms added adjacently
"Note that this symmetry has been found by the DEM whereas it is the next-lowest-energy structure; see Table 11. It may be interesting to follow what happens just after the big bang. What happens to atoms just at the big bang during a single application of a local minimization procedure depends, of course, on the particular minimization algorithm (two minimization algorithms have been used: the steepest-descent method and the modified Cholesky factorization"). We have tested, however, that the minimizing configuration (at any step of the reversing procedure including the first one) is independent of the minimization algorithm chosen. It is of great importance to apply minimization algorithms that are able to detect saddle points and safely go out from there. During the reversing procedure the (11) Gill, P. E.; Murray, W. Math. Prog. 1974, 7, 311. (12) Mackay, A. L. Acto Crystallogr. 1962, 15, 916.
4341
geometry of the system changes considerably. As an example, we report the cases of N = 7-9 argon atoms (see Table 111). We do not report here the trajectories of the individual atoms that have been collected by us for the DEM as well as for the shift method. It is interesting, however, that there is a large difference in their character in both methods, even when the final result is the same. In the DEM for smaller clusters (up to Ar,) the trajectories are smooth and have a characteristic shape shown in Figure 3a. This shape corresponds to making first (after the big bang) the cluster extremely extended, e.g., 10 times the physical size, and then to shrink it to its physical size. For larger clusters (e.g., ArI3and larger) this picture becomes complicated, because the trajectory becomes sometimes suddenly irregular (Figure 3b). On the contrary, the shift method produced trajectories which are smooth and change monotonically for any cluster size (Figure 3c). 4. Conclusions Aggregates of atoms represent a challenge for methods aiming to search for their optimum geometrical structures, because the number of local minima of the function to be minimized is astronomical, e.g., of the order of 2 million for N = 19 and for N = 55. When counting the permutational symmetry present when atoms of the same kind interact, the numbers are to be multiplied by some other huge numbers of the order of lo1' and respectively. The problem of how to get low-energy structures is until now unsolved, although there is acquired a considerable numerical experience in starting from very many configurations of such atoms based on the so-called Mackay icosahedron and minimizing energy.8-1° The DEM represented until now the only deterministic procedure to locate the lowest-energy structure of the Lennard-Jones atoms; it failed, however, for four cases out of 14 studied.' The shift method of the present paper was able to reproduce the results of all those cases studied successfully by the DEM (except two, the most time-consuming ones being out of reach of the computers currently at our disposal) as well as to get the correct structures in all of the four cases where the DEM has failed. One may conclude, therefore, that the shift method shows a good performance in finding low-energy structures of clusters of the Lennard-Jones atoms.
Acknowledgment. This work was partly supported by a grant (No. BST-38) from the University of Warsaw.
Vibrational Levels of Water by the Self-Consistent-Field Method Using Radau Coordinates Jose Ziiiiiga,* Adolfo Bastida, and Albert0 Requena* Departamento de Quimica Fisica, Universidad de Murcia, 301 00-Murcia, Spain (Received: July 8, 1991)
Highly excited vibrational energy levels of the water molecule are calculated in the framework of the self-consistent-field (SCF) approximation using Radau coordinates to describe the vibrational motions. The calculations are performed for a global Sorbie-Murrell-type potential energy surface, which is well-behaved in both near-equilibrium and dissociating regions and for which adiabatic and exact quantum variational results up to very high energies are available for comparison. The discrete variable representation (DVR) method is used to solve the coupled SCF equations due to its simplicity and efficiency in computing the coupled multidimensional SCF integrals. The results obtained present a remarkable agreement with the exact variational values, mostly for states with increasing bending excitation, and are, in general, more accurate than the energy eigenvalues derived from adiabatic treatments.
I. Introduction The theoretical characterization of vibrational states is essential to the understanding of the dynamics and spectroscopy of polyatomic molecules. In the low-energy regime, the vibrational couplings are usually weak and the familiar normal modes approach provides a simple and excellent description of the molecular
motions. At higher levels of excitation, however, vibrations become greatly delocalized, anharmonic, and coupled, and solution of the vibrational Schrixlinger equation turns into a difficult task. In recent years, a considerable amount of effort has been devoted to the development of theoretical methods to deal with this problem.'"
0022-3654/92/2096-4341$03.00/0 0 1992 American Chemical Society
