Inquiry About Per-Atom or Residue-Wise RMSD Calculation in MDAnalysis

[royyraa]

Hi,

I am unsure if a per-atom RMSD calculation tool is available in MDAnalysis. I did find that a 2D RMSD calculation for trajectories is available. However, I want to calculate RMSD values either residue-wise or atom-wise between two PDB files or two frames.

Do you have any suggestions or guidance on how to approach this? Any help would be greatly appreciated.

[orbeckst]

There's no pre-built tool for you. But it shouldn't be difficult for you to calculate it yourself: After structural superposition of a trajectory frame (e.g., with the tools in analysis.align) you have access to all positions in u.atoms.positions.

from MDAnalysis.analysis.align import alignto

u = mda.Universe(PSF, DCD)
p = u.select_atoms("protein")

# need a different Universe for the reference so that coordinates are independent
ref = u.copy()
ref_p = ref.select_atoms("protein")
ref.trajectory[0]   # use first frame of ref

u.trajectory[-1]    # use last frame of trajectory

# superimpose current frame of u
alignto(p, ref_p)

For each atom $i$, the "RMSD" at time $t$ relative to a reference value at time $t'$ is simply the distance $d_{i}(t,t') = \sqrt{\left(\mathbf{r}_i(t) - \mathbf{r}_i^\mathrm{ref}(t')\right)^2}$ and you can calculate all distances for all atoms with

atom_distances = mda.lib.distances.calc_bonds(p, ref_p)

If you want atoms by residues, you can, for example, use groupby("resids") as in p = u.select_atoms("protein"); atoms_by_res = p.atoms.groupby("resids") where atoms_by_res is a dict with an AtomGroup for each residue.

p_byres = p.groupby("resids")
ref_p_byres = ref_p.groupby("resids")

# for residue 10
from MDAnalysis.analysis.rms import rmsd
rmsd_residue_10 = rmsd(ref_p_byres[10].positions, p_byres[10].positions)

You'll have to build your own loops over frames and residues but this should give you an idea what you can do.