Change of axis for moment of inertia calculation.

.github/renovate.json

@@ -7,6 +7,8 @@
     "group:recommended"
   ],
   "rangeStrategy": "pin",
+  "schedule": ["before 8am on weekdays"],
+  "labels": ["dependencies"],
   "lockFileMaintenance": {
     "enabled": true
   },

.github/workflows/mdanalysis-compatibility.yaml

@@ -18,10 +18,10 @@ jobs:
 
     steps:
       - name: Checkout repo
-        uses: actions/checkout@v5
+        uses: actions/checkout@v6
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v6.0.0
+        uses: actions/setup-python@v6.1.0
         with:
           python-version: ${{ matrix.python-version }}
 

.github/workflows/project-ci.yaml

@@ -18,10 +18,10 @@ jobs:
         python-version: ["3.11", "3.12", "3.13", "3.14"]
     steps:
       - name: Checkout repo
-        uses: actions/checkout@v5
+        uses: actions/checkout@v6
 
       - name: Set up Python ${{ matrix.python-version }}
-        uses: actions/setup-python@v6.0.0
+        uses: actions/setup-python@v6.1.0
         with:
           python-version: ${{ matrix.python-version }}
 
@@ -40,9 +40,9 @@ jobs:
     runs-on: ubuntu-latest
     timeout-minutes: 15
     steps:
-      - uses: actions/checkout@v5
+      - uses: actions/checkout@v6
       - name: Set up Python 3.14
-        uses: actions/setup-python@v6.0.0
+        uses: actions/setup-python@v6.1.0
         with:
           python-version: 3.14
       - name: Install python dependencies
@@ -56,9 +56,9 @@ jobs:
     runs-on: ubuntu-24.04
     timeout-minutes: 15
     steps:
-      - uses: actions/checkout@v5
+      - uses: actions/checkout@v6
       - name: Set up Python 3.14
-        uses: actions/setup-python@v6.0.0
+        uses: actions/setup-python@v6.1.0
         with:
           python-version: 3.14
       - name: Install python dependencies

.github/workflows/release.yaml

@@ -18,10 +18,10 @@ jobs:
     steps:
     - name: Checkout repository
       id: repo
-      uses: actions/checkout@v5
+      uses: actions/checkout@v6
 
     - name: Set up Python
-      uses: actions/setup-python@v6.0.0
+      uses: actions/setup-python@v6.1.0
       with:
         python-version: 3.14
 
@@ -46,7 +46,7 @@ jobs:
     steps:
 
       - name: checkout
-        uses: actions/checkout@v5
+        uses: actions/checkout@v6
 
       - name: Change version in repo and CITATION.cff
         run: |
@@ -61,7 +61,7 @@ jobs:
 
       - name: send PR
         id: pr_id
-        uses: peter-evans/[email protected].8
+        uses: peter-evans/create-pull-request@84ae59a2cdc2258d6fa0732dd66352dddae2a412 # v7.0.9
         with:
           commit-message: Update version to ${{ github.event.inputs.version }}
           branch: version-update
@@ -95,7 +95,7 @@ jobs:
     runs-on: ubuntu-24.04
     steps:
       - name: Checkout repository
-        uses: actions/checkout@v5
+        uses: actions/checkout@v6
         with:
           ref: main
 
@@ -126,12 +126,12 @@ jobs:
     steps:
 
     - name: checkout
-      uses: actions/checkout@v5
+      uses: actions/checkout@v6
       with:
         ref: main
 
     - name: Set up Python
-      uses: actions/setup-python@v6.0.0
+      uses: actions/setup-python@v6.1.0
       with:
         python-version: 3.14
 

CodeEntropy/levels.py

@@ -102,7 +102,6 @@ def get_matrices(
           force_matrix (np.ndarray): Accumulated force covariance matrix.
           torque_matrix (np.ndarray): Accumulated torque covariance matrix.
         """
-
         # Make beads
         list_of_beads = self.get_beads(data_container, level)
 
@@ -229,7 +228,6 @@ def get_dihedrals(self, data_container, level):
                         + str(residue - 3)
                     )
                     atom1 = data_container.select_atoms(atom_string)
-
                     atom_string = (
                         "resindex "
                         + str(residue - 3)
@@ -389,12 +387,15 @@ def get_axes(self, data_container, level, index=0):
 
         if level == "polymer":
             # for polymer use principle axis for both translation and rotation
+            # axes are going to be sorted as MDAnalysis sorts them
+            # highest eigenvalue corresponds to first axis and so on
             trans_axes = data_container.atoms.principal_axes()
             rot_axes = data_container.atoms.principal_axes()
 
         elif level == "residue":
             # Translation
             # for residues use principal axes of whole molecule for translation
+            # principal axes are going to be sorted as MDAnalysis sorts them
             trans_axes = data_container.atoms.principal_axes()
 
             # Rotation
@@ -412,6 +413,7 @@ def get_axes(self, data_container, level, index=0):
 
             if len(atom_set) == 0:
                 # if no bonds to other residues use pricipal axes of residue
+                # principal axes are going to be sorted as MDAnalysis sorts them
                 rot_axes = residue.atoms.principal_axes()
 
             else:
@@ -426,7 +428,15 @@ def get_axes(self, data_container, level, index=0):
 
         elif level == "united_atom":
             # Translation
-            # for united atoms use principal axes of residue for translation
+            # same axes as residue level rotation
+            # residue = data_container.residues()
+            # res_index_prev = res_index - 1
+            # res_index_next = res_index + 1
+            # atom_set = data_container.select_atoms(
+            #    f"(resindex {index_prev} or resindex {index_next}) "
+            #    f"and bonded resid {index}"
+            # )
+
             trans_axes = data_container.residues.principal_axes()
 
             # Rotation
@@ -666,6 +676,7 @@ def get_weighted_torques(self, data_container, bead, rot_axes, force_partitionin
             coords_rot = (
                 data_container.atoms[atom.index].position - bead.center_of_mass()
             )
+
             coords_rot = np.matmul(rot_axes, coords_rot)
             # update local forces in rotational frame
             forces_rot = np.matmul(rot_axes, data_container.atoms[atom.index].force)
@@ -681,10 +692,32 @@ def get_weighted_torques(self, data_container, bead, rot_axes, force_partitionin
             torques += torques_local
 
         # divide by moment of inertia to get weighted torques
+        # moment of inertia axes match the level rotation axes
         # moment of inertia is a 3x3 tensor
         # the weighting is done in each dimension (x,y,z) using the diagonal
         # elements of the moment of inertia tensor
-        moment_of_inertia = bead.moment_of_inertia()
+        # moment of inertia is calculated using the rotational axes
+        # axes are already sorted
+        # Ixx = sum(m_i (y_i^2+z_i^2))
+        # Iyy = sum(m_i (x_i^2+z_i^2))
+        # Izz = sum(m_i (x_i^2+y_i^2))
+
+        moment_of_inertia_diagonals = np.zeros((3,))
+        for atom in bead.atoms:
+            mass = atom.mass
+            coords_rot = (
+                data_container.atoms[atom.index].position - bead.center_of_mass()
+            )
+            coords_rot = np.matmul(rot_axes, coords_rot)
+            moment_of_inertia_diagonals[0] += mass * (
+                coords_rot[1] * coords_rot[1] + coords_rot[2] * coords_rot[2]
+            )
+            moment_of_inertia_diagonals[1] += mass * (
+                coords_rot[0] * coords_rot[0] + coords_rot[2] * coords_rot[2]
+            )
+            moment_of_inertia_diagonals[2] += mass * (
+                coords_rot[0] * coords_rot[0] + coords_rot[1] * coords_rot[1]
+            )
 
         for dimension in range(3):
             # Skip calculation if torque is already zero
@@ -693,23 +726,23 @@ def get_weighted_torques(self, data_container, bead, rot_axes, force_partitionin
                 continue
 
             # Check for zero moment of inertia
-            if np.isclose(moment_of_inertia[dimension, dimension], 0):
+            if np.isclose(moment_of_inertia_diagonals[dimension], 0):
                 raise ZeroDivisionError(
                     f"Attempted to divide by zero moment of inertia in dimension "
                     f"{dimension}."
                 )
 
             # Check for negative moment of inertia
-            if moment_of_inertia[dimension, dimension] < 0:
+            if moment_of_inertia_diagonals[dimension] < 0:
                 raise ValueError(
                     f"Negative value encountered for moment of inertia: "
-                    f"{moment_of_inertia[dimension, dimension]} "
+                    f"{moment_of_inertia_diagonals[dimension]} "
                     f"Cannot compute weighted torque."
                 )
 
             # Compute weighted torque
             weighted_torque[dimension] = torques[dimension] / np.sqrt(
-                moment_of_inertia[dimension, dimension]
+                moment_of_inertia_diagonals[dimension]
             )
 
         logger.debug(f"Weighted Torque: {weighted_torque}")

pyproject.toml

@@ -38,15 +38,15 @@ classifiers = [
 keywords = ["entropy", "macromolecular systems", "MD simulation"]
 requires-python = ">=3.11"
 dependencies = [
-    "numpy==2.3.4",
+    "numpy==2.3.5",
     "mdanalysis==2.10.0",
     "pandas==2.3.3",
     "psutil==7.1.3",
     "PyYAML==6.0.3",
     "python-json-logger==4.0.0",
     "rich==14.2.0",
     "art==6.5",
-    "waterEntropy==1.2.2",
+    "waterEntropy==1.2.3",
     "requests==2.32.5",
 ]