RMSD Calculation

[Maria Thereza]

Hi everyone,

I am trying to calculate RMSD values for my trajectory. My simulation box contains a sugar and many molecules of water. I tried to account for PBC effects through this:

# define universe with MDAnalysis
u = mda.Universe('md.tpr', 'md.trr')

# define sugar and water molecules
sugar = u.select_atoms('resname UNL')
water = u.select_atoms('resname SOL')

transforms = [trans.unwrap(sugar),
                      trans.center_in_box(sugar, wrap=True),
                      trans.wrap(water)]

u.trajectory.add_transformations(*transforms)

And then I first calculate RMSD values taking the trajectory itself as reference state:

RMSD = mda.analysis.rms.RMSD(u, u , select='resname UNL', ref_frame=0)
RMSD.run()

# retrieve the RMSD values
rmsd_values = RMSD.rmsd.T

I have values that oscillate from 0 to 1.75. My question is: aren't the RMSD values for this case (analyzing the trajectory with respect to the trajectory itself) supposed to be 0?

Then I do pretty much the same thing, but taking as reference a .pdb file for the molecule:

# load the reference structure
ref_103 = mda.Universe("103.pdb")

# initialize an empty list to store RMSD values
rmsd_values_103 = []

# define a selection for the sugar molecule
selection = "resname UNL"

# create a RMSD analysis object
rmsd_analysis_103 = mda.analysis.rms.RMSD(u, select=selection, ref=ref_103, center=True, superposition=True)
rmsd_analysis_103.run()

# retrieve the RMSD values
rmsd_values_103 = rmsd_analysis_103.rmsd.T

And I obtain exactly the same output of the analysis of the trajectory itself. Also, whenever I try to calculate the RMSD values for other .pdb files (that correspond to other conformations of the molecule), I still get the same output.

What am I doing wrong?

Thank you very much!

Have a great day :)

[Fiona Naughton]

Hi!

> I have values that oscillate from 0 to 1.75. My question is: aren't the RMSD values for this case (analyzing the trajectory with respect to the trajectory itself) supposed to be 0?

The RMSD calculation always uses a single frame as a reference, even if the reference is itself a trajectory - so when you use the trajectory as it's own reference it's comparing each frame to a selected frame within the trajectory (in this case frame 0, since you have ref_frane=0), rather than comparing each frame to itself. 

> rmsd_analysis_103 = mda.analysis.rms.RMSD(u, select=selection, ref=ref_103, center=True, superposition=True)

> And I obtain exactly the same output of the analysis of the trajectory itself. Also, whenever I try to calculate the RMSD values for other .pdb files (that correspond to other conformations of the molecule), I still get the same output.

This is a little confusing, but the keyword for specifying a reference structure is "reference", not "ref" (see the docs here https://docs.mdanalysis.org/stable/documentation_pages/analysis/rms.html#MDAnalysis.analysis.rms.RMSD ). It looks like it accepts "ref" (rather than throwing an error) because the function is set to accept other keyword arguments, but this won't end up being used; and since then there's no input for the "reference" argument it defaults to using the trajectory itself as the reference - so you get the same result as your first example. If you change "ref" to "reference" in your code, it should work!

Hope that cleared things up, let us know if tou have any further questions! 

Best,

- Fiona

[Maria Thereza]

Hello Fiona,

Thank you very much! It worked by putting reference instead of ref, silly me. 

Have a great day :)

Best regards,

Maria