@@ -298,15 +298,27 @@ def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5)
298298 """
299299 Function to calculate the moment of inertia weighted torques for a given bead.
300300
301- Input
302- -----
303- bead : the part of the molecule to be considered
304- rot_axes : the axes relative to which the forces and coordinates are located
305- frame : the frame number from the trajectory
306-
307- Output
308- ------
309- weighted_torque : the mass weighted sum of the torques in the bead
301+ This function computes torques in a rotated frame and then weights them using
302+ the moment of inertia tensor. To prevent numerical instability, it treats
303+ extremely small diagonal elements of the moment of inertia tensor as zero
304+ (since values below machine precision are effectively zero). This avoids
305+ unnecessary use of extended precision (e.g., float128).
306+
307+ Parameters
308+ ----------
309+ data_container : object
310+ Contains atomic positions and forces.
311+ bead : object
312+ The part of the molecule to be considered.
313+ rot_axes : np.ndarray
314+ The axes relative to which the forces and coordinates are located.
315+ force_partitioning : float, optional
316+ Factor to adjust force contributions, default is 0.5.
317+
318+ Returns
319+ -------
320+ np.ndarray
321+ The mass-weighted sum of the torques in the bead.
310322 """
311323
312324 torques = np .zeros ((3 ,))
@@ -336,18 +348,19 @@ def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5)
336348 moment_of_inertia = bead .moment_of_inertia ()
337349
338350 for dimension in range (3 ):
339- # Check if the moment of inertia is valid for square root calculation
340- inertia_value = moment_of_inertia [dimension , dimension ]
341-
342- if np .isclose (inertia_value , 0 ):
343- raise ValueError (
344- f"Invalid moment of inertia value: { inertia_value } . "
345- f"Cannot compute weighted torque."
351+ if np .isclose (moment_of_inertia [dimension , dimension ], 0 ):
352+ weighted_torque [dimension ] = torques [dimension ]
353+ else :
354+ if moment_of_inertia [dimension , dimension ] < 0 :
355+ raise ValueError (
356+ f"Negative value encountered for moment of inertia: "
357+ f"{ moment_of_inertia [dimension , dimension ]} "
358+ f"Cannot compute weighted torque."
359+ )
360+ weighted_torque [dimension ] = torques [dimension ] / np .sqrt (
361+ moment_of_inertia [dimension , dimension ]
346362 )
347363
348- # compute weighted torque if moment of inertia is valid
349- weighted_torque [dimension ] = torques [dimension ] / np .sqrt (inertia_value )
350-
351364 return weighted_torque
352365
353366
0 commit comments