[WIP] Fix for LPRMSD after mdtraj merge

[schwancr]

LPRMSD no longer worked when we moved mdtraj into msmbuilder.

It's still necessary because the altindices behavior is not available in mdtraj (at least to my knowledge).

I attempted to fix it by just appending the old RMSD.TheoData class into lprmsd.py

[schwancr]

.... But, unfortunately that wasn't good enough because I'm segfaulting all over the place. Not sure why

[schwancr]

cc: @leeping @cxh

[leeping]

I tried to fix this before, but I only spent ~2 hours on it and didn't succeed. Can I take a look @ it when I come back?

[schwancr]

You ran into the same seg fault?

On Tue, Jun 17, 2014 at 8:39 PM, Lee-Ping notifications@github.com wrote:

I tried to fix this before, but I only spent ~2 hours on it and didn't succeed. Can I take a look @ it when I come back?

— Reply to this email directly or view it on GitHub https://github.com/SimTk/msmbuilder/pull/427#issuecomment-46392884.

[rmcgibbo]

Using mdtraj.lprmsd, you can trivially get the former alt_indices behavior by using superpose=True and then calculating the root mean squared deviation manually.

i.e

md.lprmsd(target, reference, index, atom_indices=atom_indices,
superpose=True)
deviation = target.xyz[:, alt_indices, :] - reference.xyz[index,
alt_indices, :]
value = np.sqrt(np.mean(deviation**2, axis=1))

On Tue, Jun 17, 2014 at 10:48 PM, Christian Schwantes < notifications@github.com> wrote:

You ran into the same seg fault?

On Tue, Jun 17, 2014 at 8:39 PM, Lee-Ping notifications@github.com wrote:

I tried to fix this before, but I only spent ~2 hours on it and didn't succeed. Can I take a look @ it when I come back?

— Reply to this email directly or view it on GitHub https://github.com/SimTk/msmbuilder/pull/427#issuecomment-46392884.

— Reply to this email directly or view it on GitHub https://github.com/SimTk/msmbuilder/pull/427#issuecomment-46398399.

[schwancr]

Ok that's an easy enough way to implement it.

@leeping @mlawrenz @cxhernandez and anyone else, what features in LPRMSD do you need. There are multiple "align to density," or "align to moments" functions. Do you need that?

[leeping]

Hi Christian,

Thanks; I don't need those other features.

When I installed msmbuilder 2.8, LPRMSD was broken and it wasn't using the MDTraj implementation, so I reverted to msmbuilder 2.7. Is the situation different now?

[schwancr]

I'm working on translating it now. It slipped through the cracks during the transfer to mdtraj. It will have this functionality:

1) Align to one set of atoms and compute the distance with another set 2) Align / calculate distance with permutable atoms

[leeping]

Cool. Can we implement code to get the actual aligned structures using LPRMSD? That was a feature we had before that I don't want to lose - I actually use it.

I think Trajectory.superpose() in MDTraj currently provides the methods for getting aligned structures using RMSD, but there are other metrics (e.g. LPRMSD) that have the ability to provide aligned structures.

What do you think is the best way to do it?

[schwancr]

So, as robert pointed out, that is in the mdtraj.lprmsd function. If you pass superpose=True then it will modify the contents of the trajectory to be aligned to the reference.

[leeping]

Oh, cool. :) This is something I could add to SaveStructures at some point.

[kyleabeauchamp]

Lee-Ping is right in pointing out that we want more alignment features to be supported and added to MDTraj. Ideally, RMSD, LPRMSD, alignment.py and rotation.py could be cleaned up to make a single API--right now it's a bit of a mess IMHO.

[leeping]

It is a bit messy but I think this could be due to the small number of metrics that can do alignments. RMSD and LPRMSD are the two I can think of - what other metrics are there?

[kyleabeauchamp]

In MDTraj, there is no concept of metric.

[schwancr]

@leeping do you have any data from the old version that I could use to test the implementation?

I need a simple test for making sure the permute groups are working and that the alt indices are working.

[schwancr]

nevermind, I came up with some simple examples. There seems to be a bug right now. I think it's in mdtraj though.

[schwancr]

There's no bug, I assumed the default for permute_groups was different than it is

[schwancr]

This test fails: https://gist.github.com/schwancr/648d2955b2315835234f

The distance is correct at zero, but the superpose isn't working because the rotation matrix cannot be converged. The warning is printed here: https://github.com/rmcgibbo/mdtraj/blob/master/MDTraj/rmsd/src/theobald_rmsd.c#L267-L319

This region is computing a rotation matrix for a given quaternion, which I think can be done many ways (though I don't know where these lines of code are coming from), since in @leeping's initial code, there was a series of nested if statements checking for this "non-convergence" and then solving the problem in a different way. See here: https://github.com/SimTk/msmbuilder/blob/998dfdaf22946dd81a8617442adbd155498a6e93/Extras/LPRMSD/src/qcprot.c#L272-L319

@rmcgibbo or @leeping does any of this sound familiar? I don't know the math to be able to figure this out at the moment.

[schwancr]

Ok, so I found the reference code and put it in a new branch in mdtraj (rmcgibbo/mdtraj#503), but it still doesn't fix this test.

[leeping]

Hi Christian,

I'm not sure if even the reference code would pass the test you proposed. There's a lot of degeneracy in the rotation matrix for aligning a system of linear atoms that isn't there for the calculations that I ran.

Thanks,

• Lee-Ping

[schwancr]

Yea, I think you're right, it seems to be working otherwise though.

[leeping]

Do you think we should fix the case of linear atoms or come up with a different test case (e.g. a cube of atoms)?

[schwancr]

I think a different test case is easiest.

Right now I fixed the issues with msmbuilder.metrics.Positions which used LPRMSD, but I haven't replaced the LPRMSD metric yet