Simulated Annealing

Consider the task of computing a thermodynamic ensemble average for a system with many degrees of freedom at low temperature (e.g., a large cluster, a polymer, or a protein).
The challenge presented by these many-body systems is that in addition to their global minumum energy configuration they usually have many local energy minima separated by high energy barriers. A Metropolis Monte Carlo computation at low temperature that starts from a configuration that is far from the minimum energy geometry usually leads to erroneous results. The reason for this is that the configurations that make the most important contributions to an ensemble average are those that are close to the minimum energy configuration and the algorithm is inefficient at sampling configurations that are beyond high potential energy barriers. Reliable Monte Carlo calculations thus require obtaining first the minimum energy configuration of the system.

The simulated annealing algorithm (by S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi Science 220 671-680, 1983) is an efficient technique to find the minimum energy configuration of the system. The algorithm can be described as follows:
Step (1): Initialize an arbitrary configuration $ \xi_0$ for the system at temperature $ T_0$, where $ T_0$ is reasonably large.
Step (2): Starting from the configuration $ \xi_0$, sample N configurations $ \xi_{1}, \xi_{2}, ..., \xi_{N}$ by implementing the Metropolis Monte Carlo method with $ P(\xi)=Z^{-1}$   exp$ (-E(\xi)/(k T_0))$.
Step (3): Go to (2), replacing $ \xi_0$ by $ \xi_{N}$ and $ T_0$ by a lower temperature.
Step (3) is repeated each time on the subsequent configuration $ \xi_N$ until the temperature of the system is equal to 0. It can be shown that the configuration that corresponds to the global minimum of $ E(\xi)$ can be reached according to such algorithm whenever the temperature decreases at a logarithmic rate (e.g., see S. Geman and D. Geman IEEE Transactions on Pattern Analysis and Machine Intelligence 6:721-741, 1984). In practice, however, a linear or even exponential temperature decrease schedule can often be implemented.

Exercise (due 04/24/03): This computational assignment has been designed and organized by Dr. Jose A. Gascon.

  1. Write a program for implementing the simulated annealing procedure and find the minimumm energy geometry of a ``cluster'' of two atoms interacting according to the 12-6 Lennard-Jones potential
    2. $\displaystyle U(r)=4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],$ (362)

    where $ \varepsilon=0.0104$    ev and $ \sigma=3.4$    Å. You can write you own code or modify the Fortran program mclj.for attached to the www page for the lecture notes of this course.
    Notes:

    1. If you decide to use the mclj.for program, you will have to edit the code and write a few lines as specified in the mclj.for file. The missing lines should specify the Metropolis Monte Carlo procedure.
    2. When running the mclj.for program you will be asked for the initial and final temperature. A reasonable value for the initial temperature is 10 K (just type 10). Since we are trying to find the global minimum the final temperature must be zero.
    3. When asked "Is the initial geometry random (yes/no)?", type "yes" to have the program select an initial guess of the geometry. You can eventually put your own initial guess in which case type "no". To create an initial geometry you must create a file called initial.xyz in the standard xyz format, where the first line is the number of atoms, the second line is left blank or with any coment, and the following lines have the atom type and coordinates as shown below for a cluster of N argon atoms.

      4. N
      comment
      Ar $ x_{1}$$ y_{1}$$ z_{1}$
      Ar $ x_{2}$$ y_{2}$$ z_{2}$
      ...
      Ar $ x_{N}$$ y_{N}$$ z_{N}$
    4. The mclj.for program reports on the screen the numbers of steps, the average interatomic distance, the energy at that step and the ratio of accepted trials out of 100.
    5. The final geometry is recorded in the output file final.xyz.
  2. Visualize the final geometry using a molecular viewer program such as Rasmol or Molden. The file movie.xyz contains snapshots of the accepted geometries along the simulation. The movie with the sequence of accepted geometries can be visualized by using Molden. Compare the minimum energy geometry (i.e., the minimum energy distance between the two Ar atoms) found by the simulated annealing algorithm with the exact value computed from the expression of the 12-6 Lennard-Jones potential.
  3. In order to appreciate the power of the simulated annealing method, find the minimum energy geometry of clusers with 3 and 13 argon atoms and report the values of the minimum energy. For the cluster with 13 atoms run the program with three different initial temperatures, 10 K, 20 K and 30 K. Compare the final results. Do the final energy and geometry depend on the initial temperature? Why, or why not?
  4. How would you compute a thermodynamic average at a constant temperature using the program for simulating annealing ?