Modern strategies to minimize a function in several unconstrained variables were defined in the 1970s. These methods fall into the general category of “pseudo-Newton” approaches and make extensive use of the gradient of the function, as the only way to retain computational efficiency with large numbers of variables. The minimizers in AMPAC™ (BFGS, EF, DFP, TRUSTE, or NEWTON for energy and TS, TRUSTG, or LTRD for gradient) are all local minimizers in that they converge toward the minimum (or critical point) that is in the most favorable position when related to the starting point. This may not be the global minimum, but only the one that is lowest in energy in a specific region of the potential energy surface. When multiple minima do exist on the potential surface, it is often difficult to locate them directly. A tedious process of trial and error is usually followed by defining a variety of starting points around the area of viability for the optimizable parameters (a grid search). Thus, an open problem in theoretical optimization, the location of the absolute minimum of a function in many unconstrained variables, has direct application to one of the most persistent and annoying problems in computational chemistry, the “multiple-minima” dilemma. In 1983, Kirkpatrick et al. suggested a heuristic answer to this problem, making use of a method of simulation in statistical mechanics coupled with Boltzmann’s law at decreasing temperatures.^[36] “Simulated annealing” was the name given the procedure.

A reasonable question for computational chemistry is not necessarily finding the minimum solution of the energy/geometry function defining the molecule, but collecting a variety of possible solutions, i.e. local geometric minima. If properly structured and implemented, the simulated annealing process can be used to collect an entire set of minima. The simulated annealing procedure is very slow in refining the exact location of a minimum, and specialized local methods (much more rapid) are used for exact location. A mixed strategy can be derived in which the approximate locations of the minima are generated by the annealing procedure and the final refinement is carried out from each starting place by one of the local optimization protocols mentioned above.

The theory and terminology behind simulated annealing comes from statistical mechanics (the behavior of systems with many degrees of freedom in thermal equilibrium at a finite temperature). The idea of annealing is to heat an object (e.g. a piece of metal) to a high temperature and then allow the object to slowly cool down. The high temperature allows the atoms to move from their original positions allowing the system to rearrange into a more stable state. As the system cools, it becomes trapped in its new configuration and the system will continue to evolve toward the lowest energy configuration in that local region. Slow cooling may also be followed by a rapid cooling (“quenching”) which removes any residual heat, effectively locking the atoms in their new configuration.

Annealing can used to move a macroscopic system toward its most stable configuration (global minimum) and so has elicited interest in how it can be exploited to solve various types of problems. (In real systems, however, multiple heating/cooling cycles may be required and one is never certain that the true optimal configuration has been achieved.) In typical macroscopic systems (solid or liquid) the density of atoms is on the order of 10²³ atoms per cubic centimeter. For such a system in thermal equilibrium, statistical mechanics shows that the behavior of the system is dominated by the most probable behavior and small fluctuations from this average behavior can be neglected. In this type of ensemble, each configuration, defined by a set of atomic positions {r_i}, of the system is weighted by its Boltzmann probability factor, exp(-E{r_i }/k_BT), where E{r_i } is the energy of the configuration, k_B is Boltzmann’s constant, and T is temperature.

This Boltzmann factor was used to develop the Metropolis step,^[37] a simple algorithm that could be used to efficiently simulate a collection of atoms in equilibrium at a given temperature. At each step, a new geometry is generated as a random displacement from the current geometry. The energy at the new geometry is computed and the difference in energy between the current and the new geometries is given as ΔE. The probability that this new geometry is accepted, P(ΔE), is given by Metropolis’ criterion:

(13.1)

Thus, if the new geometry is lower in energy than the current geometry, the step is always accepted. If the step leads to a higher energy, the step will only be accepted if a randomly generated number (between 0 and 1) is less than or equal to the Boltzmann probability factor. At high temperatures, the Boltzmann factor is close to one, so uphill steps (increasing energy) are frequently accepted allowing the system to climb up out of its local valley. At lower temperatures, the Boltzmann factor drops toward zero causing frequent rejection of uphill steps. Thus the system is prevented from escaping its local region and is progressively driven down in energy.

The simple annealing procedure just defined is useful for the global minimum problem but is not well suited for finding multiple minima (or critical points). Also the final stage of annealing converges very slowly and so is inefficient for our purposes. Traditional energy (or gradient) minimizers are much more efficient for the purpose of converging on local extrema. This suggests that the best approach is a heuristic algorithm that effectively mixes the rapid global searching capability of annealing with the efficient local optimization of traditional minimizers. In AMPAC, this mixed annealing strategy is accomplished in three distinct steps.^[37]

Simulated Annealing

Each new geometry is generated as a random displacement from the current geometry. Displacements are generated randomly but must be scaled so that they fall within a boundary hypersphere of radius, R. To do this, a random number is generated between 0 and R and the length of the generated displacement is scaled to match this value. By default, a uniform random number is used for this step. Keyword GAUSSIAN forces the use of gaussian random number generator instead. (GAUSSIAN only affects the scaling of the displacement, not the generation of the displacement.) This generated displacement is then added to the existing geometry to produce a new displaced geometry. The energy is computed at this new geometry and is accepted as the new current geometry if it passes the Metropolis criterion (described above). (Note, that the Metropolis criterion can be generalized as will be discussed under the individual simulated annealing methods given later.) When a step is rejected, a new displacement from the current geometry is generated and tested until a successful step is found.

The first step of the simulated annealing phase is to “melt” the system by performing a Metropolis walk on the system at a high initial temperature, T_init. This walk is continued until the system reaches a steady state (“thermalization”). Subsequently, the temperature is lowered and again allowed to run until thermalization is achieved. Each new temperature, T_n, is determined by multiplying the previous temperature by factor, T_dec, where 0 < T_dec < 1. Thus the n^th temperature used is given by:

(13.2)

This process of alternately lowering the temperature and running until thermalization is achieved is continued until the system “freezes” and no further changes occur. (T_init and T_dec are set by keywords TEMP=n.n and TLAW=n.n respectively.)

The boundary radius, R, used in generating displacements is adjusted at each temperature, T_n. As a rule of thumb, the radius, R_n, is chosen to give a rejection ratio of approximately 0.5, where the rejection ratio is the number of points rejected by the Metropolis criterion divided by the total number of points. This adjustment to R_n can be handled automatically but is upper bounded by some value, R_max, to avoid possible excessive steps. A minimum boundary, R_min, can also be specified as f_min * R_max. If the radius R_n falls below R_min, the system is considered “frozen” and the annealing phase is ended. (R_max and f_min are set by STEP=n.n and STEPCV=n.n respectively.)

The Metropolis walk is continued at a given temperature until “thermalization” is achieved or the maximum number of steps, N_max, is exceeded. A fairly simple thermalization criterion that can be easily computed is used for this purpose.^[38] At every N_check points, E_σ (standard deviation of energy), ē (average energy), and e_min (minimum energy) is computed for that set of points and tested against our thermalization threshold, σ_therm. This thermalization criterion is given as:

(13.3)

(N_max, N_check, and σ_therm are set by NMAX=n, NCHECK=n, and STD=n.n, respectively.)

Metropolis walks are performed at successively lower temperatures until the system is deemed to be “frozen” (trapped in a local basin) and the simulated annealing phase is ended. The system is deemed frozen if (a) the boundary radius, R_n, is less than R_min; (b) the temperature falls below a certain threshold (usually T_init/100); or (c) the final geometry of two or three Markov chains is deemed to be the “same.” By “same” here, we mean that the “distance” between the geometries falls below some threshold.

In AMPAC, simulated annealing has been generalized to include four related strategies, each with their own strengths and weaknesses. The original Metropolis step gives the probability of accepting the step in terms of the change in energy during the step. To apply the Metropolis step and the simulated annealing strategy to more general problems, we can generalize from change in energy to a change in a generalized “cost function.” For chemical systems, this opens the door to using the RMS gradient norm (gnorm) as the cost function in addition to the traditional energy (heat of formation). An additional generalization is that we can add penalty functions to the basic energy or gnorm. TSANNEAL, GANNEAL, and MANNEAL include an energy window penalty term that is directly added to the energy/gnorm used in the Metropolis step. This focuses the annealing search on energies in the window centered at F_ref and extends by S_ref in either direction. Configurations with an energy outside this window are penalized by a quadratic function with coefficient P_ref.

(13.4)

(See annealing keywords, FREF=n.n, SREF=n.n, and PREF=n.n.)

Most annealing algorithms and other multiple-minima search procedures are applied in the context of molecular mechanics (MM). This is because a large number of independent energy calculations usually need to be carried out and MM is the only chemical modeling method efficient enough to do this. To minimize the number of such calculations required using the higher quality, but much slower semiempirical methods, the researcher can apply certain additional constraints to the preliminary search (prior to the use of semiempirical energy calculations), based on his chemical understanding of the problem. These constraints are treated as penalty functions to the genuine search criterion (energy (ANNE, TSANN), gradient norm (GANN), or both (MANN)). Enormous savings and enhanced efficiency can result. For example, in conformational problems involving aliphatic rings, up to 99% of the configurations generated correspond to geometries that are so poor that SCF convergence will be difficult. These conformations are discarded by the annealing procedure at a cost of only 1% of the computer time. The constraints presently implemented are

While the annealing procedure in AMPAC is designed to somewhat automate the search for minima (and other critical points) on a potential energy surface, it is not a “black box.” The judicious use of geometry input and the various penalty functions can make the search both complete in the quantum mechanical potential energy space and efficient computationally. It is neither chemically correct nor necessary to search ALL of the potential energy surface, when only a few geometric variables define all of the meaningful differences between minima and maxima. A few points to note: