Calculating the End-to-End distance of two atoms using MDAnalysis

If you’ve come across this article, you must be a python-dev (like myself) that has come across MDAnalysis but reading over the docs you’ve discovered you know nothing about molecular dynamics. So let’s start at the basics: Take the distance between two atoms in space over a trajectory.

I’ve got two atoms in space first I’m going to load up my topology and trajectory using mda and initiate the universe (sounds so cool):

def initiate_universe(topology, trajectory):
    universe = mda.Universe(topology, trajectory)
    return universe
if __name__ == '__main__':
    # File Setup
    current_directory = os.getcwd()
    topology = os.path.join(current_directory, 'topology.psf')
    trajectory = os.path.join(current_directory, 'trajectory.xtc')
    
    initiate_universe(topology, trajectory)

Simple, lets take a look at two specific atoms as they traverse through the trajectory in my case it’s two oxygens and we’ll quickly view the atoms using the nglviewer:

oxygen_1 = universe.select_atoms('segid A and resid 1 and name O1')
oxygen_4 = universe.select_atoms('segid A and resid 10 and name O4')
merge = mda.Merge(oxygen_1, oxygen_4)
depiction = nv.show_mdanalysis(merge)
depiction

In my specific case, I have one chain that is 10 monomer units long and it uses the CHARMM selecting language of atoms and I use the merge function in mda to combine the two atoms universe for the viewer.

Note: I don’t actually know if that is the correct way to do this — but it works for me :)

Looks beautiful. So how do we calculate a distance from two atoms in space with three dimensions (x, y, z)

The End-to-End distance is calculated using the difference between the two vectors positions in the x, y, z , square it, and take a square root of the sum. Who came up with this groundbreaking math? Pythagoras!

What does the code look like?

def find_magnitude(x, y, z):
    
    # imports
    # -------
    import math
    difference_x = x[1] - x[0]
    difference_y = y[1] - y[0]
    difference_z = z[1] - z[0]
    difference_x = difference_x ** 2
    difference_y = difference_y ** 2
    difference_z = difference_z ** 2
    magnitude = math.sqrt(difference_x + difference_y + difference_z)
    return magnitude

And we have to set up our “domains of interest” which are the two atoms.

def setup_domains_of_interest(universe):
        
        domains = {}
        domains['O1'] = universe.select_atoms('segid A and resid 1 and name O1')
        domains['O4'] = universe.select_atoms('segid A and  resid 10 and name O4')
        
    return domains

I am selecting segid (Chain A), resid (Residue 1 and 10), name (Oxygen 1 and Oxygen 4). We want to store our domains as a dictionary and I believe there is some MDAnalysis funkiness about why we setup domains this way.

Note: So what I have discovered is that when iterating through a trajectory having a universe.select_atoms makes my loops super low in magnitudes in order of like 3 seconds vs 24 minutes.

Alright now we are cooking, let’s declare some variables

calculated_magnitudes = []
timestamps = []
domains = setup_domains_of_interest(universe)

And now let’s set up our for loop to iterate through the trajectory, get the positions of O1 and O4 and fetch the magnitude.

for timestep in universe.trajectory:
        
    x0 = domains['O1'].positions[0][0]
    x1 = domains['O4'].positions[0][0]
    y0 = domains['O1'].positions[0][1]
    y1 = domains['O4'].positions[0][1]
    z0 = domains['O1'].positions[0][2]
    z1 = domains['O4'].positions[0][2]
    x = [x0, x1]
    y = [y0, y1]
    z = [z0, z1]
        
    calculated_magnitudes.append(find_magnitude(x, y, z))
    timestamps.append(universe.trajectory.time)
        

If you want to be extra fancy you can wrap a progress bar so you’re not wondering if the code is running like me…

import tqdm
progressbar_trajectory = tqdm(total=len(universe.trajectory))
    
for timestep in universe.trajectory:
        
    x0 = domains['O1'].positions[0][0]
    x1 = domains['O4'].positions[0][0]
    y0 = domains['O1'].positions[0][1]
    y1 = domains['O4'].positions[0][1]
    z0 = domains['O1'].positions[0][2]
    z1 = domains['O4'].positions[0][2]
    x = [x0, x1]
    y = [y0, y1]
    z = [z0, z1]
        
    calculated_magnitudes[segid].append(find_magnitude(x, y, z))
    timestamps.append(universe.trajectory.time)
        
    progressbar_trajectory.update(1)
        
progressbar_trajectory.close()

Run it!

Looks good! Now you can have some fun plotting the data in any package you like.

Again, more big ups to Paween! Postdoc in Jana Shen’s Lab who I ask relentless questions too and always answers them.

Python Env (3.6.9)

nglview==2.7.7
tqdm==4.48.2
MDAnalysis==1.0.0