Find_Hydrogen_Bonds.py

import math
import pybel

def sqr_dist(a, b):
    ac = a.coords
    bc = b.coords
    return (ac[0]-bc[0])**2 + (ac[1]-bc[1])**2 + (ac[2]-bc[2])**2

# Definitions taken from
#    http://baoilleach.blogspot.com/2007/07/pybel-hack-that-sd-file.html
# I should fix them to incorporate GMC2007's comments on that page
HBD = pybel.Smarts("[!#6;!H0]")
HBA = pybel.Smarts("[$([$([#8,#16]);!$(*=N~O);" +
                   "!$(*~N=O);X1,X2]),$([#7;v3;" +
                   "!$([nH]);!$(*(-a)-a)])]")

HBOND_CUTOFF = 3.1
SQR_HBOND_CUTOFF = HBOND_CUTOFF**2

if __name__ == "__main__":

    protein = pybel.readfile("pdb", "1HV6_macromol.pdbqt").next()
    ligand = pybel.readfile("pdb", "1hv6_trimer_001_model_.pdb").next()
    protein_atoms = protein.atoms
    ligand_atoms = ligand.atoms

    # Protein donor, ligand acceptor
    for protein_donor in HBD.findall(protein):
        protein_atom = protein_atoms[protein_donor[0] - 1]
        for ligand_acceptor in HBA.findall(ligand):
            ligand_atom = ligand_atoms[ligand_acceptor[0] - 1]
            mysqr_dist = sqr_dist(protein_atom, ligand_atom)
            if mysqr_dist < SQR_HBOND_CUTOFF:
                print "H Bond between protein don (%d) and ligand acc (%d) (%.1f Ang)" % (
                    protein_donor[0], ligand_acceptor[0], math.sqrt(mysqr_dist))

    # Protein acceptor, ligand donor
    for protein_acceptor in HBA.findall(protein):
        protein_atom = protein_atoms[protein_acceptor[0] - 1]
        for ligand_donor in HBD.findall(ligand):
            ligand_atom = ligand_atoms[ligand_donor[0] - 1]
            mysqr_dist = sqr_dist(protein_atom, ligand_atom)
            if mysqr_dist < SQR_HBOND_CUTOFF:
                print "H Bond between protein acc (%d) and ligand don (%d) (%.1f Ang)" % (
                    protein_acceptor[0], ligand_donor[0], math.sqrt(mysqr_dist))