1"""Module calculates interactions between two molecules
2(proein-protein, protein-ligand, small-small).
3Currently following interacions are implemented:
4    * hydrogen bonds
5    * halogen bonds
6    * pi stacking (parallel and perpendicular)
7    * salt bridges
8    * hydrophobic contacts
9    * pi-cation
10    * metal coordination
11    * pi-metal
12"""
13
14import numpy as np
15from oddt.spatial import angle, angle_2v, distance
16
17__all__ = ['close_contacts',
18           'hbond_acceptor_donor',
19           'hbonds',
20           'halogenbond_acceptor_halogen',
21           'halogenbonds',
22           'pi_stacking',
23           'salt_bridge_plus_minus',
24           'salt_bridges',
25           'hydrophobic_contacts',
26           'pi_cation',
27           'acceptor_metal',
28           'pi_metal']
29
30
31def close_contacts(x, y, cutoff, x_column='coords', y_column='coords',
32                   cutoff_low=0.):
33    """Returns pairs of atoms which are within close contac distance cutoff.
34    The cutoff is semi-inclusive, i.e (cutoff_low, cutoff].
35
36    Parameters
37    ----------
38    x, y : atom_dict-type numpy array
39        Atom dictionaries generated by oddt.toolkit.Molecule objects.
40
41    cutoff : float
42        Cutoff distance for close contacts
43
44    x_column, ycolumn : string, (default='coords')
45        Column containing coordinates of atoms (or pseudo-atoms,
46        i.e. ring centroids)
47
48    cutoff_low : float (default=0.)
49        Lower bound of contacts to find (exclusive). Zero by default.
50        .. versionadded:: 0.6
51
52    Returns
53    -------
54    x_, y_ : atom_dict-type numpy array
55        Aligned pairs of atoms in close contact for further processing.
56    """
57    if len(x[x_column]) > 0 and len(x[x_column]) > 0:
58        d = distance(x[x_column], y[y_column])
59        index = np.argwhere((d > cutoff_low) & (d <= cutoff))
60        return x[index[:, 0]], y[index[:, 1]]
61    else:
62        return x[[]], y[[]]
63
64
65def hbond_acceptor_donor(mol1, mol2, cutoff=3.5, base_angle=120, tolerance=30):
66    """Returns pairs of acceptor-donor atoms, which meet H-bond criteria
67
68    Parameters
69    ----------
70    mol1, mol2 : oddt.toolkit.Molecule object
71        Molecules to compute H-bond acceptor and H-bond donor pairs
72
73    cutoff : float, (default=3.5)
74        Distance cutoff for A-D pairs
75
76    base_angle : int, (default=120)
77        Base angle determining allowed direction of hydrogen bond formation,
78        which is devided by the number of neighbors of acceptor atom
79        to establish final directional angle
80
81    tolerance : int, (default=30)
82        Range (+/- tolerance) from perfect direction (base_angle/n_neighbors)
83        in which H-bonds are considered as strict.
84
85    Returns
86    -------
87    a, d : atom_dict-type numpy array
88        Aligned arrays of atoms forming H-bond, firstly acceptors,
89        secondly donors.
90
91    strict : numpy array, dtype=bool
92        Boolean array align with atom pairs, informing whether atoms
93        form 'strict' H-bond (pass all angular cutoffs). If false,
94        only distance cutoff is met, therefore the bond is 'crude'.
95    """
96    a, d = close_contacts(mol1.atom_dict[mol1.atom_dict['isacceptor']],
97                          mol2.atom_dict[mol2.atom_dict['isdonor']],
98                          cutoff)
99    # skip empty values
100    if len(a) > 0 and len(d) > 0:
101        angle1 = angle(d['coords'][:, np.newaxis, :],
102                       a['coords'][:, np.newaxis, :],
103                       a['neighbors'])
104        a_neighbors_num = np.sum(~np.isnan(a['neighbors'][:, :, 0]), axis=-1)[:, np.newaxis]
105        angle2 = angle(a['coords'][:, np.newaxis, :],
106                       d['coords'][:, np.newaxis, :],
107                       d['neighbors'])
108        d_neighbors_num = np.sum(~np.isnan(d['neighbors'][:, :, 0]), axis=-1)[:, np.newaxis]
109        strict = (((np.nan_to_num(angle1) > (base_angle / a_neighbors_num - tolerance)) |
110                   np.isnan(angle1)) &
111                  ((np.nan_to_num(angle2) > (base_angle / d_neighbors_num - tolerance)) |
112                   np.isnan(angle2))).all(axis=-1)
113        return a, d, strict
114    else:
115        return a, d, np.array([], dtype=bool)
116
117
118def hbonds(mol1, mol2, *args, **kwargs):
119    """Calculates H-bonds between molecules
120
121    Parameters
122    ----------
123    mol1, mol2 : oddt.toolkit.Molecule object
124        Molecules to compute H-bond acceptor and H-bond donor pairs
125
126    cutoff : float, (default=3.5)
127        Distance cutoff for A-D pairs
128
129    base_angle : int, (default=120)
130        Base angle determining allowed direction of hydrogen bond formation,
131        which is devided by the number of neighbors of acceptor atom
132        to establish final directional angle
133
134    tolerance : int, (default=30)
135        Range (+/- tolerance) from perfect direction (base_angle/n_neighbors)
136        in which H-bonds are considered as strict.
137
138    Returns
139    -------
140    mol1_atoms, mol2_atoms : atom_dict-type numpy array
141        Aligned arrays of atoms forming H-bond
142
143    strict : numpy array, dtype=bool
144        Boolean array align with atom pairs, informing whether atoms
145        form 'strict' H-bond (pass all angular cutoffs). If false,
146        only distance cutoff is met, therefore the bond is 'crude'.
147    """
148    a1, d1, s1 = hbond_acceptor_donor(mol1, mol2, *args, **kwargs)
149    a2, d2, s2 = hbond_acceptor_donor(mol2, mol1, *args, **kwargs)
150    return np.concatenate((a1, d2)), np.concatenate((d1, a2)), np.concatenate((s1, s2))
151
152
153def halogenbond_acceptor_halogen(mol1,
154                                 mol2,
155                                 base_angle_acceptor=120,
156                                 base_angle_halogen=180,
157                                 tolerance=30,
158                                 cutoff=4):
159    """Returns pairs of acceptor-halogen atoms, which meet halogen bond criteria
160
161    Parameters
162    ----------
163    mol1, mol2 : oddt.toolkit.Molecule object
164        Molecules to compute halogen bond acceptor and halogen pairs
165
166    cutoff : float, (default=4)
167        Distance cutoff for A-H pairs
168
169    base_angle_acceptor : int, (default=120)
170        Base angle determining allowed direction of halogen bond formation,
171        which is devided by the number of neighbors of acceptor atom
172        to establish final directional angle
173
174    base_angle_halogen : int (default=180)
175        Ideal base angle between halogen bond and halogen-neighbor bond
176
177    tolerance : int, (default=30)
178        Range (+/- tolerance) from perfect direction (base_angle/n_neighbors)
179        in which halogen bonds are considered as strict.
180
181    Returns
182    -------
183    a, h : atom_dict-type numpy array
184        Aligned arrays of atoms forming halogen bond, firstly acceptors,
185        secondly halogens
186
187    strict : numpy array, dtype=bool
188        Boolean array align with atom pairs, informing whether atoms
189        form 'strict' halogen bond (pass all angular cutoffs). If false,
190        only distance cutoff is met, therefore the bond is 'crude'.
191    """
192    a, h = close_contacts(mol1.atom_dict[mol1.atom_dict['isacceptor']],
193                          mol2.atom_dict[mol2.atom_dict['ishalogen']],
194                          cutoff)
195    # skip empty values
196    if len(a) > 0 and len(h) > 0:
197        angle1 = angle(h['coords'][:, np.newaxis, :],
198                       a['coords'][:, np.newaxis, :],
199                       a['neighbors'])
200        angle2 = angle(a['coords'][:, np.newaxis, :],
201                       h['coords'][:, np.newaxis, :],
202                       h['neighbors'])
203        a_neighbors_num = np.sum(~np.isnan(a['neighbors'][:, :, 0]), axis=-1)[:, np.newaxis]
204        h_neighbors_num = np.sum(~np.isnan(h['neighbors'][:, :, 0]), axis=-1)[:, np.newaxis]
205        strict = (((np.nan_to_num(angle1) > (base_angle_acceptor / a_neighbors_num - tolerance)) |
206                   np.isnan(angle1)) &
207                  ((np.nan_to_num(angle2) > (base_angle_halogen / h_neighbors_num - tolerance)) |
208                   np.isnan(angle2))).all(axis=-1)
209        return a, h, strict
210    else:
211        return a, h, np.array([], dtype=bool)
212
213
214def halogenbonds(mol1, mol2, **kwargs):
215    """Calculates halogen bonds between molecules
216
217    Parameters
218    ----------
219    mol1, mol2 : oddt.toolkit.Molecule object
220        Molecules to compute halogen bond acceptor and halogen pairs
221
222    cutoff : float, (default=4)
223        Distance cutoff for A-H pairs
224
225    base_angle_acceptor : int, (default=120)
226        Base angle determining allowed direction of halogen bond formation,
227        which is devided by the number of neighbors of acceptor atom
228        to establish final directional angle
229
230    base_angle_halogen : int (default=180)
231        Ideal base angle between halogen bond and halogen-neighbor bond
232
233    tolerance : int, (default=30)
234        Range (+/- tolerance) from perfect direction (base_angle/n_neighbors)
235        in which halogen bonds are considered as strict.
236
237    Returns
238    -------
239    mol1_atoms, mol2_atoms : atom_dict-type numpy array
240        Aligned arrays of atoms forming halogen bond
241
242    strict : numpy array, dtype=bool
243        Boolean array align with atom pairs, informing whether atoms
244        form 'strict' halogen bond (pass all angular cutoffs). If false,
245        only distance cutoff is met, therefore the bond is 'crude'.
246    """
247    a1, h1, s1 = halogenbond_acceptor_halogen(mol1, mol2, **kwargs)
248    a2, h2, s2 = halogenbond_acceptor_halogen(mol2, mol1, **kwargs)
249    return np.concatenate((a1, h2)), np.concatenate((h1, a2)), np.concatenate((s1, s2))
250
251
252def pi_stacking(mol1, mol2, cutoff=5, tolerance=30):
253    """Returns pairs of rings, which meet pi stacking criteria
254
255    Parameters
256    ----------
257    mol1, mol2 : oddt.toolkit.Molecule object
258        Molecules to compute ring pairs
259
260    cutoff : float, (default=5)
261        Distance cutoff for Pi-stacking pairs
262
263    tolerance : int, (default=30)
264        Range (+/- tolerance) from perfect direction (parallel or
265        perpendicular) in which pi-stackings are considered as strict.
266
267    Returns
268    -------
269    r1, r2 : ring_dict-type numpy array
270        Aligned arrays of rings forming pi-stacking
271
272    strict_parallel : numpy array, dtype=bool
273        Boolean array align with ring pairs, informing whether rings
274        form 'strict' parallel pi-stacking. If false, only distance cutoff is met,
275        therefore the stacking is 'crude'.
276
277    strict_perpendicular : numpy array, dtype=bool
278        Boolean array align with ring pairs, informing whether rings
279        form 'strict' perpendicular pi-stacking (T-shaped, T-face, etc.).
280        If false, only distance cutoff is met, therefore the stacking is 'crude'.
281    """
282    r1, r2 = close_contacts(mol1.ring_dict,
283                            mol2.ring_dict,
284                            cutoff,
285                            x_column='centroid',
286                            y_column='centroid')
287    if len(r1) > 0 and len(r2) > 0:
288        angle1 = angle_2v(r1['vector'], r2['vector'])
289        angle2 = angle(r1['vector'] + r1['centroid'],
290                       r1['centroid'],
291                       r2['centroid'])
292        strict_parallel = (((angle1 > 180 - tolerance) | (angle1 < tolerance)) &
293                           ((angle2 > 180 - tolerance) | (angle2 < tolerance)))
294        strict_perpendicular = (((angle1 > 90 - tolerance) & (angle1 < 90 + tolerance)) &
295                                ((angle2 > 180 - tolerance) | (angle2 < tolerance)))
296        return r1, r2, strict_parallel, strict_perpendicular
297    else:
298        return r1, r2, np.array([], dtype=bool), np.array([], dtype=bool)
299
300
301def salt_bridge_plus_minus(mol1, mol2, cutoff=4):
302    """Returns pairs of plus-mins atoms, which meet salt bridge criteria
303
304    Parameters
305    ----------
306    mol1, mol2 : oddt.toolkit.Molecule object
307        Molecules to compute plus and minus pairs
308
309    cutoff : float, (default=4)
310        Distance cutoff for A-H pairs
311
312    Returns
313    -------
314    plus, minus : atom_dict-type numpy array
315        Aligned arrays of atoms forming salt bridge, firstly plus, secondly minus
316
317    """
318    m1_plus, m2_minus = close_contacts(mol1.atom_dict[mol1.atom_dict['isplus']],
319                                       mol2.atom_dict[mol2.atom_dict['isminus']],
320                                       cutoff)
321    return m1_plus, m2_minus
322
323
324def salt_bridges(mol1, mol2, *args, **kwargs):
325    """Calculates salt bridges between molecules
326
327    Parameters
328    ----------
329    mol1, mol2 : oddt.toolkit.Molecule object
330        Molecules to compute plus and minus pairs
331
332    cutoff : float, (default=4)
333        Distance cutoff for plus-minus pairs
334
335    Returns
336    -------
337    mol1_atoms, mol2_atoms : atom_dict-type numpy array
338        Aligned arrays of atoms forming salt bridges
339    """
340    m1_plus, m2_minus = salt_bridge_plus_minus(mol1, mol2, *args, **kwargs)
341    m2_plus, m1_minus = salt_bridge_plus_minus(mol2, mol1, *args, **kwargs)
342    return np.concatenate((m1_plus, m1_minus)), np.concatenate((m2_minus, m2_plus))
343
344
345def hydrophobic_contacts(mol1, mol2, cutoff=4):
346    """Calculates hydrophobic contacts between molecules
347
348    Parameters
349    ----------
350    mol1, mol2 : oddt.toolkit.Molecule object
351        Molecules to compute hydrophobe pairs
352
353    cutoff : float, (default=4)
354        Distance cutoff for hydrophobe pairs
355
356    Returns
357    -------
358    mol1_atoms, mol2_atoms : atom_dict-type numpy array
359        Aligned arrays of atoms forming hydrophobic contacts
360
361    """
362    h1, h2 = close_contacts(mol1.atom_dict[mol1.atom_dict['ishydrophobe']],
363                            mol2.atom_dict[mol2.atom_dict['ishydrophobe']],
364                            cutoff)
365    return h1, h2
366
367
368def pi_cation(mol1, mol2, cutoff=5, tolerance=30):
369    """Returns pairs of ring-cation atoms, which meet pi-cation criteria
370
371    Parameters
372    ----------
373    mol1, mol2 : oddt.toolkit.Molecule object
374        Molecules to compute ring-cation pairs
375
376    cutoff : float, (default=5)
377        Distance cutoff for Pi-cation pairs
378
379    tolerance : int, (default=30)
380        Range (+/- tolerance) from perfect direction (perpendicular)
381        in which pi-cation are considered as strict.
382
383    Returns
384    -------
385    r1 : ring_dict-type numpy array
386        Aligned rings forming pi-stacking
387
388    plus2 : atom_dict-type numpy array
389        Aligned cations forming pi-cation
390
391    strict_parallel : numpy array, dtype=bool
392        Boolean array align with ring-cation pairs, informing whether
393        they form 'strict' pi-cation. If false, only distance cutoff is met,
394        therefore the interaction is 'crude'.
395
396    """
397    r1, plus2 = close_contacts(mol1.ring_dict,
398                               mol2.atom_dict[mol2.atom_dict['isplus']],
399                               cutoff,
400                               x_column='centroid')
401    if len(r1) > 0 and len(plus2) > 0:
402        angle1 = angle_2v(r1['vector'], plus2['coords'] - r1['centroid'])
403        strict = (angle1 > 180 - tolerance) | (angle1 < tolerance)
404        return r1, plus2, strict
405    else:
406        return r1, plus2, np.array([], dtype=bool)
407
408
409def acceptor_metal(mol1, mol2, base_angle=120, tolerance=30, cutoff=4):
410    """Returns pairs of acceptor-metal atoms, which meet metal coordination criteria
411    Note: This function is directional (mol1 holds acceptors, mol2 holds metals)
412
413    Parameters
414    ----------
415    mol1, mol2 : oddt.toolkit.Molecule object
416        Molecules to compute acceptor and metal pairs
417
418    cutoff : float, (default=4)
419        Distance cutoff for A-M pairs
420
421    base_angle : int, (default=120)
422        Base angle determining allowed direction of metal coordination,
423        which is devided by the number of neighbors of acceptor atom
424        to establish final directional angle
425
426    tolerance : int, (default=30)
427        Range (+/- tolerance) from perfect direction (base_angle/n_neighbors)
428        in metal coordination are considered as strict.
429
430    Returns
431    -------
432    a, d : atom_dict-type numpy array
433        Aligned arrays of atoms forming metal coordination,
434        firstly acceptors, secondly metals.
435
436    strict : numpy array, dtype=bool
437        Boolean array align with atom pairs, informing whether atoms
438        form 'strict' metal coordination (pass all angular cutoffs).
439        If false, only distance cutoff is met, therefore the interaction
440        is 'crude'.
441    """
442    a, m = close_contacts(mol1.atom_dict[mol1.atom_dict['isacceptor']],
443                          mol2.atom_dict[mol2.atom_dict['ismetal']],
444                          cutoff)
445    # skip empty values
446    if len(a) > 0 and len(m) > 0:
447        angle1 = angle(m['coords'][:, np.newaxis, :],
448                       a['coords'][:, np.newaxis, :],
449                       a['neighbors'])
450        a_neighbors_num = np.sum(~np.isnan(a['neighbors'][:, :, 0]), axis=-1)[:, np.newaxis]
451        strict = ((np.nan_to_num(angle1) > (base_angle / a_neighbors_num - tolerance)) |
452                  np.isnan(angle1)).all(axis=-1)
453        return a, m, strict
454    else:
455        return a, m, np.array([], dtype=bool)
456
457
458def pi_metal(mol1, mol2, cutoff=5, tolerance=30):
459    """Returns pairs of ring-metal atoms, which meet pi-metal criteria
460
461    Parameters
462    ----------
463    mol1, mol2 : oddt.toolkit.Molecule object
464        Molecules to compute ring-metal pairs
465
466    cutoff : float, (default=5)
467        Distance cutoff for Pi-metal pairs
468
469    tolerance : int, (default=30)
470        Range (+/- tolerance) from perfect direction (perpendicular)
471        in which pi-metal are considered as strict.
472
473    Returns
474    -------
475    r1 : ring_dict-type numpy array
476        Aligned rings forming pi-metal
477
478    m : atom_dict-type numpy array
479        Aligned metals forming pi-metal
480
481    strict_parallel : numpy array, dtype=bool
482        Boolean array align with ring-metal pairs, informing whether
483        they form 'strict' pi-metal. If false, only distance cutoff is met,
484        therefore the interaction is 'crude'.
485
486    """
487    r1, m = close_contacts(mol1.ring_dict,
488                           mol2.atom_dict[mol2.atom_dict['ismetal']],
489                           cutoff,
490                           x_column='centroid')
491    if len(r1) > 0 and len(m) > 0:
492        angle1 = angle_2v(r1['vector'], m['coords'] - r1['centroid'])
493        strict = (angle1 > 180 - tolerance) | (angle1 < tolerance)
494        return r1, m, strict
495    else:
496        return r1, m, np.array([], dtype=bool)
497