Compute bond lengths along paths

This example demonstrates how to compute minimum interatomic distance for each frame in each path and plot the result. In this example we consider all the bonds, but you can choose some specific bonds 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_bond_length_stats(bond_indexer):
    '''
    Compute min, max, and avg length of bonds
    Return results in angstrom units.
    '''
    if len(bond_indexer) == 0:
        return []

    min_bl = max_bl = bond_indexer[0].getLength()   # initialize min and max
    avg_bl = SBQuantity.length(0.0)
    for b in bond_indexer:               # loop over bonds in the bond indexer
        bl = b.getLength().angstrom      # get the bond length in Angstroms
        avg_bl = avg_bl + bl
        if bl < min_bl : min_bl = bl
        if bl > max_bl : max_bl = bl

    avg_bl = avg_bl / len(bond_indexer)

    return [max_bl.angstrom.value, avg_bl.angstrom.value, min_bl.angstrom.value]

def compute_bond_length_stats_per_path(path, bond_indexer):
    '''
    Compute min, max, and avg length of bonds along the path
    '''

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

    bond_length_along_path = []

    # loop over path frames
    for i in range(path.numberOfSteps):
        # set the current step of the path
        path.currentStep = i
        # compute bond length stats for the current step
        bond_length_along_path.append(compute_bond_length_stats(bond_indexer))
        # If you want to see in the viewport how the path changes, force SAMSON to process events.
        # Else comment it for bigger systems, since it may slow down the computations
        SAMSON.processEvents()

    # restore the current step
    path.currentStep = current_step

    return bond_length_along_path

# get an indexer of all paths in the active document
path_indexer = SAMSON.getNodes('node.type path')
# get an indexer of all bonds in the active document
bond_indexer = SAMSON.getNodes('node.type bond')

# Compute bond length stats for each path in the document
bond_length_along_paths = []
# loop over paths in the path indexer
for path in path_indexer:
    res = compute_bond_length_stats_per_path(path, bond_indexer)
    bond_length_along_paths.append( res )

# Plot results using matplotlib
for i, path in enumerate(path_indexer):
    plt.figure(i)
    plt.plot(bond_length_along_paths[i])
    plt.title(path.name)
    plt.legend(['max','mean','min'])
    plt.ylabel('bond length, Å')
    plt.xlabel('frame')
    plt.show()