Previous discussion suggest that minimization typically involves repetitive evaluation of values and derivatives of possibly rather complex multi-variable functions. Such task are well suited for computers. Writing an efficient minimizer for problems commonly encountered in quantum chemical calculations can be considered a challenging programming task even for expert programmers. In this tutorial we will illustrate how to program minimizers for a simple quadratic function without thinking much about the memory efficiency or the execution speed of the code.
It is relatively straightforward to write a computer program that minimizes one-dimensional functions using the steepest descent or Newton-Rhapson methods. In general, such a program must have the following parts:
The code below illustrates the Python implementation of the steepest descent minimizer for the one-dimensional harmonic potential. While not as fast as C, FORTRAN or assembler languages, Python is easy to learn and allow to write very human-readable code. Furthermore, because Python interpreters are freely available for most modern computer platforms, the code can be trivially ported from one system to another. Please note that on Linux systems, the first line is needed to locate the Python interpreter.
#!/usr/bin/env python
# The previous line is required on Linux systems to find the Python interpreter
#
# Steepest Descent minimizer: Harmonic Potential
# Dr. Kalju Kahn, Chem 126 (UCSB), January 2006
#
# Level one: short, ugly, and rigid
# Potential Function Parameters
k = 2743.0 # Harmonic force constant
r_eq = 1.1283 # Equilibrium distance
# Initial Search Point
r_ini = 1.55
# Steepest Descent search
r = r_ini
step = 0.0001
E = 0.5 * k * (r - r_eq)**2 # Energy
F = k * (r - r_eq) # Force (first derivative)
print "Step Distance Energy Force Hessian "
print "0", r, E, F # Initial values
for i in range(1,50): # Iterate up to 50 times
r = r - step*F # Steepest Descent update
E = 0.5 * k * (r - r_eq)**2
F = k * (r - r_eq)
print i, r, E, F
if (abs(F) < 0.01): # Are we at the stationary point already?
break
You should copy this code to a new text file, and save it as HarmOsc_SD1.py. If under Unix, set the permission of this file to 755 (executable) and run the program by typing on the Unix shell
chmod 755 HarmOsc_SD1.py ./HarmOsc_SD1.py
Examine the output and verify that the minimum (1.1283 Angstroms, of course) has been found.
The code above illustrates well the core ideas of the steepest descent minimizer. Implementation in other languages would, in general, follow a similar logic but the syntax would be different. For example, the steepest descent minimizer can be implemented in Mathematica© as shown below. Notice that this program takes advantage of Mathematica©'s ability to find of functions. This feature can come in handy with more complicated functions that would be tedious to differentiate by hand or error-prone to code into a program.
Remove["Global`*"]
Off [General::spell]
(* Steepest Descent minimizer: Harmonic Potential;
Dr. Kalju Kahn,Chem 126 (UCSB),January 2006
Level one: short, ugly, and rigid;
*)
(* Potential Function Parameters *)
k = 2743.0 ;
Req = 1.1283;
(* Initial Search Point *)
IniR = 1.55;
(* Steepest Descent search: ER and FR are energy and force at point R *)
Energy = 0.5 k (R - Req)^2;
Step = 0.0001;
ER = Energy /. {R -> IniR};
FR = D[Energy, R] /. {R -> IniR};
Print[" R ", " Energy " , "Force"]
Print[IniR, " ", ER, " ", FR]
NextR = IniR;
For[i=1, i<50, i=i+1,
NextR = NextR - Step*FR;
ER = Energy /. {R -> NextR};
FR = D[Energy, R] /. {R -> NextR};
Print[NextR," ", ER " ",FR]
If [(Abs[FR] < 0.01), Break[] ]
]