4.2.4. Calculating root mean square quantities — MDAnalysis.analysis.rms
¶
 Author
Oliver Beckstein, David L. Dotson, John Detlefs
 Year
2016
 Copyright
GNU Public License v2
New in version 0.7.7.
Changed in version 0.11.0: Added RMSF
analysis.
Changed in version 0.16.0: Refactored RMSD to fit AnalysisBase API
The module contains code to analyze root mean square quantities such
as the coordinat root mean square distance (RMSD
) or the
perresidue root mean square fluctuations (RMSF
).
This module uses the fast QCP algorithm [Theobald2005] to calculate
the root mean square distance (RMSD) between two coordinate sets (as
implemented in
MDAnalysis.lib.qcprot.CalcRMSDRotationalMatrix()
).
When using this module in published work please cite [Theobald2005].
See also
MDAnalysis.analysis.align
aligning structures based on RMSD
MDAnalysis.lib.qcprot
implements the fast RMSD algorithm.
4.2.4.1. Example applications¶
4.2.4.1.1. Calculating RMSD for multiple domains¶
In this example we will globally fit a protein to a reference structure and investigate the relative movements of domains by computing the RMSD of the domains to the reference. The example is a DIMS trajectory of adenylate kinase, which samples a large closedtoopen transition. The protein consists of the CORE, LID, and NMP domain.
superimpose on the closed structure (frame 0 of the trajectory), using backbone atoms
calculate the backbone RMSD and RMSD for CORE, LID, NMP (backbone atoms)
The trajectory is included with the test data files. The data in
RMSD.rmsd
is plotted with matplotlib.pyplot.plot()
:
import MDAnalysis
from MDAnalysis.tests.datafiles import PSF,DCD,CRD
u = MDAnalysis.Universe(PSF,DCD)
ref = MDAnalysis.Universe(PSF,DCD) # reference closed AdK (1AKE) (with the default ref_frame=0)
#ref = MDAnalysis.Universe(PSF,CRD) # reference open AdK (4AKE)
import MDAnalysis.analysis.rms
R = MDAnalysis.analysis.rms.RMSD(u, ref,
select="backbone", # superimpose on whole backbone of the whole protein
groupselections=["backbone and (resid 129 or resid 60121 or resid 160214)", # CORE
"backbone and resid 122159", # LID
"backbone and resid 3059"], # NMP
filename="rmsd_all_CORE_LID_NMP.dat")
R.run()
import matplotlib.pyplot as plt
rmsd = R.rmsd.T # transpose makes it easier for plotting
time = rmsd[1]
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.plot(time, rmsd[2], 'k', label="all")
ax.plot(time, rmsd[3], 'k', label="CORE")
ax.plot(time, rmsd[4], 'r', label="LID")
ax.plot(time, rmsd[5], 'b', label="NMP")
ax.legend(loc="best")
ax.set_xlabel("time (ps)")
ax.set_ylabel(r"RMSD ($\AA$)")
fig.savefig("rmsd_all_CORE_LID_NMP_ref1AKE.pdf")
4.2.4.2. Functions¶

MDAnalysis.analysis.rms.
rmsd
(a, b, weights=None, center=False, superposition=False)[source]¶ Returns RMSD between two coordinate sets a and b.
a and b are arrays of the coordinates of N atoms of shape \(N times 3\) as generated by, e.g.,
MDAnalysis.core.groups.AtomGroup.positions()
.Note
If you use trajectory data from simulations performed under periodic boundary conditions then you must make your molecules whole before performing RMSD calculations so that the centers of mass of the mobile and reference structure are properly superimposed.
 Parameters
a (array_like) – coordinates to align to b
b (array_like) – coordinates to align to (same shape as a)
weights (array_like (optional)) – 1D array with weights, use to compute weighted average
center (bool (optional)) – subtract center of geometry before calculation. With weights given compute weighted average as center.
superposition (bool (optional)) – perform a rotational and translational superposition with the fast QCP algorithm [Theobald2005] before calculating the RMSD; implies
center=True
.
 Returns
rmsd – RMSD between a and b
 Return type
Notes
The RMSD \(\rho(t)\) as a function of time is calculated as
\[\rho(t) = \sqrt{\frac{1}{N} \sum_{i=1}^N w_i \left(\mathbf{x}_i(t)  \mathbf{x}_i^{\text{ref}}\right)^2}\]It is the Euclidean distance in configuration space of the current configuration (possibly after optimal translation and rotation) from a reference configuration divided by \(1/\sqrt{N}\) where \(N\) is the number of coordinates.
The weights \(w_i\) are calculated from the input weights weights \(w'_i\) as relative to the mean:
\[w_i = \frac{w'_i}{\langle w' \rangle}\]Example
>>> u = Universe(PSF,DCD) >>> bb = u.select_atoms('backbone') >>> A = bb.positions.copy() # coordinates of first frame >>> u.trajectory[1] # forward to last frame >>> B = bb.positions.copy() # coordinates of last frame >>> rmsd(A, B, center=True) 3.9482355416565049
4.2.4.3. Analysis classes¶

class
MDAnalysis.analysis.rms.
RMSD
(atomgroup, reference=None, select='all', groupselections=None, filename='rmsd.dat', weights=None, tol_mass=0.1, ref_frame=0, **kwargs)[source]¶ Class to perform RMSD analysis on a trajectory.
The RMSD will be computed for two groups of atoms and all frames in the trajectory belonging to atomgroup. The groups of atoms are obtained by applying the selection selection select to the changing atomgroup and the fixed reference.
Note
If you use trajectory data from simulations performed under periodic boundary conditions then you must make your molecules whole before performing RMSD calculations so that the centers of mass of the selected and reference structure are properly superimposed.
Run the analysis with
RMSD.run()
, which stores the results in the arrayRMSD.rmsd
. Parameters
atomgroup (AtomGroup or Universe) – Group of atoms for which the RMSD is calculated. If a trajectory is associated with the atoms then the computation iterates over the trajectory.
reference (AtomGroup or Universe (optional)) – Group of reference atoms; if
None
then the current frame of atomgroup is used.select (str or dict or tuple (optional)) –
The selection to operate on; can be one of:
any valid selection string for
select_atoms()
that produces identical selections in atomgroup and reference; ora dictionary
{'mobile': sel1, 'reference': sel2}
where sel1 and sel2 are valid selection strings that are applied to atomgroup and reference respectively (theMDAnalysis.analysis.align.fasta2select()
function returns such a dictionary based on a ClustalW or STAMP sequence alignment); ora tuple
(sel1, sel2)
When using 2. or 3. with sel1 and sel2 then these selection strings are applied to atomgroup and reference respectively and should generate groups of equivalent atoms. sel1 and sel2 can each also be a list of selection strings to generate a
AtomGroup
with defined atom order as described under Ordered selections).groupselections (list (optional)) –
A list of selections as described for select, with the difference that these selections are always applied to the full universes, i.e.,
atomgroup.universe.select_atoms(sel1)
andreference.universe.select_atoms(sel2)
. Each selection describes additional RMSDs to be computed after the structures have been superimposed according to select. No additional fitting is performed.The output contains one additional column for each selection.Note
Experimental feature. Only limited error checking implemented.
filename (str (optional)) –
write RMSD into file with
RMSD.save()
weights ({“mass”,
None
} or array_like (optional)) – choose weights. With"mass"
uses masses as weights; withNone
weigh each atom equally. If a float array of the same length as atomgroup is provided, use each element of the array_like as a weight for the corresponding atom in atomgroup.tol_mass (float (optional)) – Reject match if the atomic masses for matched atoms differ by more than tol_mass.
ref_frame (int (optional)) – frame index to select frame from reference
verbose (bool (optional)) – Show detailed progress of the calculation if set to
True
; the default isFalse
.
 Raises
SelectionError – If the selections from atomgroup and reference do not match.
TypeError – If weights is not of the appropriate type; see also
MDAnalysis.lib.util.get_weights()
ValueError – If weights are not compatible with atomgroup (not the same length) or if it is not a 1D array (see
MDAnalysis.lib.util.get_weights()
). AValueError
is also raised if weights are not compatible with groupselections: only equal weights (weights=None
) or massweighted (weights="mass"
) are supported for additional groupselections.
Notes
The root mean square deviation \(\rho(t)\) of a group of \(N\) atoms relative to a reference structure as a function of time is calculated as
\[\rho(t) = \sqrt{\frac{1}{N} \sum_{i=1}^N w_i \left(\mathbf{x}_i(t)  \mathbf{x}_i^{\text{ref}}\right)^2}\]The weights \(w_i\) are calculated from the input weights weights \(w'_i\) as relative to the mean of the input weights:
\[w_i = \frac{w'_i}{\langle w' \rangle}\]The selected coordinates from atomgroup are optimally superimposed (translation and rotation) on the reference coordinates at each time step as to minimize the RMSD. Douglas Theobald’s fast QCP algorithm [Theobald2005] is used for the rotational superposition and to calculate the RMSD (see
MDAnalysis.lib.qcprot
for implementation details).The class runs various checks on the input to ensure that the two atom groups can be compared. This includes a comparison of atom masses (i.e., only the positions of atoms of the same mass will be considered to be correct for comparison). If masses should not be checked, just set tol_mass to a large value such as 1000.
See also
New in version 0.7.7.
Changed in version 0.8: groupselections added
Changed in version 0.16.0: Flexible weighting scheme with new weights keyword.
Deprecated since version 0.16.0: Instead of
mass_weighted=True
(removal in 0.17.0) use newweights='mass'
; refactored to fit with AnalysisBase APIChanged in version 0.17.0: removed deprecated mass_weighted keyword; groupselections are not rotationally superimposed any more.
Deprecated since version 0.19.0: filename will be removed in 1.0

rmsd
¶ Contains the time series of the RMSD as an N×3
numpy.ndarray
array with content[[frame, time (ps), RMSD (A)], [...], ...]
.

run
(start=None, stop=None, step=None, verbose=None)¶ Perform the calculation

save
(**kwds)¶ save is deprecated!
Save RMSD from
RMSD.rmsd
to text file filename. Parameters
filename (str (optional)) – if no filename is given the default provided to the constructor is used.
Deprecated since version 0.19.0: You can instead use
np.savetxt(filename, RMSD.rmsd)
. save will be removed in release 1.0.0.

class
MDAnalysis.analysis.rms.
RMSF
(atomgroup, **kwargs)[source]¶ Calculate RMSF of given atoms across a trajectory.
Note
No RMSDsuperposition is performed; it is assumed that the user is providing a trajectory where the protein of interest has been structurally aligned to a reference structure (see the Examples section below). The protein also has be whole because periodic boundaries are not taken into account.
Run the analysis with
RMSF.run()
, which stores the results in the arrayRMSF.rmsf
. Parameters
atomgroup (AtomGroup) – Atoms for which RMSF is calculated
start (int (optional)) – starting frame, default None becomes 0.
stop (int (optional)) – Frame index to stop analysis. Default: None becomes n_frames. Iteration stops before this frame number, which means that the trajectory would be read until the end.
step (int (optional)) – step between frames, default None becomes 1.
verbose (bool (optional)) – Show detailed progress of the calculation if set to
True
; the default isFalse
.
 Raises
ValueError – raised if negative values are calculated, which indicates that a numerical overflow or underflow occured
Notes
The root mean square fluctuation of an atom \(i\) is computed as the time average
\[\rho_i = \sqrt{\left\langle (\mathbf{x}_i  \langle\mathbf{x}_i\rangle)^2 \right\rangle}\]No mass weighting is performed.
This method implements an algorithm for computing sums of squares while avoiding overflows and underflows [Welford1962].
Examples
In this example we calculate the residue RMSF fluctuations by analyzing the \(\text{C}_\alpha\) atoms. First we need to fit the trajectory to the average structure as a reference. That requires calculating the average structure first. Because we need to analyze and manipulate the same trajectory multiple times, we are going to load it into memory using the
MemoryReader
. (If your trajectory does not fit into memory, you will need to write out intermediate trajectories to disk or generate an inmemory universe that only contains, say, the protein):import MDAnalysis as mda from MDAnalysis.analysis import align from MDAnalysis.tests.datafiles import TPR, XTC u = mda.Universe(TPR, XTC, in_memory=True) protein = u.select_atoms("protein") # 1) need a step to center and make whole: this trajectory # contains the protein being split across periodic boundaries # # TODO # 2) fit to the initial frame to get a better average structure # (the trajectory is changed in memory) prealigner = align.AlignTraj(u, select="protein and name CA", in_memory=True).run() # 3) reference = average structure reference_coordinates = u.trajectory.timeseries(asel=protein).mean(axis=1) # make a reference structure (need to reshape into a 1frame "trajectory") reference = mda.Merge(protein).load_new( reference_coordinates[:, None, :], order="afc")
We created a new universe
reference
that contains a single frame with the averaged coordinates of the protein. Now we need to fit the whole trajectory to the reference by minimizing the RMSD. We useMDAnalysis.analysis.align.AlignTraj
:aligner = align.AlignTraj(u, reference, select="protein and name CA", in_memory=True).run()
The trajectory is now fitted to the reference (the RMSD is stored as aligner.rmsd for further inspection). Now we can calculate the RMSF:
from MDAnalysis.analysis.rms import RMSF calphas = protein.select_atoms("name CA") rmsfer = RMSF(calphas, verbose=True).run()
and plot:
import matplotlib.pyplot as plt plt.plot(calphas.resnums, rmsfer.rmsf)
References
 Welford1962
B. P. Welford (1962). “Note on a Method for Calculating Corrected Sums of Squares and Products.” Technometrics 4(3):419420.
New in version 0.11.0.
Changed in version 0.16.0: refactored to fit with AnalysisBase API
Deprecated since version 0.16.0: the keyword argument quiet is deprecated in favor of verbose.
Changed in version 0.17.0: removed unused keyword weights

rmsf
¶ Results are stored in this Nlength
numpy.ndarray
array, giving RMSFs for each of the given atoms.

run
(start=None, stop=None, step=None, verbose=None)¶ Perform the calculation