Calculates centroids and points above and below a series of rings in a molecule. Useful for NICS calculations.

rings2nics.py

#!/usr/bin/env python

import sys
from pybel import readfile
import numpy as np

if __name__ == "__main__":
	if len(sys.argv) < 2:
		sys.exit("Expecting one argument. Abort.")

	for molecule in readfile("xyz", sys.argv[1]):
		while True:
			typed = raw_input("Atom #s: ").strip()
			if typed == "exit":
				print "Bye!"
				break

			list_of_atoms = [int(n) for n in typed.split()]

			atoms = []
			for i in list_of_atoms:
				atoms.append(np.array(molecule.atoms[i - 1].coords))

			n = np.cross(atoms[2] - atoms[0], atoms[1] - atoms[0])
			n = n/np.linalg.norm(n)

			centr = sum(atoms)/len(atoms)
			upper = centr + n
			lower = centr - n

			print "X	%f	%f	%f" % (upper[0], upper[1], upper[2])
			print "X	%f	%f	%f" % (centr[0], centr[1], centr[2])
			print "X	%f	%f	%f" % (lower[0], lower[1], lower[2])