rmsdvkpdb.py

#!/usr/bin/python


import sys
sys.path.append("/home/visvaldas/software/openbabel/lib")
import math
import pybel

def squared_distance(coordsA, coordsB):
    """Find the squared distance between two 3-tuples"""
    sqrdist = sum( (a-b)**2 for a, b in zip(coordsA, coordsB) )
    return sqrdist

def rmsd(allcoordsA, allcoordsB):
    """Find the RMSD between two lists of 3-tuples"""
    deviation = sum(squared_distance(atomA, atomB) for
                    (atomA, atomB) in zip(allcoordsA, allcoordsB))
    return math.sqrt(deviation / float(len(allcoordsA)))

# Make sure we have a filename
try:
    filename1 = sys.argv[1]
except:
    print "Usage: python energy.py filename1 filename2"
    sys.exit(1)
try:
    filename2 = sys.argv[2]
except:
    print "Usage: python energy.py filename1 filename2"
    sys.exit(1)


if __name__ == "__main__":
    # Read crystal pose
    crystal = next(pybel.readfile("pdb", filename1))

    # Find automorphisms involving only non-H atoms
    mappings = pybel.ob.vvpairUIntUInt()
    bitvec = pybel.ob.OBBitVec()
    lookup = []
    for i, atom in enumerate(crystal):
        if not atom.OBAtom.IsHydrogen():
            bitvec.SetBitOn(i+1)
            lookup.append(i)
    success = pybel.ob.FindAutomorphisms(crystal.OBMol, mappings, bitvec)

    # Find the RMSD between the crystal pose and each docked pose
    xtalcoords = [atom.coords for atom in crystal if not atom.OBAtom.IsHydrogen()]
    for i, dockedpose in enumerate(pybel.readfile("pdb", filename2)):
        posecoords = [atom.coords for atom in dockedpose if not atom.OBAtom.IsHydrogen()]
        minrmsd = 999999999999
        for mapping in mappings:
            automorph_coords = [None] * len(xtalcoords)
            for x, y in mapping:
                automorph_coords[lookup.index(x)] = xtalcoords[lookup.index(y)]
            mapping_rmsd = rmsd(posecoords, automorph_coords)
            if mapping_rmsd < minrmsd:
                minrmsd = mapping_rmsd
        print("The RMSD is %.3f for pose %d" % (minrmsd, (i+1)))