rmsd

mdtools.structure.rmsd(refpos, selpos, weights=None, center=False, inplace=False, xyz=False, box=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 the refpos and selpos by their (weighted) center, respectively, before calculating the RMSD.
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 must be orthogonal and 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. Works currently only for orthogonal boxes.

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^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^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^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.])
>>> 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. ]])
>>> 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. ]])
>>> 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. ]]])