Compute interatomic distances along path

This example demonstrates how to compute the minimum interatomic distance for each frame in each path and plot the result. In this example we consider all the atoms, but you can choose some specific atoms using NSL.

First, open a trajectory in SAMSON or fetch one from RCSB PDB, e.g. SAMSON.fetchPDB('1d3z').

import numpy as np
import matplotlib.pyplot as plt

def compute_min_interatomic_distance(atom_indexer, neighbor_search, rcutoff = SBQuantity.angstrom(2)):
    '''
    Compute the minimal interatomic distances between atoms in atom_indexer
    Return result in angstrom units.
    '''
    neighbor_search.updateNeighborLists()

    # initialize the minimum distance to rcutoff
    min_distance = rcutoff
    # go through all the atoms
    for i, atom in enumerate(atom_indexer):
        neighbors = neighbor_search.getNeighborVector(i)
        for neighbor_atom in neighbors:
            distance = (atom.getPosition() - neighbor_atom.getPosition()).norm()
            if distance < min_distance:
                min_distance = distance

    return min_distance.angstrom.value

def compute_min_interatomic_distance_per_path(path, atom_indexer, rcutoff = SBQuantity.angstrom(2.0)):
    '''
    Compute interatomic distances between atoms in atom_indexer along the path
    '''
    # Create a particle system from the atom indexer
    particle_system = SBParticleSystem(atom_indexer)
    # Create the neighbor search for this particle system
    neighbor_search = SBNeighborSearchParticleSystemGrid(particle_system, rcutoff)
    neighbor_search.initializeNeighborLists()

    # store the current step
    current_step = path.currentStep
    path.currentStep = 0

    interatomic_length = []
    # loop over path steps
    for i in range(path.numberOfSteps):
        # set the current step of the path
        path.currentStep = i
        interatomic_length.append( compute_min_interatomic_distance(atom_indexer, neighbor_search, rcutoff) )
        # If you want to see in the viewport how the path changes, force SAMSON to process events
        # Comment it for bigger systems, since it may slow down the computations
        SAMSON.processEvents()

    # restore the current step
    path.currentStep = current_step

    return interatomic_length

# get an indexer of all paths in the active document
path_indexer = SAMSON.getNodes('node.type path')
# get an indexer of all atoms in the active document
atom_indexer = SAMSON.getNodes('node.type atom')

# set the radius of cutoff
rcutoff = SBQuantity.angstrom(2.0)

# Compute min interatomic distance for each path in the document
interatomic_length_per_path = []
# loop over paths in the path indexer
for path in path_indexer:
    value = compute_min_interatomic_distance_per_path(path, atom_indexer, rcutoff)
    interatomic_length_per_path.append( value )

# Plot results using matplotlib
legend = []
for i, path in enumerate(path_indexer):
    plt.plot(interatomic_length_per_path[i])
    legend.append(path.name)

plt.legend(legend)
plt.ylabel('min interatomic distance, Å')
plt.xlabel('frame')
plt.tight_layout()
plt.show()