 Next Article in Journal / Special Issue
From Metal Cluster to Metal Nanowire: A Topological Analysis of Electron Density and Band Structure Calculation
Previous Article in Journal / Special Issue
Structural transitions in biomolecules - a numerical comparison of two approaches for the study of phase transitions in small systems
Article

# Global Optimization by Adiabatic Switching

School of Physical Sciences, Jawaharlal Nehru University, New Delhi 110 067, India
*
Author to whom correspondence should be addressed.
Int. J. Mol. Sci. 2002, 3(1), 30-37; https://0-doi-org.brum.beds.ac.uk/10.3390/i3010030
Received: 14 September 2001 / Accepted: 10 October 2001 / Published: 31 January 2002

## Abstract

We apply a recently introduced method for global optimization to determine the ground state energy and configuration for model metallic clusters. The global minimum for a given N–atom cluster is found by following the damped dynamics of the N particle system on an evolving potential energy surface. In this application, the time dependent interatomic potential interpolates adiabatically between the Lennard–Jones (LJ) and the Sutton–Chen (SC) forms. Starting with an ensemble of initial conditions corresponding to the ground state configuration of the Lennard–Jones cluster, the system asymptotically reaches the ground state of the Sutton–Chen cluster. We describe the method and present results for specific cluster size N = 15, when the ground state symmetry of LJN and SCN differ.

## 1 Introduction

Determination of the lowest energy configuration for a cluster of N atoms is a nontrivial task [1, 2, 3, 4, 5]. The complexity arises in part from the exponentially (in N) large number of minima in the potential energy surface . Furthermore, the geometry of the potential energy landscape itself can make it computationally hard.
The problem is, however, simple to describe: Find the lowest energy minimum of a N body potential energy surface, $V ( r → )$,
$V ( r → ) = ∑ i , j V ( r i j )$
where $r →$ are the atomic coordinates, and $V ( r i j )$ is the interatomic potential of interaction between atoms i and j. For small N one can hope to enumerate all possible minima and decide the lowest of these, but even for moderate N and for the simplest V such as the Lennard–Jones (LJ) potential, this becomes difficult. A case in point is the 38 atom LJ cluster which has the so–called “double funnel” structure; the global minimum, which has octahedral symmetry, is marginally lower than the first excited state which has icosahedral symmetry. These were respectively found by the basin hopping technique  and a genetic algorithm method [8, 9].
A number of techniques of global optimization have been applied to this problem [1, 10, 11] and by now there are extensive compilations of global minima for a number of different clusters  , notably those described by two–body or many–body potentials which are commonly applied in atomistic simulations. A major difficulty is in ensuring that the algorithms reach the global minimum without being trapped in local minima. One method of overcoming such trapping  is by transforming the PES, broadening the thermodynamic transitions so as to increase the probability of finding the global minimum at temperatures where the free energy barriers are almost unsurmountable. For example, addition of a linear term to the PES provides a compressing effect which has been shown to be successful in locating the true minima in multiple funneled global structures . Locatelli and Schoen  used such transformations to locate the global minimum for non-icosahedral clusters.
We have recently proposed a new method of global optimization wherein time–dependence is introduced in the potential energy landscape . The evolving landscape is designed in a manner such that asymptotically, the potential energy surface develops into the surface of interest. A number of other techniques can be used to follow the evolving minima. Applications have been made to determining the ground state configurations of simple cluster systems .
In the present paper we apply this method to determine the ground state confi tions and energies of atomic clusters described by the many–body Sutton–Chen potential [17, 18] by switching from a known ground state. Initially, the interaction is chosen to be the Lennard–Jones potential,
$V 0 ( r → ) = ∑ i j V L J ( r i j ) = ∑ i j 4 ∊ [ ( σ / r i j ) 12 − ( σ / r i j ) 6 ]$
while the surface of interest is the potential
$V f ( r → ) = ∊ ∑ i [ 1 2 ∑ j ≠ i ( a r i j ) n − c ρ i ] , ρ i = ∑ j ≠ i ( a r i j ) m$
One choice for the time dependent potential energy surface is 
$V ( t ) = V 0 ( r → ) g ( t ) + V f ( r → ) h ( t )$
with g(t) an adiabatically varying switching function that interpolates between 1 and 0, and h(t) doing the reverse.
In the next section we describe the method as applied to the problem of ground state energy determination for Sutton–Chen clusters. Detailed results are presented for one cluster size, while the more general application and results are indicated in brief. This is followed by a discussion and summary.
It is a pleasure to dedicate this article to Steve Berry who has directly and indirectly influenced much of the development in the area of cluster studies over the past few decades. We have learned a lot from him, both in conversation as well as through his many articles and reviews .

The adiabatic optimization method , is a heuristic technique for locating minima. The essential idea is as follows.
Time dependence is introduced into the potential energy landscape directly by the incorporation of slowly varying terms as discussed in Eq. (4). A given choice is made for the switching functions g(t) and h(t), though in practice the choice does not affect the results greatly. A similar application of the adiabatic principle to determine semiclassical ground states of multidimensional systems  has noted the insensitivity of the technique to the precise form of the switching function, so long as the induced variation of the potential energy surface is slow enough. We note, parenthetically, that the switching principle has wide applicability, and in recent work has been used in the computation of the free–energy of finite clusters .
Location of the evolving minima can be done by any of a number of techniques. The simplest procedure is to introduce damping into the equations of motion and allow the system to evolve to a position of rest in a potential minimum; by starting with an ensemble of initial configurations and varying the available parameters, a number of minima can be located, and the putative global minimum can be recognized. Elsewhere  we have suggested the conjugate gradient  or simulated annealing (SA)  as other possible methods for locating the minima. It is likely that of these, the conjugate gradient technique will be more efficient as compared to SA though some SA variants  may also provide a suitable method for following the evolving minima.
The overall procedure can be summarized as follows:
• Take the initial configuration of the N atom cluster to be the ground state for the LJN cluster .
• Choose some switching function, say g(t). Similarly choose h(t), and the simplest choice, which we make here, is h(t) = 1 − g(t). We have explored a large variety of switching functions and in the present application we use g(t) = cos2(3πζt) exp(ζt), where ζ is the adiabaticity parameter.
• Perform molecular dynamics simulations for this N–particle cluster with forces deriving from Eq. (4), with an additional damping term, namely the equations of motion
$m r → k ̈ + γ r → k ̇ + ∂ V ∂ r → k = 0 , k = 1 , 1 … , N ,$
where $r → k$ is the position vector for the kth particle, m is its mass and γ is the damping constant.
• Vary ζ and γ, keeping in mind the natural timescales of the problem. Evolve to a minimum energy configuration, namely when the particle velocities become zero; the lowest energy found in an ensemble of simulations is the ground state energy predicted by the present method.

#### 2.1 Results for Sutton–Chen global minima

Here we attempt to switch from the minimum of the LJN system to the minimum of the SCN system. Both sets of minima have been extensively studied earlier and are tabulated in the Cambridge Cluster Database . A point of interest is that for the Sutton-Chen 9-6 family of potentials,[17, 18] the symmetries of the global minimum configurations are frequently different from the symmetries of the Lennard–Jones minima, so that in the adiabatic switching process, the cluster atoms must also move so as to adopt a different symmetry.
We present detailed results for the cluster size N = 15, though we have applied this technique to larger clusters and obtained results in agreement with the current standards . For the 15 atom LJ cluster, the ground state has the point group symmetry C2v while for the SC cluster the symmetry is D6d. Shown in Fig. 1 is a plot of the potential energy versus time for a particular choice of g(t), ζ and γ. Also shown is the effect of instantaneous switching, namely taking the limit g(t) = 0, where it can be seen that the system finds the nearest available local minimum from which it does not move. The time–dependence in the potential effectively permits the system to explore the multidimensional potential energy landscape of the SC cluster in an efficient manner. Finding a local minimum does not trap the system since there is always kinetic energy until the adiabatic switch is essentially over. Inset in the Figure is a schematic of the cluster configuration at different times during the process, showing how the cluster both contracts as well as rearranges to eventually reach the minimum of the SC surface.
It should be added that we have performed simulations for a variety of cluster sizes and in all cases we find that the procedure successfully finds the tabulated minima of SC clusters; these are not presented here since the details are repetitive. As we have emphasized elsewhere , the present method is heuristic, and thus some exploration of different switching functions, variation in the adiabaticity and damping parameters, and indeed the choice of initial potential, $V 0 ( r → )$ is necessary.

## 3 Summary and Discussion

In this paper we have presented the outline of a general procedure for global optimization with specific application to the problem of cluster ground state geometry determination. The application here, to the determination of the minimum of model metallic (Sutton-Chen) clusters by adiabatically deforming the potential energy surface relevant to model rare–gas (Lennard Jones) clusters is meant to be illustrative rather than exhaustive: the method introduced here is one of a class of techniques that employs time–dependence in the potential energy surface to enhance the exploration of phase space in contrast to other means of achieving the same objective .
A multiplicity of techniques is needed to approach hard problems such as global optimization. Few rigorous results are available, and application of most techniques is not guaranteed, with few possible exceptions, to give reliable (or certifiable) results. The present adiabatic switching method locally solves the optimization for an evolving surface, and thus mimics other methods of making large scale excursions in configuration or phase space.
We are presently studying this technique in detail with respect to the variation of parameters as well as to functional variations. One of the main issues of concern, and one that we are addressing in current work, is the relative efficiency of this method in comparison to other global optimization techniques. In a number of applications, we find that this method gives very encouraging results, and permits the determination of fairly reliable minimum energy configurations for a wide variety of cluster systems. The flexibility of choice of a number of starting potentials including the free
Figure 1. Plot of the potential energy versus iteration number, switching from the Lennard Jones potential to the Sutton Chen. The switching function used is g(t) = exp(−ζt) cos2(3πζt) with ζ = 0.4 and γ = 0.1. The time step is 0.01 in units natural to the LJ cluster, for which we also take = σ = 1. At time t = 0 the cluster has the C2v symmetry. At different times, as indicated, the cluster configuration is shown, and asymptotically, the configuration reached is the 9–6 Sutton Chen global minimum, with D6d symmetry. The parameters used for this latter model are taken from . The dashed line shows the result of the simulation in the absence of switching, namely when the LJ potential is suddenly transformed to the Sutton Chen potential.
Figure 1. Plot of the potential energy versus iteration number, switching from the Lennard Jones potential to the Sutton Chen. The switching function used is g(t) = exp(−ζt) cos2(3πζt) with ζ = 0.4 and γ = 0.1. The time step is 0.01 in units natural to the LJ cluster, for which we also take = σ = 1. At time t = 0 the cluster has the C2v symmetry. At different times, as indicated, the cluster configuration is shown, and asymptotically, the configuration reached is the 9–6 Sutton Chen global minimum, with D6d symmetry. The parameters used for this latter model are taken from . The dashed line shows the result of the simulation in the absence of switching, namely when the LJ potential is suddenly transformed to the Sutton Chen potential. particle case  as well as the flexible choice of switching functions and parameters, and finally the flexibility in the dynamical evolution all combine to suggest that while the method is heuristic, it holds promise.

## Acknowledgment:

This work is supported by a grant from the Department of Science and Technology. We thank Subir Sarkar for discussions.

## References

1. Wales, D.J.; Doye, J. P. K.; Miller, M. A.; Mortenson, P. N.; Walsh, T. R. Adv. Chem. Phys. 2000, 115, 1.
2. Wales, D. J. Atomic Clusters and Nanoparticles; Les Houches Session LXXIII; Guet, C., Ed.; Springer-Verlag: Heidelberg, 2001. [Google Scholar]
3. Wales, D. J.; Miller., M. A.; Walsh., T. R. Nature. 1998, 758, 394.
4. Doye, J. P. K. to appear in Global Optimization-Selected Case Studies; Pinter, J.D., Ed.; Kluwer: Dordrecht, 2001. [Google Scholar]
5. Wales, D. J.; Scheraga, H. A. Science. 1999, 285, 1368.
6. Stillinger, F. H. Phys. Rev. E 1999, 48, 59.
7. Doye, J. P. K.; Miller, M. A.; Wales, D. J. J. Chem. Phys. 1999, 110, 6896.
8. Deaven, D. M.; Tit, N.; Morris, J. R.; Ho, K. M. Chem. Phys. Lett. 1996, 195, 256.
9. Deaven, D. M.; Ho, K. M. Phys. Rev. Lett. 1995, 75, 288. [CrossRef] [PubMed]
10. Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A 1997, 101, 5111.
11. Hartke, B. J. Comp. Chem. 1999, 20, 1752.
12. Cambridge Cluster Database: http://brian.ch.cam.ac.uk/CCD.html.
13. Doye, J. P. K.; Wales, D. J. Phys. Rev. Lett. 1998, 80, 1357. [CrossRef]
14. Doye, J. P. K. Phys. Rev. E 2000, 62, 8753.
15. Locatelli, M.; Schoen, F. Comput. Optim. Appl. to be published. 2001.
16. Hunjan, J. S.; Sarkar, S.; Ramaswamy, R. to be published.
17. Sutton, A. P.; Chen, J. Phil. Mag. Lett 1990, 61, 139. [CrossRef]
18. Deyirmenjian, V. B.; Heine, V.; Payne, M. C.; Milman, V.; Lynden-Bell, R. M.; Finnis, M. W. Phys. Rev. B 1995, 52, 15191. [CrossRef]
19. See e.g. Berry, R. S. Nature. 1998, 393, 212.
20. Johnson, B. R. J. Chem. Phys. 1985, 83, 1204.Skodje, R.; Borondo, F.; Reinhardt, W. J. Chem. Phys. 1985, 82, 4611.
21. Miller, M. A.; Reinhardt, W. J. Chem. Phys. 2000, 113, 7035.
22. See e.g. Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B.P. Numerical Recipes; Cambridge University Press: Cambridge, 1992. [Google Scholar]
23. Kirkpatrick, S.; Gellat, C. D.; Vecchi, M. P. Science. 1983, 220, 671.
24. Wenzel, W; Hamacher, K. Phys. Rev. Lett. 1999, 82, 3003.Dittes, F.-M. Phys. Rev. Lett. 1996, 76, 4651.