A python program for reading coordinates and gradients from a GAMESS output file

Some people liked the previous post on python, so I dusted off this python program for reading locating coordinates and gradients in a GAMESS log file and writing out and .xyz+vib file (to create this post).  I have only used it on this output file, but it should work for other files as well.

Even if you have no interest in .xyz+vib files, the program might be useful as an example of how to locate specific information in a file and read it in. If you have similar programs (or questions) please feel free to share them in the comments

	"""\
	This python program extracts coordinates and gradients from a GAMESS output
	file and creates a .xzy+vib file. Note that the sign of the gradient is chaged to produce the force
	To use it type (assuming opt.py and the output file is in the same directory):
	python opt.py x xx > xxx.xyz
	where x is the number of atoms in the molecule and xx is the name of the GAMESS output file. xxx is whatever you want to call your xyz file.
	Written by Jan H. Jensen, released under the GNU GPL license
	"""
	import string,sys
	atoms = eval(sys.argv[1]) # convert x to number
	filename = sys.argv[2]
	file = open(filename,'r')
	coord = []
	grad = []
	while 1:
	line = file.readline() # read a line from he file
	if not line: break # if end-of-file: quit
	words = string.split(line) # split line into words
	# equilibrium coordinates are printed out a second time and should not be read
	if len(words) > 1 and words[1] == 'EQUILIBRIUM': break
	# find and read coordniates
	if len(words) > 0 and words[0] == 'COORDINATES': # find line with EQUILIRBIUM as the 1st word
	file.readline() # skip a line
	file.readline() # skip a line
	for i in range(atoms): # read in x, y, z coordinates
	line = file.readline()
	words = string.split(line)
	coord.append((words[0],words[2],words[3],words[4]))
	# find and read gradient
	if len(words) > 0 and words[0] == 'NSERCH=': # find line with NSERCH= as the 1st word
	for i in range(6):
	file.readline() # skip 6 lines
	for i in range(atoms): # read x, y, z component of gradient
	line = file.readline()
	words = string.split(line)
	grad.append((eval(words[3]),eval(words[4]),eval(words[5]))) # gradients to be multiplied by -1, so convert to number
	# print in .xyz+vib format and change sign of gradient
	scale = -1.
	for i in range( len(coord) ):
	if i%atoms == 0: print atoms
	if i%atoms == 0: print 'title'
	tog = (coord[i][0],coord[i][1],coord[i][2],coord[i][3],scale*grad[i][0],scale*grad[i][1],scale*grad[i][2])
	print '%s %s %s %s %f %f %f' % tog

view raw opt.py hosted with ❤ by GitHub

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

MOPAC/PCM interface in GAMESS

A new version of GAMESS was released recently.  Among the many new features is +Casper Steinmann's recent work on interfacing the PCM solvation method with the semiempirical methods AM1 and PM3.  So you can now simulate molecules in solution with AM1 and PM3.  

Casper also parallelized AM1 and PM3 (both in gas phase and solution), though molecules have to be pretty big before you see an appreciable speed up, as you can see from the figure below (on small proteins).

The work has been published in the open access journal PLoS ONE. +Jimmy Charnley Kromann is hard at work on implementing PM6.

This work is licensed under a Creative Commons Attribution 3.0 

Computational Chemistry Highlights: June issue

The June issue of Computational Chemistry Highlights is out.

CCH is an overlay journal that identifies the most important papers in computational and theoretical chemistry published in the last 1-2 years. CCH is not affiliated with any publisher: it is a free resource run by scientists for scientists. You can read more about it here.

Table of content for this issue features contributions from CCH editors Steven Bachrach, Jonathan Goodman, and Jan Jensen:

Molecularspace.org

Three interesting recent Angew. Chem. papers

Unraveling the Enigmatic Mechanism of L‐Asparaginase II with QM/QM calculations

Interested in more? There are many ways to subscribe to CCH updates.