RMSD throws divide-by-zero error if a position is exactly 0.

[jdetle]

Expected behaviour

A proper rmsd should be outputted from two sets of XYZ coordinates.

Actual behaviour

If any entry in either of the two arrays has XYZ coordinates are all truly zero, the code FastCalcRMSDandRotation code will throw a divide by zero error.

Code to reproduce the behaviour

import MDAnalysis.analysis.rms as rms
a_mda = np.ones((100, 3))
b_mda = np.zeros((100, 3))

for i in range(100):
    a_mda[i, :] += i
    b_mda[i, :] += i

b_mda[0, :] = 0
x_mda = rms.rmsd(a_mda, b_mda, superposition=True)

Changing b_mda[0, :] = 0to 0.0001 throws no error.

Currently version of MDAnalysis:

Develop branch

There are a few hacky ways I can imagine to go about solving this problem, I look forward to hearing your thoughts on the matter.

[orbeckst]

Although I can reproduce the error with your example, this does not appear to be true in general: The following code works for me:

import numpy as np
import MDAnalysis.analysis.rms as rms

a = 10 * np.random.randn(10, 3)
a[0, :] = [0, 0, 0]
b = 10 * np.random.randn(10, 3)

x = rms.rmsd(a, b, superposition=True)
print(x)

and even

a = np.zeros((10, 3))
b = 10 * np.random.randn(10, 3)
x = rms.rmsd(a, b, superposition=True)
print(x)

works for me.

Not really sure where the exact problem lies... perhaps some actual code digging is order.

However, can you check that my examples work for you?

[jdetle]

Will do, as usual I am busier than I thought I'd be but I'll try to get to this soon.

[orbeckst]

This is a weird bug.

import numpy as np
from MDAnalysis.analysis import  rms
a = np.ones((100, 3))
b = np.zeros((100, 3))
for i in range(100):
    a[i, :] += i
    b[i, :] += i

assert np.allclose(b[0], [0, 0, 0])

try:
     rms.rmsd(a, b, superposition=True)
except ZeroDivisionError as err:
     print("It fracking failed, just as @jdetle said: {}".format(err))

def check_rms(a, b):
   for ix in range(len(b)):
       q = b.copy()
       q[0, :] = q[ix, :]
       q[ix, :] = 0.
       try:
           rr = rms.rmsd(a, q, superposition=True)
       except ZeroDivisionError as err:
           print("a, q: [0,0,0] @ index {}: {}".format(ix, err))
       try:
           rr = rms.rmsd(q, a, superposition=True)
       except ZeroDivisionError as err:
           print("q, a: [0,0,0] @ index {}: {}".format(ix, err))

check_rms(a, b)

gives

It fracking failed, just as @jdetle said: float division
a, q: [0,0,0] @ index 0: float division
q, a: [0,0,0] @ index 0: float division

so we only get this weird error with [0,0,0] in the first position and when using these highly specific input data a and b (because using general random numbers as in https://github.com/MDAnalysis/mdanalysis/issues/945#issuecomment-240806096 does not trigger the bug).

If I change a:

a *= 0.5

or (!)

a[0, :] = 0

then no errors.

Apparently one only need a very small change:

check_rms(a+1e-13, b)

passes but

check_rms(a+1e-14, b)

fails.

If anyone has any ideas... I would have liked to close this issue but it remains as a potentially obscure bug, but nevertheless a bug.

[kain88-de]

Do you know where the exception is thrown?

[orbeckst]

In [127]: rms.rmsd(b, a, superposition=True)
---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-127-79653630aa63> in <module>()
----> 1 rms.rmsd(b, a, superposition=True)

~/anaconda3/envs/mda/lib/python3.6/site-packages/MDAnalysis/analysis/rms.py in rmsd(a, b, weights, center, superposition)
    222
    223     if superposition:
--> 224         return qcp.CalcRMSDRotationalMatrix(a, b, N, None, weights)
    225     else:
    226         if weights is not None:

MDAnalysis/lib/qcprot.pyx in MDAnalysis.lib.qcprot.CalcRMSDRotationalMatrix (MDAnalysis/lib/qcprot.c:2764)()

MDAnalysis/lib/qcprot.pyx in MDAnalysis.lib.qcprot.FastCalcRMSDAndRotation (MDAnalysis/lib/qcprot.c:3625)()

ZeroDivisionError: float division

[orbeckst]

Possibly related https://stackoverflow.com/q/63736861

[AMKCam]

Was this issue every fixed? I'm also having the same problem for the RMSD analysis.

[laurenalex77]

I've also been getting this error with the following code:

u.trajectory[0]      
ref.trajectory[0]

aligned = mda.analysis.rms.RMSD(u, ref,
                           select='backbone',
                           weights='mass', verbose=True)
R = aligned.run()

[laurenalex77]

The issue in my case was with using certain selection commands on a dms file - it seems that the dms file format doesn't have the information that MDA uses to classify some of the selection commands (ie. selection commands such as 'backbone' and 'name CA' don't work but more general ones such as 'protein' do work).

I'm not sure if there's a better workaround, but ultimately I switched to a .gro topology file for the input to run backbone RMSD. It might also be helpful to add to the error message to show when a selection command fails due to lack of info from the input file.