780.20: 1094 Session 14

Handouts: printouts of multimin_sa_compare.cpp and multimin_sim.cpp, excerpt on simulated annealing

In this session, we'll consider multidimensional minimization with standard downhill routines for local minima and then an application of Monte Carlo techniques to find the global minimum: simulated annealing.

Simulated Annealing

Standard optimization methods are very good in finding local minima near where the minimization was started, but not good at finding the global minimum. Finding the global minimum of a function (such as the energy) is often (but not always) the goal. One strategy using conventional minimizers is to start the minimization at many different places in the parameter space (perhaps chosen at random) and keep the best minimum found.

Here we'll adapt our Monte Carlo algorithm for generating a canonical Boltzmann distribution of configurations at a temperature T to mimic how physical systems find their ground states (i.e., the energy minimum at T=0). At high temperature (compared to characteristic energy spacings), the equilibrium distribution will include many states. If the system is cooled slowly, then it will have the opportunity to settle into the lowest energy state as T goes to zero. This is called annealing. If the system is cooled quickly, it can get stuck in a state that is not the minimum ("quenching"); this is analogous to the routines we looked at in Session 10, which rapidly go "downhill" to find local minima.

The strategy here is to simulate the annealing process by treating the function to be minimized as an energy (it might actually be an energy!), introducing an artificial temperature T, generating a sequence of states in a canonical distribution via the Metropolis algorithm. Then we lower the temperature according to a "schedule" and let the system settle into (hopefully!) a configuration that minimizes the energy.

Global vs. Local Minimum

We first consider the simple example of a damped sine wave in one dimension, which has many local minima.
f(x) = e^-(x-1)2sin(8x)

(This example is from the GSL documentation.) Suppose we start with an initial guess of 15.5; let's see how "downhill" minimization methods compare to simulated annealing.

1. Plot the function in Mathematica so we know what to expect ahead of time (remember that in the general case of interest we won't be able to do this!). Use FindMinimum starting near 15.5 to find the local minimum. Then find the global minimum (using the graph to identify a good starting value). Answers:
   
   
   
2. The code multimin_sa_compare.cpp (the multimin_sa_compare project should also include random_seed.cpp) applies the GSL routines we used in a previous session to find the minimum and then applies a simple simulated annealing algorithm. Run the program and look at the code to figure out where the simulated annealing control parameters (starting and finishing temperatures, maximum step size, cooling rate, # of iterations at each temperature) are set. How well does each method work?
   
   
   
   
3. Now adjust the simulated annealing control parameters to give it a chance! The step size should be large enough (or small enough) so that successes are neither assured nor happen too infrequently (50% would be ideal!). The initial temperature should be considerably larger than the largest typical delta_E normally encountered. You shouldn't cool to rapidly. Make these adjustments and rerun the code. What did you do to get it to work best?
   
   
   
   
4. Just to check that the Metropolis algorithm is doing something, once the code is working, change the line with "random < exp()" to "random > exp()" and verify that it no longer works.

Multidimensional Minimization with GSL Routines II

In Session 10, you tried the gsl routines for multidimensional minimization. Here we finish up that exploration by considering a real physics problem: minimizing the potential energy of molecular configurations built from Na and Cl ions to find stable configurations (this problem is described in a handout from Pang). This will work without complications for small configurations, but bigger clusters will be tricky.

1. Consider how you would modify the multidimensional minimization code from Session 10 to study the problem described in the Pang handout: determining the stable structures of molecules of Na and Cl ions.
   • Read the Pang handout. What is the physics of the interaction energy (i.e., what does each term represent?)? What is the minimization problem to be solved?
     
     
     
     
   • What information to you need to find r_ij for every possible pair of ions? What information do you want to have when the program is finished? (I.e., what should the output be?)
     
     
     
     
   • Does the vector you pass to the minimizer include the x,y,z coordinates of every ion? Write down what is passed for the cases of NaCl, Na₂Cl, and Na₂Cl₂.
     
     
     
     
   • What will the functions "my_f" and "my_df" look like? Remember that they do their calculation based on the vector passed to them, which will not include every x,y,z coordinate (why not?).
     
     
     
     
   • Discuss with your partner a plan to solve the problem using the gsl minimization routines.
2. Now consider the code multimin_nacl.cpp, which I wrote to address this problem. (Warning: This is not a great code.) Create a multimin_nacl project (with random_seed.cpp).
   • Compare your plan to my solution. Explain why your plan is better. :)
     
     
   • In what ways can you check if the problem is coded correctly? List some of the tests you would make. [Hint 1: Think of some things you can change and still expect to get the same answer from the minimization.] [Hint 2: Can you check a simple case a different way?] What symmetry can you use to (partially) check Na₂Cl?
     
     
     
     
   • Are the code's answers for two (1 Na and 1 Cl) and three atoms (2 Na's and 1 Cl) physically reasonable? Here are some questions that address this point: Should the final energy be positive or negative? What is a reasonable number of eV? What should the distance be roughly in angstroms? (E.g., what are typical atomic spacings?) Should the atoms be closer or further apart in NaCl compared to Na₂Cl?
     
     
     
     
   • Now try to find the structure of Na₂Cl₂. Keep in mind that the program identifies local minima. There may be more stable configurations. How would you tell which is more stable?
   • Now consider Na₃Cl₂⁺ and Na₃Cl₃. Can you reproduce the claimed minimum energy structures found in Fig. 4.2? (Note: Make sure the solution has converged to a minimum; you may need to increase the number of iterations allowed, adjust the initial step size, or adjust the tolerance.)
     
     
     
     
   • To check for better local minima, try different starting points (for the x-vector) and try different algorithms (it makes a difference!!!). Finding a global minimum (and being sure you've found it) is a tough problem! To address it, we turn to simulated annealing.

Applying Simulated Annealing to Find Molecular Configurations

Now we'll apply simulated annealing to the problem of the shape of molecules built from Na and Cl atoms. Warning: This method (at least as implemented here) is far from automatic. You have to carefully explore the control parameters.

1. The code "multimin_sim.cpp" is basically the "multimin_nacl.cpp" code with a simulated annealing step inserted between two "downhill" minimizations. The starting point is still the arbitrary set of coordinates from the original code. Look at the code and find where the simulated annealing control parameters are set. Note that the starting values mean that the simulated annealing part is skipped.
   
   
   
   
2. Create a multimin_sim project and run the code. Record the best answers for three Na's and two and three Cl's.
   
   
   
   
3. Now play with the control parameters and try to do better. Check with your neighbors to see who can get the lowest energies. Good luck!