rmsd

mdtools.structure.rmsd(refpos, selpos, weights=None, center=False, inplace=False, xyz=False, box=None, out_tmp=None)[source]

Calculate the Root Mean Square Deviation (RMSD) between two sets of positions.

Todo

Implement a “superposition” functionality like in MDAnalysis.analysis.rms.rmsd().

Parameters:

refpos, selpos (array_like) – Reference and candidate set of positions. Position arrays must be of shape (3,), (n, 3) or (k, n, 3) where n is the number of particles and k is the number of frames. refpos and selpos must contain the same number of particles n, but can have different numbers of frames k. If refpos and selpos do not have the same shape, they must be broadcastable to a common shape (which becomes the shape of the output). The user is responsible to provide inputs that result in a physically meaningful broadcasting!
weights (None or array_like, optional) – Array of shape (n,) containing the weight of each particle contained in the position arrays. If weights is None, all particles are assumed to have a weight equal to one.
center (bool, optional) – If True, shift refpos and selpos by their (weighted) center, respectively, before calculating the RMSD. Note that no coordinate wrapping is done while calculating the (weighted) centers even if box is supplied.
inplace (bool, optional) – If True, subtract the weighted center from the reference and candidate positions in place (i.e. the input arrays will be changed). Note that refpos and selpos must have an appropriate dtype in this case.
xyz (bool, optional) – If True, return the x, y and z component of the RMSD separately instead of summing all components. Note however, that in this case the square root is not taken, because you cannot extract the square root of each summand individually.
box (None or array_like, optional) – The unit cell dimensions of the system, which can be orthogonal or triclinic and must provided in the same format as returned by MDAnalysis.coordinates.base.Timestep.dimensions: [lx, ly, lz, alpha, beta, gamma]. box can also be an array of boxes of shape (k, 6), where k must match the number of frames in refpos. If given, the minimum image convention is taken into account when calculating the distance between the reference and candidate positions.
out_tmp (None or numpy.ndarray, optional) – Preallocated array of the given dtype into which temporary result are stored. If provided, it must have the same shape as the result of selpos - refpos. Providing out_tmp can speed up the calculated if this function is called many times.

Raises:

ValueError : – If weights is not of shape (n,) or if it sums up to zero.

See also

MDAnalysis.analysis.rms.RMSD: Class to perform RMSD analysis on a trajectory
MDAnalysis.analysis.rms.rmsd(): Calculate the RMSD between two coordinate sets

Notes

The root mean square deviation is calculated by

\[RMSD = \sqrt{\frac{1}{N} \frac{1}{W} \sum_{i=1}^N w_i \left( \mathbf{r}_i - \mathbf{r}_i^{ref} \right)^2}\]

where \(N\) is the number of particles, \(w_i\) is the weight of the \(i\)-th particle, \(W = \sum_{i=1}^N w_i\) is the sum of particle weights, \(\mathbf{r}_i\) are the candidate positions and \(\mathbf{r}_i^{ref}\) are the reference positions.

If xyz is True, each component of the RMSD is returned individually without taking the square root. For instance, the \(x\)-component is

\[\langle \Delta x^2 \rangle = \frac{1}{N} \frac{1}{W} \sum_{i=1}^N w_i \left( x_i - x_i^{ref} \right)^2.\]

Examples

Shape of refpos and selpos is (3,):

>>> refpos = np.array([ 3,  4,  5])
>>> selpos = np.array([ 1,  7, -1])
>>> selpos - refpos
array([-2,  3, -6])
>>> mdt.strc.rmsd(refpos, selpos)
7.0
>>> mdt.strc.rmsd(refpos, selpos, xyz=True)
array([ 4.,  9., 36.])
>>> box = np.array([5, 3, 3, 90, 90, 90])
>>> mdt.strc.rmsd(refpos, selpos, xyz=True, box=box)
array([4., 0., 0.])
>>> mdt.strc.rmsd(refpos, selpos, box=box)
2.0
>>> mdt.strc.rmsd(refpos, selpos, weights=[3])
7.0
>>> mdt.strc.rmsd(refpos, selpos, weights=[3], center=True)
0.0
>>> mdt.strc.rmsd(refpos, selpos, center=True, inplace=True)
0.0
>>> refpos
array([0, 0, 0])
>>> selpos
array([0, 0, 0])

Shape of refpos and selpos is (n, 3):

>>> refpos = np.array([[ 3.,  4.,  5.],
...                    [ 1., -1.,  1.]])
>>> selpos = np.array([[ 1.,  7., -1.],
...                    [-1.,  0., -1.]])
>>> selpos - refpos
array([[-2.,  3., -6.],
       [-2.,  1., -2.]])
>>> rmsd = mdt.strc.rmsd(refpos, selpos)
>>> rmsd_expected = np.sqrt(0.5 * (49 + 9))
>>> np.isclose(rmsd, rmsd_expected, rtol=0)
True
>>> mdt.strc.rmsd(refpos, selpos, xyz=True)
array([ 4.,  5., 20.])
>>> box = np.array([5, 3, 3, 90, 90, 90])
>>> mdt.strc.rmsd(refpos, selpos, xyz=True, box=box)
array([4. , 0.5, 0.5])
>>> weights = np.array([0.575, 0.425])
>>> rmsd = mdt.strc.rmsd(refpos, selpos, weights=weights)
>>> np.isclose(rmsd, 4, rtol=0)
True
>>> mdt.strc.rmsd(refpos, selpos, weights=weights, center=True)
1.5632498200863483
>>> mdt.strc.rmsd(
...     refpos, selpos, xyz=True, center=True, inplace=True
... )
array([0., 1., 4.])
>>> refpos
array([[ 1. ,  2.5,  2. ],
       [-1. , -2.5, -2. ]])
>>> selpos
array([[ 1. ,  3.5,  0. ],
       [-1. , -3.5,  0. ]])

Shape of refpos and selpos is (k, n, 3):

>>> refpos = np.array([[[ 3.,  4.,  5.],
...                     [ 1., -1.,  1.]],
...
...                    [[ 5.,  0.,  1.],
...                     [ 1.,  2., -3.]]])
>>> selpos = np.array([[[ 1.,  7., -1.],
...                     [-1.,  0., -1.]],
...
...                    [[ 3.,  4.,  5.],
...                     [ 1., -1.,  1.]]])
>>> selpos - refpos
array([[[-2.,  3., -6.],
        [-2.,  1., -2.]],

       [[-2.,  4.,  4.],
        [ 0., -3.,  4.]]])
>>> rmsd = mdt.strc.rmsd(refpos, selpos)
>>> rmsd_expected = np.sqrt(0.5 * np.array([49 + 9, 36 + 25]))
>>> np.allclose(rmsd, rmsd_expected, rtol=0)
True
>>> mdt.strc.rmsd(refpos, selpos, xyz=True)
array([[ 4. ,  5. , 20. ],
       [ 2. , 12.5, 16. ]])
>>> mdt.strc.rmsd(refpos, selpos, xyz=True, box=box)
array([[4. , 0.5, 0.5],
       [2. , 0.5, 1. ]])
>>> box = np.array([[5, 3, 3, 90, 90, 90],
...                 [3, 5, 5, 90, 90, 90]])
>>> mdt.strc.rmsd(refpos, selpos, xyz=True, box=box)
array([[4. , 0.5, 0.5],
       [0.5, 2.5, 1. ]])
>>> weights = np.array([0.575, 0.425])
>>> rmsd = mdt.strc.rmsd(refpos, selpos, weights=weights)
>>> rmsd_expected = np.sqrt([16.    , 15.6625])
>>> np.allclose(rmsd, rmsd_expected, rtol=0)
True
>>> mdt.strc.rmsd(refpos, selpos, weights=weights, center=True)
array([1.56324982, 2.54478634])
>>> mdt.strc.rmsd(
...     refpos, selpos, xyz=True, center=True, inplace=True
... )
array([[ 0.  ,  1.  ,  4.  ],
       [ 1.  , 12.25,  0.  ]])
>>> refpos
array([[[ 1. ,  2.5,  2. ],
        [-1. , -2.5, -2. ]],

       [[ 2. , -1. ,  2. ],
        [-2. ,  1. , -2. ]]])
>>> selpos
array([[[ 1. ,  3.5,  0. ],
        [-1. , -3.5,  0. ]],

       [[ 1. ,  2.5,  2. ],
        [-1. , -2.5, -2. ]]])

Shape of refpos is (n, 3) and shape of selpos is (k, n, 3):

>>> refpos = np.array([[ 3.,  4.,  5.],
...                    [ 1., -1.,  1.]])
>>> selpos = np.array([[[ 1.,  7., -1.],
...                     [-1.,  0., -1.]],
...
...                    [[ 3.,  4.,  5.],
...                     [ 1., -1.,  1.]]])
>>> selpos - refpos
array([[[-2.,  3., -6.],
        [-2.,  1., -2.]],

       [[ 0.,  0.,  0.],
        [ 0.,  0.,  0.]]])
>>> rmsd = mdt.strc.rmsd(refpos, selpos)
>>> rmsd_expected = np.sqrt(0.5 * np.array([49 + 9, 0]))
>>> np.allclose(rmsd, rmsd_expected, rtol=0)
True
>>> mdt.strc.rmsd(refpos, selpos, xyz=True)
array([[ 4.,  5., 20.],
       [ 0.,  0.,  0.]])
>>> box = np.array([5, 3, 3, 90, 90, 90])
>>> mdt.strc.rmsd(refpos, selpos, xyz=True, box=box)
array([[4. , 0.5, 0.5],
       [0. , 0. , 0. ]])
>>> weights = np.array([0.575, 0.425])
>>> rmsd = mdt.strc.rmsd(refpos, selpos, weights=weights)
>>> rmsd_expected = np.array([4., 0.])
>>> np.allclose(rmsd, rmsd_expected, rtol=0)
True
>>> mdt.strc.rmsd(refpos, selpos, weights=weights, center=True)
array([1.56324982, 0.        ])
>>> mdt.strc.rmsd(
...    refpos, selpos, xyz=True, center=True, inplace=True
... )
array([[0., 1., 4.],
       [0., 0., 0.]])
>>> refpos
array([[ 1. ,  2.5,  2. ],
       [-1. , -2.5, -2. ]])
>>> selpos
array([[[ 1. ,  3.5,  0. ],
        [-1. , -3.5,  0. ]],

       [[ 1. ,  2.5,  2. ],
        [-1. , -2.5, -2. ]]])