reduce number of duplicate hbond interactions when ignoring hydrogens

Acellera · Jan 10, 2024 · e772289 · e772289
1 parent 46515d4
commit e772289
Show file tree

Hide file tree

Showing 2 changed files with 13 additions and 5 deletions.
diff --git a/moleculekit/interactions/hbonds/hbonds.pyx b/moleculekit/interactions/hbonds/hbonds.pyx
@@ -57,7 +57,10 @@ def calculate(
             for a in range(n_acceptors):
                 a_idx = acceptors[a]
                 d_idx_d = donors[d, 0]
-                d_idx_h = donors[d, 1]
+                d_idx = d_idx_d
+                if not ignore_hs:
+                    d_idx_h = donors[d, 1]
+                    d_idx = d_idx_h
 
                 # Don't calculate h-bonds with identical donor-acceptor atom
                 if a_idx == d_idx_d:
@@ -74,9 +77,6 @@ def calculate(
 
                 dist2_a = 0
                 dist2_b = 0
-                d_idx = d_idx_h
-                if ignore_hs: # If we ignore the hydrogen just calculate heavy to heavy distance
-                    d_idx = d_idx_d
                 # Calculate donor hydrogen (or heavy if ignore_hs) to acceptor vector
                 for i in range(3):
                     val = coords[a_idx, i, f] - coords[d_idx, i, f] 

diff --git a/moleculekit/interactions/interactions.py b/moleculekit/interactions/interactions.py
@@ -6,7 +6,7 @@
 import networkx as nx
 import numpy as np
 import unittest
-from moleculekit.tools.moleculechecks import isProteinProtonated
+from moleculekit.tools.moleculechecks import isProteinProtonated, proteinHasBonds
 import logging
 
 logger = logging.getLogger(__name__)
@@ -95,6 +95,10 @@ def get_donors_acceptors(mol, exclude_water=True, exclude_backbone=False):
         raise RuntimeError(
             "The protein seems to not be protonated. You must provide a protonated system for H-bonds to be detected."
         )
+    if not proteinHasBonds(mol):
+        raise RuntimeError(
+            "The protein is missing bonds. Cannot calculate donors/acceptors"
+        )
 
     donor_elements = ("N", "O")
     acceptor_elements = ("N", "O")
@@ -366,6 +370,10 @@ def hbonds_calculate(
     donors = donors[np.all(np.isin(donors, sel_idx), axis=1)]
     acceptors = acceptors[np.isin(acceptors, sel_idx)]
 
+    if ignore_hs:
+        # if we are ignoring hydrogens, reduce the donors list to only heavy atoms
+        donors = np.unique(donors[:, 0])[:, None]
+
     hb = hbonds.calculate(
         donors,
         acceptors,