Illustrating energy minimization: a simple python program

Here is a simple python program I wrote to illustrate energy minimization.  The program uses steepest descent and a force field to minimize the energy of a water molecule in internal coordinates.

If you have a Mac or Linux machine you already have python installed (Windows users please Google and/or leave a comment).  You can also get an invite for koding.com and run it there. You can run it by opening a terminal and typing "python Emin.py".

Things to play around with (once you have it running): change the starting geometry, step size (c), number of steps (n_steps), try printing out the geometry, energy and gradient for each step.

You can also try to write a similar program for other molecules (like methane), although once you have van der Waals interactions you will probably have to switch to Cartesian coordinates and things get complicated.

	# A simple energy minimization program that uses steepest descent
	# and a force field to minimize the energy of water in internal coordinates
	# Written by Jan H. Jensen, 2013, released in the the GNU GPL license
	kOH = 50.0
	rOHe = 0.95
	kHOH = 50.0
	thetaHOHe = 104.5
	c = 0.005
	n_steps = 20
	#starting geometry
	rOH = 10.0
	thetaHOH = 180.0
	def Eandg(rOH,thetaHOH):
	E = 2*kOH*(rOH-rOHe)**2 + kHOH*(thetaHOH-thetaHOHe)**2
	grOH = 2*kOH*(rOH-rOHe)
	gthetaHOH = 2*kHOH*(thetaHOH-thetaHOHe)
	return (E,grOH,gthetaHOH)
	for i in range(n_steps):
	(E,grOH,gthetaHOH) = Eandg(rOH,thetaHOH)
	if (abs(grOH) >0.001/c or abs(gthetaHOH) > 0.01/c ):
	rOH = rOH - c*grOH
	thetaHOH = thetaHOH - c*gthetaHOH
	if (abs(grOH) >0.001/c or abs(gthetaHOH) > 0.01/c ):
	print "not converged"
	else:
	print "converged"
	print E,rOH,thetaHOH

view raw Emin.py hosted with ❤ by GitHub

Related blog posts
The force is strong in this one
The autoopt tool in Avogadro