Skip to content

Conversation

@jkalayan
Copy link
Collaborator

@jkalayan jkalayan commented Feb 2, 2026

Summary

CodeEntropy has been developed from code that calculated protein entropies, this PR aligns with the entropy methods used in the 2019 pure liquids paper, which compared calculated entropies with experimental values of ~50 pure liquids. To ensure CodeEntropy is comparative with these experimental results, this PR gives the options to use the custom axes and combined force-torque covariance matrix methods used in the 2019 paper.

The 2019 paper used the bonded axes to calculate Srovib at the united-atom (UA) level only, the bonded axes are dependent on what a united atom is bonded to. It also used a variant of the principal axes, where the heavy atom positions and combined united atom masses are used to calculate principal axes. Lastly, the 2019 paper used the combined force-torque covariance matrix to calculate trans- and ro-vibrational entropies at the molecule/highest level only.

This differs from the current CodeEntropy, which uses separate force/torque covariance matrices for all levels and uses the MDA principal_axes() function for axes across all length-scales/levels to rotate forces.

Changes

Update Covariance Matrices :

  • Option to use a combined (6,6) force-torque matrix for the highest length scale level, use the combined_forcetorque argument in the config to activate this combined force-torque matrix, instead of using the default case of separate force/torque matrices across all levels.
  • Resolve small bug that didn't populate the first force and torque submatrices
  • If this option is activated, the Svib terms outputted are FTmat-Transvibrational and FTmat-Rovibrational, otherwise the default Transvibrational and Rovibrational (using separate force/trque matrices) are displayed on the splashscreen

Custom Axes:

  • Option for custom axes for Svib calculations, use the customised_axes argument in the config to activate custom axes, instead of the default principal axes
  • Add an AxesManager class to calculate axes to rotate forces around
  • Calculate axes from bonds for UA-rovib and with a tweaked principal axes for all other Svib
  • Use of custom axes to get the moment of inertia and torques does give more zero and negative values, so ValueError and warnings have been switched off
  • Update default principal axes to use the MDA moment_of_inertia() function, gives similar but slightly different Svib values compared to current main

Dimensions in Universe :

  • Include dimensions in universe to account for periodic boundary conditions for vector calculations used to calculate custom bonded axes

Impact

  • Forms a basis for testing CodeEntropy against experimental values, ensures future code changes will either match or improve calculated entropies compared with experiments.
  • Aligns entropy results with the 2019 paper, for the seven molecules tested so far (WM= residue/whole molecule level). Gives <0.25 J/K/mol difference with Jon's C++ code version 37, when run for 1 frame for all 500 molecules in the simulation (same trajectories as 2019 paper, GAFF)
Screenshot 2026-02-02 at 11 12 34

@jkalayan jkalayan self-assigned this Feb 2, 2026
heavy_bonded, light_bonded = self.find_bonded_atoms(atom.index, system)
UA = atom + light_bonded
# UA_all = atom + heavy_bonded + light_bonded

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are these commented out lines to be removed?

dimensions,
)
# case4, not used in Jon's code, use case5 instead
# if len(heavy_bonded) == 2:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should these be removed too?

eigenvalues, _eigenvectors = np.linalg.eig(moment_of_inertia)
# sort eigenvalues of moi tensor by largest to smallest magnitude
order = abs(eigenvalues).argsort()[::-1] # decending order
# principal_axes = principal_axes[order] # PI already ordered correctly
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If principal axes is going to be a user set option, would this being commented out cause issues?

Copy link
Member

@jimboid jimboid left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like the changes are good, I placed some questions on the PR source asking about some commented lines of code. Should we keep them or remove them, there was 1 line for principal axes that is commented out, if we are keeping those as a user requested option then is that one commented out line going to cause issues?

We also know from the numerical tests that this is much more accurate than the code it is going to fix.

@jkalayan
Copy link
Collaborator Author

jkalayan commented Feb 2, 2026

Looks like the changes are good, I placed some questions on the PR source asking about some commented lines of code. Should we keep them or remove them, there was 1 line for principal axes that is commented out, if we are keeping those as a user requested option then is that one commented out line going to cause issues?

We also know from the numerical tests that this is much more accurate than the code it is going to fix.

Thanks James, removed redundant commented lines and functions that were only used to test the custom axes method

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants