Principal inertia axes calculation for rigid bodies is buggy

[tanmoy7989]

The function IMP.core.get_initial_reference_frame() calculates a body-centric reference frame for the set of particles given in its argument, where the three axes of the frame point along the principal axes of inertia of the specified set of particles. However, it fails to calculate this reference frame consistently across different structures which are affine transforms of one another.

Minimum working example

System description

• before_alignment.pdb is a binary complex with chains A and B.
• reference.pdb is an slightly different conformation of this complex.
• after_alignment.pdb is a re-oriented conformation of before_alignment.pdb where the chain A has been aligned with chain A of reference.pdb. Natually, chain B in this structure is also very close (though not exact) to chain B of rererence.pdb.
• view1.png and view2.png visually compares the before and after alignment structures (grey and green respectively) with the reference (red). view2.png shows the alignment of chain A.

In each case, IMP.core.get_initial_reference_frame() is used to calculate a body centric reference frame, which is reported as a IMP.algebra.ReferenceFrame3D object. (see here) From this reference frame, the unit vectors along the three axes are extracted and reported.

Expectation The unit vectors for chain A (which is rigid body 0) should be exactly equal (in all their components) for reference.pdb and after_alignment.pdb, since coordinates of chain A was used to perform the alignment. Unit vectors for chain B (which is rigid body 1) between these two pdb files should be very close.

Output Run the python script (mwe.py). It produces:

Body centric axes before alignment:
rigid body 0, i.e. chain A: ((0.841722, -0.369063, -0.394078), (0.394398, 0.918763, -0.0180361), (0.36872, -0.140242, 0.9189))
rigid body 1, i.e. chain B: ((0.786562, 0.614593, 0.059965), (-0.451695, 0.638846, -0.622773), (-0.421061, 0.462764, 0.780101))

Body centric axes after alignment:
rigid body 0, i.e. chain A: ((-0.872915, -0.373216, 0.314213), (0.441882, -0.877798, 0.184962), (0.206785, 0.300301, 0.93116))
rigid body 1, i.e. chain B: ((0.701693, 0.0802605, 0.707944), (-0.527335, -0.609667, 0.591797), (0.479108, -0.788584, -0.385475))

Body centric axes of reference:
rigid body 0, i.e. chain A: ((0.863825, 0.387039, -0.322502), (0.458039, -0.869921, 0.182857), (-0.209779, -0.305675, -0.928739))
rigid body 1, i.e. chain B: ((0.732408, 0.167787, 0.659868), (0.517754, -0.766637, -0.379736), (0.442164, 0.619771, -0.648363))

There are differences in sign between the components (although the magnitudes are somewhat close). I have found other examples where the magnitudes are off too by small amounts for synthetic examples like this one where I expect 0 difference in magnitude between the reference chain and the aligned chain.

bcc_bug.tar.gz

[tanmoy7989]

The problem most probably stems from the fact that in https://github.com/salilab/imp/blob/develop/modules/core/src/rigid_bodies.cpp#L1017-L1020, the real eigenvectors are never ordered by eigenvalue. Typically in eigendecomposition problems like this, the principal axes are ordered largest to smallest, and for rigid body kinematics, the eigenvector corresponding to the largest eigenvalue is taken as the X axis.

[sethaxen]

By construction, the inertia tensor is symmetric positive semidefinite, so that all eigenvalues are non-negative. I'm guessing Eigen uses LAPACK to do eigendecomposition, in which case the eigenvalues are ordered least to greatest. Ordering the eigenvalues from greatest to least is the convention of SVD, and for symmetric matrices, SVD and eigendecomposition only differ by the sorting convention of the eigenvalues, so you could try using BDCSVD here, but I'm not certain if it internally uses symmetric eigendecomposition or not.

Looks like IMP is use using a standard eigen-solver here, while it should be using a SelfAdjointEigenSolver, which may be much more efficient and also void the need to call .real().

[tanmoy7989]

@sethaxen , yep, I'm trying to switch to SelfAdjointSolver and test before making the PR. Incidentally, turns out that the default solver IMP currently uses which is Eigen::EigenSolver doesn't really do a consistent ordering. I observed least->greatest->middle type of arrangements too. Perhaps because it produces both real and imaginary eigenvalues and you can't really compare the usual "greater than" between a real and an imaginary number, the Eigen folks avoid doing it altogether instead of doing it shoddily. This is of course based on a very superficial read through the Eigen::EigenSolver module.

On the other hand SelfAdjointSolver guarantees ascending order of the eigenvalues and eigenvectors are produced in that order too. So hopefully, this will just be a single line change in core/rigid_bodies.cpp. Fingers crossed :)

[tanmoy7989]

closed by #1058