Residue rotation and UA translation axes for terminal residues

[ioanaapapa]

Summary

Terminal residue axes are no longer dependent on atom indices. The edge atom (i.e. heavy atom bonded to neighbour residue) is chosen alongside the average position of heavy atoms bonded to the edge atom and average position of all other heavy in the residue as the three points used for determining the axes. The centre of rotation (at res level)/translation (at UA level) and axes are then computed from the three points similar to the case of non-terminal residues.

Changes

Terminal residue treatment was changed

• Residue rotation and UA translation axes for terminal residues no longer use atom indices to emulate the backbone for non-terminal residues.
• The edge atom (E1 - red), defined as the heavy atom bonded to a neighbour residue is identified. All heavy atoms bonded to the edge atom (purple) are identified and their average position (C) is used as the third point for obtaining the centre of rotation. All other heavy atoms in the residue (green) are identified and used as the equivalent to the second edge(E2).

• The centre of rotation (O - blue) is the point where the perpendicular from the C (purple) meets the E1 (red)-E2 (green) vector.
• x-axis is along the centre of rotation (O) - first edge(E1) vector
• y axis is along the centre of rotation (O) - average position of bonded atoms (C) vector
• z-axis is perpendicular to the two

Edge cases have been introduced

• If there are no other heavy atoms in a terminal residue, other than the edge atom and those bonded to the edge atom, the centre of rotation is the edge atom and axes are defined as: x-axis is along the edge-bonded vector, y axis is an arbitrary perpendicular axis and z-axis is perpendicular to the two. (This uses the get_custom_axes function and is the same as case2 in get_bonded_axes for UAs)
• If the two edge atoms in a non-terminal residue are bonded (there is no backbone atom between them), the centre of rotation is located at their average position and axes are defined as: x-axis is along the centre-of-rotation-first edge vector, y-axis is arbitrary perpendicular and z-axis is perpendicular to the two.

Tests for terminal residues and edge cases

• Unit tests for terminal residues have been updated to reflect changes.
• New tests for the edge cases covered have been introduced

Documentation

• Vibrational entropy in theory section of the documentation was modified to reflect changes.

get_residue_axes and get_UA_axes structure change

• New functions get_terminal_axes and get_non_terminal_axes are introduced and called in get_residue_axes and get_UA_axes to improve code flow.

Impact

• This closes issue [Feature]: Axes for terminal residues #399.
• Terminal residue treatment is no longer dependent on atom indexing, which is arbitrary. The same result would be obtained using topologies of the same molecule with different indexes.
• Residue rovibrational results are now slightly higher compared to total entropy results obtained using the index-based backbone. No signficant difference observed for UA transvibrational results. (only treatment of terminal residues differs in table below)

[harryswift01]

Thank you for sending these changes through! I think this is a really good improvement, and it’s great to see the terminal residue handling becoming more robust and independent of atom indexing. I’ve just left a few comments around some small technical and maintainability points. If you’re able to have a look through them and let me know what you think, that would be great. Overall, I think this is looking really good and is very close to being ready to go in!

[ioanaapapa]

@harryswift01 Thank you, Harry! I will make the changes.

[harryswift01]

Thanks, @harryswift01! Is your last comment about the test mocking the wrong part of get_custom_axes what you were referring to as an issue in multiple tests or is there a second issue, as well?

Hi @ioanaapapa, thanks for double checking! Yes, this would be the same issue, not a second one, I just spotted it in a few other tests too. If you're able to apply the same fix there as well, that'd be great.

[ioanaapapa]

Thank you! I'll change the tests

[harryswift01]

Thanks again for making these changes @ioanaapapa it's looking really good now! Just a couple of small things I noticed, mostly around the docstrings: there are a few repeated lines in there that would be great if we tidied up, I've left comments pointing to each one. I also spotted one thing in get_residue_custom_axes it might be worth adding a guard there to catch zero-length axis vectors, similar to what get_custom_axes already does, just so we don't end up with any NaN values slipping through.

Other than that, this is very close, just a few minor polishing points and I think we'll be good to merge.

[harryswift01]

These two lines also repeat here as well, lines 347 and 348

[harryswift01]

Whilst I was having a look, I saw that get_custom_axes guards against zero-length axis vectors and raises a clear ValueError if they end up with zero length, but get_residue_custom_axes doesn't have the same check. I think it would be nice to put that raise check in here too, in case any NaN values creep in.

[ioanaapapa]

@harryswift01 I've tidied up the docstrings and added a ValueError raise in get_residue_custom_axes, alongside a unit test for that. Hope it's okay now!

[harryswift01]

Hi @ioanaapapa, thank you for all your work on this! This is all looking good from a technical side now and I can't see any reason why this can't be merged in. Just before we do merge this in, as this is a change in the numerics/science, if @skfegan or @jkalayan or @jimboid could have also have a look at this as well just to make sure there isn't anything I've missed from the science side, then after their approval I'd be happy to see this merged in, this is all great work, thank you again!

[jimboid]

Improvements to the way terminal residues are treated will help to fix some edge case issues reported on the issue tracker. This has been discussed at length at the all hands meetings with Richard, happy to see this merged @ioanaapapa.