If you reached this page after writing successfully your own version, you are probably ready to tackle level two or even level three assignments. If you reached this page without being able to write a minimizer of your own, examine the program below carefully and consider answering level one or two assignments on the next page.
#!/usr/bin/env python
#
# Steepest Descent minimizer: Harmonic Potential
# Dr. Kalju Kahn, Chem 126 (UCSB), January 2006
# DeepestDescent function partially based on Rick Muller's Minimizers.py
# See http://pyquante.sourceforge.net/ for PyQuante 1.5
#
# Level two: longer, nicer, more flexible
# a) Much more code for doing the same task!
# b) We define functions and output formatting
# c) We do not fix inputs in the program code but read these from files or command line
import sys
# Check if all the input is provided
if len(sys.argv ) !=3 :
sys.exit ( "Usage: HarmOsc_SD2.py param_file initial_r" )
#
# Potential parameters come from the text file (we can use the same program for CO and HF)
#
par_file = open( sys.argv[ 1 ] )
r_ini = float(sys.argv[ 2 ]) # Initial distance from the command line (e.g. 1.4)
parameters = par_file.readlines()
for parameter in parameters:
fields = parameter.split()
mol = (fields [0] )
k = float (fields [1] ) # Harmonic force constant
r_eq = float (fields [2] ) # Potential minimum distance
#
# We now define some frequently used functions
# This makes it easier to read and change the program
#
# Hamiltonian: Harmonic potential
def Hamiltonian(r):
pot = 0.5 * k * (r - r_eq)**2
return pot
# Return the force via analytic derivative
def Force(r):
grad = k * (r - r_eq)
return grad
# Steepest Descent search with forces evaluated at each step
def SteepestDescent(r):
step = 0.0001
E = Hamiltonian(r)
F = Force(r)
print "-----------------------------------------"
print "Step Distance Energy Force "
print "-----------------------------------------"
print " 0 %6.4f %8.3f %8.3f " % (r, E, F)
for i in range(1,50): # We do not allow user more than 50 steps
r -= step*F # Change of distance is proportional to the force
E = Hamiltonian(r)
F = Force(r)
print "%3d %6.4f %8.3f %8.3f " % (i, r, E, F)
if (F < 0.01): # We do not allow user to change convergence criteria
break
return
# Go find the minimum
SteepestDescent(r_ini)