From 67c55a0c6383a0ee0694ae33341a9cd75a3071da Mon Sep 17 00:00:00 2001 From: harryswift01 Date: Wed, 12 Mar 2025 13:18:08 +0000 Subject: [PATCH 1/3] Revert weighted torque calculation and improve documentation: - Restored the original method for computing weighted torque - Clarified handling of small moment of inertia values in the docstring, explaining why values below machine precision are treated as zero - No functional changes beyond reverting to the previous calculation method --- CodeEntropy/GeometricFunctions.py | 51 +++++++++++++++++++------------ 1 file changed, 32 insertions(+), 19 deletions(-) diff --git a/CodeEntropy/GeometricFunctions.py b/CodeEntropy/GeometricFunctions.py index fb7a0b9..7cd30aa 100644 --- a/CodeEntropy/GeometricFunctions.py +++ b/CodeEntropy/GeometricFunctions.py @@ -298,15 +298,27 @@ def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5) """ Function to calculate the moment of inertia weighted torques for a given bead. - Input - ----- - bead : the part of the molecule to be considered - rot_axes : the axes relative to which the forces and coordinates are located - frame : the frame number from the trajectory - - Output - ------ - weighted_torque : the mass weighted sum of the torques in the bead + This function computes torques in a rotated frame and then weights them using + the moment of inertia tensor. To prevent numerical instability, it treats + extremely small diagonal elements of the moment of inertia tensor as zero + (since values below machine precision are effectively zero). This avoids + unnecessary use of extended precision (e.g., float128). + + Parameters + ---------- + data_container : object + Contains atomic positions and forces. + bead : object + The part of the molecule to be considered. + rot_axes : np.ndarray + The axes relative to which the forces and coordinates are located. + force_partitioning : float, optional + Factor to adjust force contributions, default is 0.5. + + Returns + ------- + np.ndarray + The mass-weighted sum of the torques in the bead. """ torques = np.zeros((3,)) @@ -336,18 +348,19 @@ def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5) moment_of_inertia = bead.moment_of_inertia() for dimension in range(3): - # Check if the moment of inertia is valid for square root calculation - inertia_value = moment_of_inertia[dimension, dimension] - - if np.isclose(inertia_value, 0): - raise ValueError( - f"Invalid moment of inertia value: {inertia_value}. " - f"Cannot compute weighted torque." + if np.isclose(moment_of_inertia[dimension, dimension], 0): + weighted_torque[dimension] = torques[dimension] + else: + if moment_of_inertia[dimension, dimension] < 0: + raise ValueError( + f"Negative value encountered for moment of inertia: " + f"{moment_of_inertia[dimension, dimension]} " + f"Cannot compute weighted torque." + ) + weighted_torque[dimension] = torques[dimension] / np.sqrt( + moment_of_inertia[dimension, dimension] ) - # compute weighted torque if moment of inertia is valid - weighted_torque[dimension] = torques[dimension] / np.sqrt(inertia_value) - return weighted_torque From efc42645d16c1a65815167df086ba415bab8e9d2 Mon Sep 17 00:00:00 2001 From: harryswift01 Date: Wed, 12 Mar 2025 13:33:43 +0000 Subject: [PATCH 2/3] Optimize weighted torque calculation to reduce crashes: - Added a check to skip the calculation if the torque is already zero, highlighted in issue #55, reducing unnecessary calculations - Retained the divide-by-zero check to ensure numerical stability when computing weighted torques - Improved docstring to reflect the updated logic --- CodeEntropy/GeometricFunctions.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/CodeEntropy/GeometricFunctions.py b/CodeEntropy/GeometricFunctions.py index 7cd30aa..346ea7a 100644 --- a/CodeEntropy/GeometricFunctions.py +++ b/CodeEntropy/GeometricFunctions.py @@ -304,6 +304,9 @@ def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5) (since values below machine precision are effectively zero). This avoids unnecessary use of extended precision (e.g., float128). + Additionally, if the computed torque is already zero, the function skips + the division step, reducing unnecessary computations and potential errors. + Parameters ---------- data_container : object @@ -348,6 +351,11 @@ def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5) moment_of_inertia = bead.moment_of_inertia() for dimension in range(3): + # Skip calculation if torque is already zero + if np.isclose(torques[dimension], 0): + weighted_torque[dimension] = 0 + continue + if np.isclose(moment_of_inertia[dimension, dimension], 0): weighted_torque[dimension] = torques[dimension] else: From 456b9e33adc56f8c2fdf194aca7be41ed32cfce8 Mon Sep 17 00:00:00 2001 From: harryswift01 Date: Thu, 13 Mar 2025 10:54:10 +0000 Subject: [PATCH 3/3] Prevent divide-by-zero error by raising ZeroDivisionError when moment of inertia is zero instead of proceeding with unweighted torque calculation. --- CodeEntropy/GeometricFunctions.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/CodeEntropy/GeometricFunctions.py b/CodeEntropy/GeometricFunctions.py index 346ea7a..67cc28a 100644 --- a/CodeEntropy/GeometricFunctions.py +++ b/CodeEntropy/GeometricFunctions.py @@ -356,19 +356,26 @@ def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5) weighted_torque[dimension] = 0 continue + # Check for zero moment of inertia if np.isclose(moment_of_inertia[dimension, dimension], 0): - weighted_torque[dimension] = torques[dimension] - else: - if moment_of_inertia[dimension, dimension] < 0: - raise ValueError( - f"Negative value encountered for moment of inertia: " - f"{moment_of_inertia[dimension, dimension]} " - f"Cannot compute weighted torque." - ) - weighted_torque[dimension] = torques[dimension] / np.sqrt( - moment_of_inertia[dimension, dimension] + 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: + raise ValueError( + f"Negative value encountered for moment of inertia: " + f"{moment_of_inertia[dimension, dimension]} " + f"Cannot compute weighted torque." ) + # Compute weighted torque + weighted_torque[dimension] = torques[dimension] / np.sqrt( + moment_of_inertia[dimension, dimension] + ) + return weighted_torque