hbonds_like_molstar.py

prot_dons = np.array(pocket.GetSubstructMatches(HDonorSmarts)).flatten()
prot_accs = np.array(pocket.GetSubstructMatches(HAcceptorSmarts)).flatten()


lig_dons = np.array(ligand.GetSubstructMatches(HDonorSmarts)).flatten()
lig_accs = np.array(ligand.GetSubstructMatches(HAcceptorSmarts)).flatten()

def unit_vector(vector):
    """ Returns the unit vector of the vector.  """
    return vector / np.linalg.norm(vector)

def angle_between(v1, v2):
    """ Returns the angle in radians between vectors 'v1' and 'v2'::

            >>> angle_between((1, 0, 0), (0, 1, 0))
            1.5707963267948966
            >>> angle_between((1, 0, 0), (1, 0, 0))
            0.0
            >>> angle_between((1, 0, 0), (-1, 0, 0))
            3.141592653589793
    """
    v1_u = unit_vector(v1)
    v2_u = unit_vector(v2)
    return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0)) * (180/np.pi)
  
def hbond_checker(pos1, pos2, atom1, atom2):
    """
    Atom1 is donor (rdkit atom)
    Atom2 is acceptor (rdkit atom)
    pos is np.array of [natoms, 3] coordinates.
    """
    
    ideal_angles = {Chem.HybridizationType.SP3: 109.4721,
            Chem.HybridizationType.SP2: 120,
            Chem.HybridizationType.SP: 180}
    #octahedral would be 90.
    
    atom1_xyz = pos1[atom1.GetIdx()] #donor
    atom2_xyz = pos2[atom2.GetIdx()] #acceptor
    
    #print(atom1.GetSymbol(), atom1.GetHybridization(), 'at1:', atom1.GetIdx(), 'at2:', atom2.GetIdx())

    donor_neighbours = [at for at in atom1.GetNeighbors() if at.GetSymbol()!='H']

    ##
    ##First, sort out whether the acceptor atom lies within an acceptable
    ##area of a hydrogen of the donor.
    ##
    angles = []
    for neighb in donor_neighbours:
        angle = angle_between(
            atom1_xyz-atom2_xyz,
            atom1_xyz-pos1[neighb.GetIdx()]
        )
        angles.append(angle)
        
    ideal_angle = ideal_angles[atom1.GetHybridization()]

    if any(abs(ideal_angle-angle)>45 for angle in angles):
        return False
    
    ##
    ##Second, we have an additional cut-off for being out-of-plane
    ##for trigonal planar atoms. 
    ##
    if atom1.GetHybridization()==Chem.HybridizationType.SP2:
        #must have two neighbours:
        if len(donor_neighbours)<2:
            donor_neighbours.append(
                donor_neighbours[0].GetNeighbors()[0]
            )
        p2, p3 = [pos1[a.GetIdx()] for a in donor_neighbours]

        ##finding the normal vector to the plane
        # Get the difference between p1 and p2, and p1 and p3
        p1p2 = p2 - atom1_xyz
        p1p3 = p3 - atom1_xyz
        # Use numpy.cross to get the cross product of p1p2 and p1p3
        normal = np.cross(p1p2, p1p3)
        normal = normal/np.linalg.norm(normal)
        angle = 90-angle_between(atom1_xyz-atom2_xyz, normal)
        #print('tautomer angle:', angle)
        if angle>45:
            return False
        

    acceptor_neighbours = [at for at in atom2.GetNeighbors() if at.GetSymbol()!='H']


    angles = []
    for neighb in acceptor_neighbours:
        angle = angle_between(
            atom2_xyz-atom1_xyz,
            atom2_xyz-pos2[neighb.GetIdx()]
        )

          
        angles.append(angle)
        
    ideal_angle = ideal_angles[atom2.GetHybridization()]

    if any([(ideal_angle-angle)>45 for angle in angles]):

        return False
    
    ##
    ##Second, we have an additional cut-off for being out-of-plane
    ##for trigonal planar atoms. 
    ##
    if atom2.GetHybridization()==Chem.HybridizationType.SP2:
        #must have two neighbours:
        if len(acceptor_neighbours)<2:
            acceptor_neighbours.append(
                acceptor_neighbours[0].GetNeighbors()[0]
            )
        p2, p3 = [pos2[a.GetIdx()] for a in acceptor_neighbours]

        ##finding the normal vector to the plane
        # Get the difference between p1 and p2, and p1 and p3
        p1p2 = p2 - atom2_xyz
        p1p3 = p3 - atom2_xyz
        # Use numpy.cross to get the cross product of p1p2 and p1p3
        normal = np.cross(p1p2, p1p3)
        normal = normal/np.linalg.norm(normal)
        angle = 90-angle_between(atom2_xyz-atom1_xyz, normal)
        
        if angle>90:
            return False
        
    return True

prot_atoms = np.array([at for at in pocket.GetAtoms()])
lig_atoms = np.array([at for at in ligand.GetAtoms()])

def compareDonorsAcceptors(pos1, pos2, atoms1, atoms2, dons, accs):
    
    hbonds = np.zeros([pos1.shape[0], pos2.shape[0]]).astype(bool)

    dist = cdist(pos1[dons], pos2[accs])
    m = (dist>1) & (dist<3.5)

    don_map = np.arange(pos1.shape[0])[dons]
    acc_map = np.arange(pos2.shape[0])[accs]
    
    a,b = np.where(m)
    for xidx, yidx in zip(a,b):
        x = don_map[xidx]
        y = acc_map[yidx]
        #print(x, y)
        res = hbond_checker(pos1, pos2, atoms1[x], atoms2[y])
        if res:
            hbonds[x,y]=True
    return hbonds
  
  
intra_hbonds = compareDonorsAcceptors(protpos, 
                                      protpos, 
                                      prot_atoms, 
                                      prot_atoms, 
                                      prot_dons, 
                                      prot_accs)

inter_hbonds = compareDonorsAcceptors(protpos, 
                                      ligpos, 
                                      prot_atoms, 
                                      lig_atoms, 
                                      prot_dons, 
                                      lig_accs)

inter_hbonds += compareDonorsAcceptors(ligpos, 
                                      protpos, 
                                      lig_atoms, 
                                      prot_atoms, 
                                      lig_dons, 
                                      prot_accs).T



view = viewProt(Chem.RemoveHs(pocket),color='green')
view = viewProt(Chem.RemoveHs(ligand), color='orange', view=view)

for x,y in zip(*np.where(intra_hbonds)):
    xyz1 = protpos[x]
    xyz2 = protpos[y]
    
    view.addCylinder({'color':'red',
                        'start':{'x':xyz1[0],'y':xyz1[1],'z':xyz1[2]},
                        'end':{'x':xyz2[0],'y':xyz2[1],'z':xyz2[2]},
                    'radius':0.1, 'dashed':True});
    
for x,y in zip(*np.where(inter_hbonds)):
    xyz1 = protpos[x]
    xyz2 = ligpos[y]
    
    view.addCylinder({'color':'red',
                        'start':{'x':xyz1[0],'y':xyz1[1],'z':xyz1[2]},
                        'end':{'x':xyz2[0],'y':xyz2[1],'z':xyz2[2]},
                    'radius':0.1, 'dashed':True});