Skip to content

Instantly share code, notes, and snippets.

@ghutchis
Last active April 20, 2024 16:30
Show Gist options
  • Save ghutchis/0980e8faf51be56780af469ccb0e3311 to your computer and use it in GitHub Desktop.
Save ghutchis/0980e8faf51be56780af469ccb0e3311 to your computer and use it in GitHub Desktop.
Open Babel / Pybel script to report all bonds, angles, torsions, and non-bonded interactions
#!/usr/bin/env python
import sys, os
from openbabel import openbabel as ob
from openbabel import pybel
# syntax:
# molml.py [files]
def atomType(mol, atomIdx):
# get the atomic type given an atom index
return mol.OBMol.GetAtom(atomIdx).GetType()
# repeat through all the files on the command-line
# we can change this to use the glob module as well
# e.g., find all the files in a set of folders
for argument in sys.argv[1:]:
filename, extension = os.path.splitext(argument)
# read the molecule from the supplied file
mol = next(pybel.readfile(extension[1:], argument))
print(mol.energy) # in kcal/mol
# ideally, we should turn this into an atomization energy
# iterate through all atoms
# .. this is commented out because Bag Of Bonds doesn't use atomic charges
# for atom in mol.atoms:
# print "Atom %d, %8.4f" % (atom.type, atom.partialcharge)
# iterate through all bonds
bonds = []
for bond in ob.OBMolBondIter(mol.OBMol):
begin = atomType(mol, bond.GetBeginAtomIdx())
end = atomType(mol, bond.GetEndAtomIdx())
if (end < begin):
# swap them for lexographic order
begin, end = end, begin
bonds.append("Bond %s-%s, %8.4f" % (begin, end, bond.GetLength()) )
print(bonds[-1])
# iterate through all angles
angles = []
for angle in ob.OBMolAngleIter(mol.OBMol):
a = (angle[0] + 1)
b = mol.OBMol.GetAtom(angle[1] + 1)
c = (angle[2] + 1)
aType = atomType(mol, a)
cType = atomType(mol, c)
if (cType < aType):
# swap them for lexographic order
aType, cType = cType, aType
angles.append( "Angle %s-%s-%s, %8.3f" % (aType, b.GetType(), cType, b.GetAngle(a, c)) )
print(angles[-1])
# iterate through all torsions
torsions = []
for torsion in ob.OBMolTorsionIter(mol.OBMol):
a = (torsion[0] + 1)
b = (torsion[1] + 1)
c = (torsion[2] + 1)
d = (torsion[3] + 1)
aType = atomType(mol, a)
bType = atomType(mol, b)
cType = atomType(mol, c)
dType = atomType(mol, d)
# output in lexographic order
if (aType < dType):
torsions.append( "Torsion %s-%s-%s-%s, %8.3f" % (aType, bType, cType, dType, mol.OBMol.GetTorsion(a, b, c, d)) )
else:
torsions.append( "Torsion %s-%s-%s-%s, %8.3f" % (dType, cType, bType, aType, mol.OBMol.GetTorsion(a, b, c, d)) )
print(torsions[-1])
nb = []
for pair in ob.OBMolPairIter(mol.OBMol):
(first, second) = pair
begin = atomType(mol, first)
end = atomType(mol, second)
if (end < begin):
# swap them for lexographic order
begin, end = end, begin
dist = mol.OBMol.GetAtom(first).GetDistance(second)
nb.append("NB %s-%s, %8.4f" % (begin, end, dist))
print(nb[-1])
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment