Get neighbors#

Finding neighboring atoms or residue groups within a distance threshold or by nearest count.

The function molsysmt.structure.get_neighbors() finds surrounding spatial neighbors for a target set of atoms or groups, either by specifying the \(k\) nearest neighbors (n_neighbors) or a distance cutoff (threshold).

Added in version 1.0.0.

API documentation

Follow this link for a detailed description of the input arguments, raised errors, and returned objects of this function: molsysmt.structure.get_neighbors().

Basic usage#

Let’s show how to query spatial neighbors for a target set of atoms (illustrated using \(C\alpha\) atoms in TcTIM, PDB ID: 1TCD):

import molsysmt as msm
import matplotlib.pyplot as plt
import numpy as np

molsys = msm.convert(msm.systems['TcTIM']['1tcd.h5msm'])

Note

Demo Systems Catalog

This tutorial uses demonstration datasets provided by MolSysMT. To explore the full catalog of bundled systems, forms, and file paths, visit the Demo Systems guide.

Querying fixed nearest neighbors#

We query the 3 closest spatial neighbors for each \(C\alpha\) atom using n_neighbors=3:

ca_atoms = msm.select(molsys, selection='atom_name=="CA"')
neighbors, distances = msm.structure.get_neighbors(molsys, selection=ca_atoms, n_neighbors=3)

print('Neighbors array shape:', neighbors.shape)
print('Distances array shape:', distances.shape)
print('Closest 3 neighbors of first CA atom:', neighbors[0, 0])
print('Distances to closest 3 neighbors:', distances[0, 0])

Neighbors array shape: (1, 497, 3)
Distances array shape: (1, 497, 3)
Closest 3 neighbors of first CA atom: [  1   2 226]
Distances to closest 3 neighbors: [0.37897700341900714 0.6656387853654624 0.7216471108770198] nanometer

Cutoff-based neighbor search#

Instead of a fixed number of neighbors, we can query all atoms located within a spherical cutoff distance (e.g. threshold='0.6 nm'):

neighbors_cut, distances_cut = msm.structure.get_neighbors(molsys, selection=ca_atoms, threshold='0.6 nm')
n_neighbors_per_ca = [len(n_list) for n_list in neighbors_cut[0]]
print(f'Average neighbors within 0.6 nm per CA atom: {np.mean(n_neighbors_per_ca):.1f}')

Average neighbors within 0.6 nm per CA atom: 5.6

Visualizing neighbor count distribution#

We plot the distribution of neighbor counts across elements to analyze local spatial packing density:

plt.figure(figsize=(9, 4.5))
plt.hist(n_neighbors_per_ca, bins=range(min(n_neighbors_per_ca), max(n_neighbors_per_ca) + 2), color='royalblue', edgecolor='black', alpha=0.75, rwidth=0.85)
plt.xlabel('Number of Neighbors within Cutoff')
plt.ylabel('Element Count')
plt.title('Distribution of Local Packing Density')
plt.grid(axis='y', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

../../../../_images/ff7a20cf33257fba8c48556d5ea05363304679de440e14e628915d1018ac692c.png

Neighbor search between different atom selections#

We can find neighbors in an external target selection (selection_2) around a reference query selection (illustrated here using interface contacts between Chain 0 and Chain 1):

ca_chain0 = msm.select(molsys, selection='chain_index==0 and atom_name=="CA"')
ca_chain1 = msm.select(molsys, selection='chain_index==1 and atom_name=="CA"')
interface_neigh, interface_dist = msm.structure.get_neighbors(molsys, selection=ca_chain0, selection_2=ca_chain1, threshold='0.8 nm')

contacting_residues = sum(1 for n in interface_neigh[0] if len(n) > 0)
print(f'Reference elements with at least one neighbor in target set (<0.8 nm): {contacting_residues} / {len(ca_chain0)}')

Reference elements with at least one neighbor in target set (<0.8 nm): 34 / 248