780.20: 2082 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:
Multidimensional Minimization with GSL Routines II
In Session 10, you tried the gsl routines for
Here we finish up that exploration, then consider a real physics problem:
minimizing the potential energy of molecular configurations built from
Na and Cl ions to find stable configurations.
This will work without complications for small configurations, but
bigger clusters will be tricky.
- [Move on to the next part if you completed this
in Session 10.] Returning to "multimin_test.cpp",
what algorithm is used initially to do the minimization?
Modify the code to try some other algorithms (some major changes
would be necessary to use the simplex algorithm; why?).
Which is "best" for this problem?
- Consider how you would modify the code 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
- What information to you need to find rij 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, Na2Cl, and Na2Cl2.
- 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?).
- Sketch out a plan to solve the
problem using the gsl minimization routines.
If at all possible, devise your plan in collaboration
with a classmate.
- Now consider the code multimin_nacl.cpp, which I wrote
to address this problem.
(Warning: This is not a great code.)
Compile and link with make_multimin_nacl.
- 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
What symmetry can you use to (partially) check Na2Cl?
- Are the 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
- Now try to find the structure of Na2Cl2.
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 Na3Cl2+
and Na3Cl3. 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
- 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! (We'll return to this problem later when
discussing 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 21, 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
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; how do "downhill" minimization methods compare to
- 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!). Using FindMinimum starting near 15.5 to find the local
minimum. Then find the global minimum.
- The code multimin_sa_compare.cpp (compile and link with
make_multimin_sa_compare) 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?
- 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. Can you get it to work?
- 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.
Applying Simulated Annealing to Find Molecular Configurations
Now we'll apply essentially the same approach to the problem of the
shape of molecules built from Na and Cl atoms.
- The code "multimin_sim.cpp" is basically the "multimin_nacl.cpp"
code from an earlier session 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.
What Metropolis algorithm is used
in this code compared to the "multimin_sa_compare.cpp" code?
- Compile, link, and run the code (use "make_multimin_sim") and
record the best answers for three Na's and two and three Cl's.
- Now play with the control parameters and try to do better.
780.20: 2082 Session 14.
Last modified: 07:16 pm, February 21, 2005.