Velocity of residue movement calculation question

[KyleStiers]

Hello all,

I am trying to develop a script to show a loop moves much faster in one trajectory than it does in another. I currently am thinking the best way to show this is on a per-residue level. So I have written the below, which gives me an output CSV matrix of distance moved per (Ca) residue per frame - analogous to RMSF but not dependent on an average structure (because I'm interested in speed of change).

However, all the values are nearly identical. It seems like from the log that everything is going correctly, i.e. current and previous coordinates are updating correctly. Am I missing something obvious?

It would also be interesting and informative if I could find a center-of-mass for a group of residues and monitor the distance traversed per frame to compare between two trajectories - but I'm not sure yet how to do that.

NOTES: you may wonder why I wrote my own distance calculation, I had trouble getting dist() to work for some reason and 562 is the number of Ca residues over 1002 frames in each trajectory.

Would appreciate the help or suggestions on how to improve the methodology!

Code to reproduce the behaviour

import MDAnalysis
from MDAnalysis.analysis.distances import dist
import numpy.linalg
import numpy
import csv
import pandas as pd
import math

outpath = ""
u = MDAnalysis.Universe("frame.pdb", "traj_fit.xtc")

counter = 0
ca = u.select_atoms('name CA')
ca_dist = numpy.empty((1001,562))

def save_traj_csv(numpy_array, name):
    df = pd.DataFrame(numpy_array)
    df.to_csv(outpath+name)

def diff(curr, prev):
        return(math.sqrt((curr[0]-prev[0])**2 + (curr[1]-prev[1])**2 + (curr[2]-prev[2])**2))

for ts in u.trajectory:
    if(counter != 0):
        for i in range(0,561):
            ca_dist[counter,i] = diff(ca.positions[i], ca.positions[i-1])
            print("Ca1: "+str(ca.positions[i])+" - CA2:"+str(ca.positions[i-1])+" = "+str(ca_dist[counter,i]))
        counter += 1
    else:
        counter += 1

save_traj_csv(ca_dist, "ca_dist_vec.csv")

....

[KyleStiers]

#1053 somewhat related to this

[KyleStiers]

Confirmed it wasn't my distance calculation (that would've been embarrassing). Changing the calculation to using:

ca_dist[counter,i] = numpy.linalg.norm(ca.positions[i]-ca.positions[i-1])

yielded the same numbers.

[richardjgowers]

Hey @KyleStiers

I'm a little confused at what you're trying to calculate, are you trying to show that the Ca positions are moving faster over time? The line diff(ca.positions[i], ca.positions[i-1]) will compare one Ca to the next Ca at the same timestep. I think if you want to do differences over time you'll want something more like:

for counter, ts in enumerate(u.trajectory):
    if counter != 0:
        ca_dist[counter, :] = multi_diff(ca.positions - old_positions)
    old_positions = ca.positions

Where you're comparing the positions in the current timestep to the last

[KyleStiers]

Hi @richardjgowers

Doh. You are absolutely right. I actually had been using something exactly like what you posted, albeit less concise, then at some point changed it to my clearly incorrect function.

That and then including the nested loop over each residue inside got it, thanks for seeing my rookie mistake. Now on to trying to figure out the center_of_mass() stuff 👍

for counter, ts in enumerate(u.trajectory):
    if counter != 0:
        for i in range(0,562):
            ca_dist[counter,i] = numpy.linalg.norm(ca.positions[i] - old_positions[i])
    old_positions = ca.positions