diff --git a/CodeEntropy/axes.py b/CodeEntropy/axes.py
new file mode 100644
index 0000000..d698b48
--- /dev/null
+++ b/CodeEntropy/axes.py
@@ -0,0 +1,507 @@
+import logging
+
+import numpy as np
+from MDAnalysis.lib.mdamath import make_whole
+
+logger = logging.getLogger(__name__)
+
+
+class AxesManager:
+    """
+    Manages the structural and dynamic levels involved in entropy calculations. This
+    includes selecting relevant levels, computing axes for translation and rotation,
+    and handling bead-based representations of molecular systems. Provides utility
+    methods to extract averaged positions, convert coordinates to spherical systems,
+    compute weighted forces and torques, and manipulate matrices used in entropy
+    analysis.
+    """
+
+    def __init__(self):
+        """
+        Initializes the AxesManager with placeholders for level-related data,
+        including translational and rotational axes, number of beads, and a
+        general-purpose data container.
+        """
+        self.data_container = None
+        self._levels = None
+        self._trans_axes = None
+        self._rot_axes = None
+        self._number_of_beads = None
+
+    def get_residue_axes(self, data_container, index):
+        """
+        The translational and rotational axes at the residue level.
+
+        Args:
+          data_container (MDAnalysis.Universe): the molecule and trajectory data
+          index (int): residue index
+
+        Returns:
+          trans_axes : translational axes (3,3)
+          rot_axes : rotational axes (3,3)
+          center: center of mass (3,)
+          moment_of_inertia: moment of inertia (3,)
+        """
+        # TODO refine selection so that it will work for branched polymers
+        index_prev = index - 1
+        index_next = index + 1
+        atom_set = data_container.select_atoms(
+            f"(resindex {index_prev} or resindex {index_next}) "
+            f"and bonded resid {index}"
+        )
+        residue = data_container.select_atoms(f"resindex {index}")
+        center = residue.atoms.center_of_mass(unwrap=True)
+
+        if len(atom_set) == 0:
+            # No bonds to other residues
+            # Use a custom principal axes, from a MOI tensor
+            # that uses positions of heavy atoms only, but including masses
+            # of heavy atom + bonded hydrogens
+            UAs = residue.select_atoms("mass 2 to 999")
+            UA_masses = self.get_UA_masses(residue)
+            moment_of_inertia_tensor = self.get_moment_of_inertia_tensor(
+                center, UAs.positions, UA_masses, data_container.dimensions[:3]
+            )
+            rot_axes, moment_of_inertia = self.get_custom_principal_axes(
+                moment_of_inertia_tensor
+            )
+            trans_axes = (
+                rot_axes  # set trans axes to same as rot axes as per Jon's code
+            )
+        else:
+            # if bonded to other residues, use default axes and MOI
+            make_whole(data_container.atoms)
+            trans_axes = data_container.atoms.principal_axes()
+            rot_axes, moment_of_inertia = self.get_vanilla_axes(residue)
+            center = residue.center_of_mass(unwrap=True)
+
+        return trans_axes, rot_axes, center, moment_of_inertia
+
+    def get_UA_axes(self, data_container, index):
+        """
+        The translational and rotational axes at the united-atom level.
+
+        Args:
+          data_container (MDAnalysis.Universe): the molecule and trajectory data
+          index (int): residue index
+
+        Returns:
+          trans_axes : translational axes (3,3)
+          rot_axes : rotational axes (3,3)
+          center: center of mass (3,)
+          moment_of_inertia: moment of inertia (3,)
+        """
+
+        index = int(index)
+
+        # use the same customPI trans axes as the residue level
+        UAs = data_container.select_atoms("mass 2 to 999")
+        UA_masses = self.get_UA_masses(data_container.atoms)
+        center = data_container.atoms.center_of_mass(unwrap=True)
+        moment_of_inertia_tensor = self.get_moment_of_inertia_tensor(
+            center, UAs.positions, UA_masses, data_container.dimensions[:3]
+        )
+        trans_axes, _moment_of_inertia = self.get_custom_principal_axes(
+            moment_of_inertia_tensor
+        )
+
+        # look for heavy atoms in residue of interest
+        heavy_atoms = data_container.select_atoms("prop mass > 1.1")
+        heavy_atom_indices = []
+        for atom in heavy_atoms:
+            heavy_atom_indices.append(atom.index)
+        # we find the nth heavy atom
+        # where n is the bead index
+        heavy_atom_index = heavy_atom_indices[index]
+        heavy_atom = data_container.select_atoms(f"index {heavy_atom_index}")
+
+        center = heavy_atom.positions[0]
+        rot_axes, moment_of_inertia = self.get_bonded_axes(
+            data_container, heavy_atom[0], data_container.dimensions[:3]
+        )
+
+        logger.debug(f"Translational Axes: {trans_axes}")
+        logger.debug(f"Rotational Axes: {rot_axes}")
+        logger.debug(f"Center: {center}")
+        logger.debug(f"Moment of Inertia: {moment_of_inertia}")
+
+        return trans_axes, rot_axes, center, moment_of_inertia
+
+    def get_bonded_axes(self, system, atom, dimensions):
+        """
+        For a given heavy atom, use its bonded atoms to get the axes
+        for rotating forces around. Few cases for choosing united atom axes,
+        which are dependent on the bonds to the atom:
+
+        ::
+
+            X -- H = bonded to zero or more light atom/s (case1)
+
+            X -- R = bonded to one heavy atom (case2)
+
+            R -- X -- H = bonded to one heavy and at least one light atom (case3)
+
+            R1 -- X -- R2 = bonded to two heavy atoms (case4)
+
+            R1 -- X -- R2 = bonded to more than two heavy atoms (case5)
+                  |
+                  R3
+
+        Note that axis2 is calculated by taking the cross product between axis1 and
+        the vector chosen for each case, dependent on bonding:
+
+        - case1: if all the bonded atoms are hydrogens, use the principal axes.
+
+        - case2: use XR vector as axis1, arbitrary axis2.
+
+        - case3: use XR vector as axis1, vector XH to calculate axis2
+
+        - case4: use vector XR1 as axis1, and XR2 to calculate axis2
+
+        - case5: get the sum of all XR normalised vectors as axis1, then use vector
+        R1R2 to calculate axis2
+
+        axis3 is always the cross product of axis1 and axis2.
+
+        Args:
+            system: mdanalysis instance of all atoms in current frame
+            atom: mdanalysis instance of a heavy atom
+            dimensions: dimensions of the simulation box (3,)
+
+        Returns:
+            custom_axes: custom axes for the UA, (3,3) array
+            custom_moment_of_inertia
+        """
+        # check atom is a heavy atom
+        if not atom.mass > 1.1:
+            return None
+        # set default values
+        custom_moment_of_inertia = None
+        custom_axes = None
+
+        # find the heavy bonded atoms and light bonded atoms
+        heavy_bonded, light_bonded = self.find_bonded_atoms(atom.index, system)
+        UA = atom + light_bonded
+        UA_all = atom + heavy_bonded + light_bonded
+
+        # now find which atoms to select to find the axes for rotating forces:
+        # case1
+        if len(heavy_bonded) == 0:
+            custom_axes, custom_moment_of_inertia = self.get_vanilla_axes(UA_all)
+        # case2
+        if len(heavy_bonded) == 1 and len(light_bonded) == 0:
+            custom_axes = self.get_custom_axes(
+                atom.position, [heavy_bonded[0].position], np.zeros(3), dimensions
+            )
+        # case3
+        if len(heavy_bonded) == 1 and len(light_bonded) >= 1:
+            custom_axes = self.get_custom_axes(
+                atom.position,
+                [heavy_bonded[0].position],
+                light_bonded[0].position,
+                dimensions,
+            )
+        # case4, not used in Jon's 2019 paper code, use case5 instead
+        # case5
+        if len(heavy_bonded) >= 2:
+            custom_axes = self.get_custom_axes(
+                atom.position,
+                heavy_bonded.positions,
+                heavy_bonded[1].position,
+                dimensions,
+            )
+
+        if custom_moment_of_inertia is None:
+            # find moment of inertia using custom axes and atom position as COM
+            custom_moment_of_inertia = self.get_custom_moment_of_inertia(
+                UA, custom_axes, atom.position, dimensions
+            )
+
+        # get the moment of inertia from the custom axes
+        if custom_axes is not None:
+            # flip axes to face correct way wrt COM
+            custom_axes = self.get_flipped_axes(
+                UA, custom_axes, atom.position, dimensions
+            )
+
+        return custom_axes, custom_moment_of_inertia
+
+    def find_bonded_atoms(self, atom_idx: int, system):
+        """
+        for a given atom, find its bonded heavy and H atoms
+
+        Args:
+            atom_idx: atom index to find bonded heavy atom for
+            system: mdanalysis instance of all atoms in current frame
+
+        Returns:
+            bonded_heavy_atoms: MDAnalysis instance of bonded heavy atoms
+            bonded_H_atoms: MDAnalysis instance of bonded hydrogen atoms
+        """
+        bonded_atoms = system.select_atoms(f"bonded index {atom_idx}")
+        bonded_heavy_atoms = bonded_atoms.select_atoms("mass 2 to 999")
+        bonded_H_atoms = bonded_atoms.select_atoms("mass 1 to 1.1")
+        return bonded_heavy_atoms, bonded_H_atoms
+
+    def get_vanilla_axes(self, molecule):
+        """
+        Compute the principal axes and sorted moments of inertia for a molecule.
+
+        This method computes the translationally invariant principal axes and
+        corresponding moments of inertia for a molecular selection using the
+        default MDAnalysis routines. The molecule is first made whole to ensure
+        correct handling of periodic boundary conditions.
+
+        The moments of inertia are obtained by diagonalising the moment of inertia
+        tensor and are returned sorted from largest to smallest magnitude.
+
+        Args:
+            molecule (MDAnalysis.core.groups.AtomGroup):
+                AtomGroup representing the molecule or bead for which the axes
+                and moments of inertia are to be computed.
+
+        Returns:
+            Tuple[np.ndarray, np.ndarray]:
+                A tuple containing:
+
+                - principal_axes (np.ndarray):
+                    Array of shape ``(3, 3)`` whose rows correspond to the
+                    principal axes of the molecule.
+                - moment_of_inertia (np.ndarray):
+                    Array of shape ``(3,)`` containing the moments of inertia
+                    sorted in descending order.
+        """
+        moment_of_inertia = molecule.moment_of_inertia(unwrap=True)
+        make_whole(molecule.atoms)
+        principal_axes = molecule.principal_axes()
+
+        eigenvalues, _eigenvectors = np.linalg.eig(moment_of_inertia)
+
+        # Sort eigenvalues from largest to smallest magnitude
+        order = np.argsort(np.abs(eigenvalues))[::-1]
+        moment_of_inertia = eigenvalues[order]
+
+        return principal_axes, moment_of_inertia
+
+    def get_custom_axes(
+        self, a: np.ndarray, b_list: list, c: np.ndarray, dimensions: np.ndarray
+    ):
+        r"""
+        For atoms a, b_list and c, calculate the axis to rotate forces around:
+
+        - axis1: use the normalised vector ab as axis1. If there is more than one bonded
+        heavy atom (HA), average over all the normalised vectors calculated from b_list
+        and use this as axis1). b_list contains all the bonded heavy atom
+        coordinates.
+
+        - axis2: use the cross product of normalised vector ac and axis1 as axis2.
+        If there are more than two bonded heavy atoms, then use normalised vector
+        b[0]c to cross product with axis1, this gives the axis perpendicular
+        (represented by |_ symbol below) to axis1.
+
+        - axis3: the cross product of axis1 and axis2, which is perpendicular to
+        axis1 and axis2.
+
+        Args:
+            a: central united-atom coordinates (3,)
+            b_list: list of heavy bonded atom positions (3,N)
+            c: atom coordinates of either a second heavy atom or a hydrogen atom
+            if there are no other bonded heavy atoms in b_list (where N=1 in b_list)
+            (3,)
+            dimensions: dimensions of the simulation box (3,)
+
+        ::
+
+            a          1 = norm_ab
+           / \         2 = |_ norm_ab and norm_ac (use bc if more than 2 HAs)
+          /   \        3 = |_ 1 and 2
+        b       c
+
+        Returns:
+            custom_axes: (3,3) array of the axes used to rotate forces
+        """
+        unscaled_axis1 = np.zeros(3)
+        # average of all heavy atom covalent bond vectors for axis1
+        for b in b_list:
+            ab_vector = self.get_vector(a, b, dimensions)
+            unscaled_axis1 += ab_vector
+        if len(b_list) >= 2:
+            # use the first heavy bonded atom as atom a
+            ac_vector = self.get_vector(c, b_list[0], dimensions)
+        else:
+            ac_vector = self.get_vector(c, a, dimensions)
+
+        unscaled_axis2 = np.cross(ac_vector, unscaled_axis1)
+        unscaled_axis3 = np.cross(unscaled_axis2, unscaled_axis1)
+
+        unscaled_custom_axes = np.array(
+            (unscaled_axis1, unscaled_axis2, unscaled_axis3)
+        )
+        mod = np.sqrt(np.sum(unscaled_custom_axes**2, axis=1))
+        scaled_custom_axes = unscaled_custom_axes / mod[:, np.newaxis]
+
+        return scaled_custom_axes
+
+    def get_custom_moment_of_inertia(
+        self,
+        UA,
+        custom_rotation_axes: np.ndarray,
+        center_of_mass: np.ndarray,
+        dimensions: np.ndarray,
+    ):
+        """
+        Get the moment of inertia (specifically used for the united atom level)
+        from a set of rotation axes and a given center of mass
+        (COM is usually the heavy atom position in a UA).
+
+        Args:
+            UA: MDAnalysis instance of a united-atom
+            custom_rotation_axes: (3,3) arrray of rotation axes
+            center_of_mass: (3,) center of mass for collection of atoms N
+
+        Returns:
+            custom_moment_of_inertia: (3,) array for moment of inertia
+        """
+        translated_coords = self.get_vector(center_of_mass, UA.positions, dimensions)
+        custom_moment_of_inertia = np.zeros(3)
+        for coord, mass in zip(translated_coords, UA.masses):
+            axis_component = np.sum(
+                np.cross(custom_rotation_axes, coord) ** 2 * mass, axis=1
+            )
+            custom_moment_of_inertia += axis_component
+
+        # Remove lowest MOI degree of freedom if UA only has a single bonded H
+        if len(UA) == 2:
+            order = custom_moment_of_inertia.argsort()[::-1]  # decending order
+            custom_moment_of_inertia[order[-1]] = 0
+
+        return custom_moment_of_inertia
+
+    def get_flipped_axes(self, UA, custom_axes, center_of_mass, dimensions):
+        """
+        For a given set of custom axes, ensure the axes are pointing in the
+        correct direction wrt the heavy atom position and the chosen center
+        of mass.
+
+        Args:
+            UA: MDAnalysis instance of a united-atom
+            custom_axes: (3,3) array of the rotation axes
+            center_of_mass: (3,) array for center of mass (usually HA position)
+            dimensions: (3,) array of system box dimensions.
+        """
+        # sorting out PIaxes for MoI for UA fragment
+
+        # get dot product of Paxis1 and CoM->atom1 vect
+        # will just be [0,0,0]
+        RRaxis = self.get_vector(UA[0].position, center_of_mass, dimensions)
+
+        # flip each Paxis if its pointing out of UA
+        custom_axis = np.sum(custom_axes**2, axis=1)
+        custom_axes_flipped = custom_axes / custom_axis**0.5
+        for i in range(3):
+            dotProd1 = np.dot(custom_axes_flipped[i], RRaxis)
+            custom_axes_flipped[i] = np.where(
+                dotProd1 < 0, -custom_axes_flipped[i], custom_axes_flipped[i]
+            )
+        return custom_axes_flipped
+
+    def get_vector(self, a: np.ndarray, b: np.ndarray, dimensions: np.ndarray):
+        """
+        For vector of two coordinates over periodic boundary conditions (PBCs).
+
+        Args:
+            a: (N,3) array of atom cooordinates
+            b: (3,) array of atom cooordinates
+            dimensions: (3,) array of system box dimensions.
+
+        Returns:
+            delta_wrapped: (N,3) array of the vector
+        """
+        delta = b - a
+        delta -= dimensions * np.round(delta / dimensions)
+
+        return delta
+
+    def get_moment_of_inertia_tensor(
+        self,
+        center_of_mass: np.ndarray,
+        positions: np.ndarray,
+        masses: list,
+        dimensions: np.array,
+    ) -> np.ndarray:
+        """
+        Calculate a custom moment of inertia tensor.
+        E.g., for cases where the mass list will contain masses of UAs rather than
+        individual atoms and the postions will be those for the UAs only
+        (excluding the H atoms coordinates).
+
+        Args:
+            center_of_mass: a (3,) array of the chosen center of mass
+            positions: a (N,3) array of point positions
+            masses: a (N,) list of point masses
+
+        Returns:
+            moment_of_inertia_tensor: a (3,3) moment of inertia tensor
+        """
+        r = self.get_vector(center_of_mass, positions, dimensions)
+        r2 = np.sum(r**2, axis=1)
+        moment_of_inertia_tensor = np.eye(3) * np.sum(masses * r2)
+        moment_of_inertia_tensor -= np.einsum("i,ij,ik->jk", masses, r, r)
+
+        return moment_of_inertia_tensor
+
+    def get_custom_principal_axes(
+        self, moment_of_inertia_tensor: np.ndarray
+    ) -> tuple[np.ndarray, np.ndarray]:
+        """
+        Principal axes and centre of axes from the ordered eigenvalues
+        and eigenvectors of a moment of inertia tensor. This function allows for
+        a custom moment of inertia tensor to be used, which isn't possible with
+        the built-in MDAnalysis principal_axes() function.
+
+        Args:
+            moment_of_inertia_tensor: a (3,3) array of a custom moment of
+             inertia tensor
+
+        Returns:
+            principal_axes: a (3,3) array for the principal axes
+            moment_of_inertia: a (3,) array of the principal axes center
+        """
+        eigenvalues, eigenvectors = np.linalg.eig(moment_of_inertia_tensor)
+        order = abs(eigenvalues).argsort()[::-1]  # decending order
+        transposed = np.transpose(eigenvectors)  # turn columns to rows
+        moment_of_inertia = eigenvalues[order]
+        principal_axes = transposed[order]
+
+        # point z axis in correct direction, as per Jon's code
+        cross_xy = np.cross(principal_axes[0], principal_axes[1])
+        dot_z = np.dot(cross_xy, principal_axes[2])
+        if dot_z < 0:
+            principal_axes[2] *= -1
+
+        return principal_axes, moment_of_inertia
+
+    def get_UA_masses(self, molecule) -> list[float]:
+        """
+        For a given molecule, return a list of masses of UAs
+        (combination of the heavy atoms + bonded hydrogen atoms. This list is used to
+        get the moment of inertia tensor for molecules larger than one UA.
+
+        Args:
+            molecule: mdanalysis instance of molecule
+
+        Returns:
+            UA_masses: list of masses for each UA in a molecule
+        """
+        UA_masses = []
+        for atom in molecule:
+            if atom.mass > 1.1:
+                UA_mass = atom.mass
+                bonded_atoms = molecule.select_atoms(f"bonded index {atom.index}")
+                bonded_H_atoms = bonded_atoms.select_atoms("mass 1 to 1.1")
+                for H in bonded_H_atoms:
+                    UA_mass += H.mass
+                UA_masses.append(UA_mass)
+            else:
+                continue
+        return UA_masses
diff --git a/CodeEntropy/config/arg_config_manager.py b/CodeEntropy/config/arg_config_manager.py
index e733952..813c6a3 100644
--- a/CodeEntropy/config/arg_config_manager.py
+++ b/CodeEntropy/config/arg_config_manager.py
@@ -85,6 +85,18 @@
         "help": "How to group molecules for averaging",
         "default": "molecules",
     },
+    "combined_forcetorque": {
+        "type": bool,
+        "help": """Use combined force-torque matrix for residue
+            level vibrational entropies""",
+        "default": True,
+    },
+    "customised_axes": {
+        "type": bool,
+        "help": """Use bonded axes to rotate forces for UA
+            level vibrational entropies""",
+        "default": True,
+    },
 }
 
 
diff --git a/CodeEntropy/entropy.py b/CodeEntropy/entropy.py
index c5e72c1..9b7016e 100644
--- a/CodeEntropy/entropy.py
+++ b/CodeEntropy/entropy.py
@@ -122,7 +122,7 @@ def execute(self):
         else:
             nonwater_groups.update(water_groups)
 
-        force_matrices, torque_matrices, frame_counts = (
+        force_matrices, torque_matrices, forcetorque_matrices, frame_counts = (
             self._level_manager.build_covariance_matrices(
                 self,
                 reduced_atom,
@@ -133,6 +133,8 @@ def execute(self):
                 step,
                 number_frames,
                 self._args.force_partitioning,
+                self._args.combined_forcetorque,
+                self._args.customised_axes,
             )
         )
 
@@ -155,6 +157,7 @@ def execute(self):
             nonwater_groups,
             force_matrices,
             torque_matrices,
+            forcetorque_matrices,
             states_ua,
             states_res,
             frame_counts,
@@ -230,6 +233,7 @@ def _compute_entropies(
         groups,
         force_matrices,
         torque_matrices,
+        forcetorque_matrices,
         states_ua,
         states_res,
         frame_counts,
@@ -309,6 +313,7 @@ def _compute_entropies(
                         f"Level: {level}",
                     )
                     highest = level == levels[groups[group_id][0]][-1]
+                    forcetorque_matrix = None
 
                     if level == "united_atom":
                         self._process_united_atom_entropy(
@@ -326,6 +331,8 @@ def _compute_entropies(
                         )
 
                     elif level == "residue":
+                        if highest:
+                            forcetorque_matrix = forcetorque_matrices["res"][group_id]
                         self._process_vibrational_entropy(
                             group_id,
                             mol,
@@ -334,6 +341,7 @@ def _compute_entropies(
                             level,
                             force_matrices["res"][group_id],
                             torque_matrices["res"][group_id],
+                            forcetorque_matrix,
                             highest,
                         )
 
@@ -347,6 +355,8 @@ def _compute_entropies(
                         )
 
                     elif level == "polymer":
+                        if highest:
+                            forcetorque_matrix = forcetorque_matrices["poly"][group_id]
                         self._process_vibrational_entropy(
                             group_id,
                             mol,
@@ -355,6 +365,7 @@ def _compute_entropies(
                             level,
                             force_matrices["poly"][group_id],
                             torque_matrices["poly"][group_id],
+                            forcetorque_matrix,
                             highest,
                         )
 
@@ -530,6 +541,7 @@ def _process_vibrational_entropy(
         level,
         force_matrix,
         torque_matrix,
+        forcetorque_matrix,
         highest,
     ):
         """
@@ -548,21 +560,43 @@ def _process_vibrational_entropy(
         # Find the relevant force and torque matrices and tidy them up
         # by removing rows and columns that are all zeros
 
-        force_matrix = self._level_manager.filter_zero_rows_columns(force_matrix)
+        if forcetorque_matrix is not None:
+            forcetorque_matrix = self._level_manager.filter_zero_rows_columns(
+                forcetorque_matrix
+            )
 
-        torque_matrix = self._level_manager.filter_zero_rows_columns(torque_matrix)
+            S_FTtrans = ve.vibrational_entropy_calculation(
+                forcetorque_matrix, "forcetorqueTRANS", self._args.temperature, highest
+            )
+            S_FTrot = ve.vibrational_entropy_calculation(
+                forcetorque_matrix, "forcetorqueROT", self._args.temperature, highest
+            )
 
-        # Calculate the vibrational entropy
-        S_trans = ve.vibrational_entropy_calculation(
-            force_matrix, "force", self._args.temperature, highest
-        )
-        S_rot = ve.vibrational_entropy_calculation(
-            torque_matrix, "torque", self._args.temperature, highest
-        )
+            self._data_logger.add_results_data(
+                group_id, level, "FTmat-Transvibrational", S_FTtrans
+            )
+            self._data_logger.add_results_data(
+                group_id, level, "FTmat-Rovibrational", S_FTrot
+            )
 
-        # Print the vibrational entropy for the molecule group
-        self._data_logger.add_results_data(group_id, level, "Transvibrational", S_trans)
-        self._data_logger.add_results_data(group_id, level, "Rovibrational", S_rot)
+        else:
+            force_matrix = self._level_manager.filter_zero_rows_columns(force_matrix)
+
+            torque_matrix = self._level_manager.filter_zero_rows_columns(torque_matrix)
+
+            # Calculate the vibrational entropy
+            S_trans = ve.vibrational_entropy_calculation(
+                force_matrix, "force", self._args.temperature, highest
+            )
+            S_rot = ve.vibrational_entropy_calculation(
+                torque_matrix, "torque", self._args.temperature, highest
+            )
+
+            # Print the vibrational entropy for the molecule group
+            self._data_logger.add_results_data(
+                group_id, level, "Transvibrational", S_trans
+            )
+            self._data_logger.add_results_data(group_id, level, "Rovibrational", S_rot)
 
         residue_group = "_".join(
             sorted(set(res.resname for res in mol_container.residues))
@@ -899,6 +933,7 @@ def vibrational_entropy_calculation(self, matrix, matrix_type, temp, highest_lev
         """
         # N beads at a level => 3N x 3N covariance matrix => 3N eigenvalues
         # Get eigenvalues of the given matrix and change units to SI units
+        logger.debug(f"matrix_type: {matrix_type}")
         lambdas = la.eigvals(matrix)
         logger.debug(f"Eigenvalues (lambdas) before unit change: {lambdas}")
 
@@ -939,6 +974,11 @@ def vibrational_entropy_calculation(self, matrix, matrix_type, temp, highest_lev
             else:
                 S_vib_total = sum(S_components[6:])
 
+        elif matrix_type == "forcetorqueTRANS":  # three lowest are translations
+            S_vib_total = sum(S_components[:3])
+        elif matrix_type == "forcetorqueROT":  # three highest are rotations
+            S_vib_total = sum(S_components[3:])
+
         else:  # torque covariance matrix - we always take all values into account
             S_vib_total = sum(S_components)
 
diff --git a/CodeEntropy/levels.py b/CodeEntropy/levels.py
index 9d43670..1426245 100644
--- a/CodeEntropy/levels.py
+++ b/CodeEntropy/levels.py
@@ -1,6 +1,7 @@
 import logging
 
 import numpy as np
+from MDAnalysis.lib.mdamath import make_whole
 from rich.progress import (
     BarColumn,
     Progress,
@@ -9,6 +10,8 @@
     TimeElapsedColumn,
 )
 
+from CodeEntropy.axes import AxesManager
+
 logger = logging.getLogger(__name__)
 
 
@@ -84,6 +87,7 @@ def get_matrices(
         force_matrix,
         torque_matrix,
         force_partitioning,
+        customised_axes,
     ):
         """
         Compute and accumulate force/torque covariance matrices for a given level.
@@ -96,6 +100,8 @@ def get_matrices(
           to.
           force_partitioning (float): Factor to adjust force contributions,
           default is 0.5.
+          customised_axes (bool): Whether to use customised axes for rotating
+          forces.
 
         Returns:
           force_matrix (np.ndarray): Accumulated force covariance matrix.
@@ -115,11 +121,27 @@ def get_matrices(
         # Calculate forces/torques for each bead
         for bead_index in range(number_beads):
             bead = list_of_beads[bead_index]
+
             # Set up axes
             # translation and rotation use different axes
             # how the axes are defined depends on the level
-            trans_axes = data_container.atoms.principal_axes()
-            rot_axes = np.real(bead.principal_axes())
+            axes_manager = AxesManager()
+            if level == "united_atom" and customised_axes:
+                trans_axes, rot_axes, center, moment_of_inertia = (
+                    axes_manager.get_UA_axes(data_container, bead_index)
+                )
+            elif level == "residue" and customised_axes:
+                trans_axes, rot_axes, center, moment_of_inertia = (
+                    axes_manager.get_residue_axes(data_container, bead_index)
+                )
+            else:
+                make_whole(data_container.atoms)
+                make_whole(bead)
+                trans_axes = data_container.atoms.principal_axes()
+                rot_axes = np.real(bead.principal_axes())
+                eigenvalues, _ = np.linalg.eig(bead.moment_of_inertia(unwrap=True))
+                moment_of_inertia = sorted(eigenvalues, reverse=True)
+                center = bead.center_of_mass(unwrap=True)
 
             # Sort out coordinates, forces, and torques for each atom in the bead
             weighted_forces[bead_index] = self.get_weighted_forces(
@@ -130,7 +152,12 @@ def get_matrices(
                 force_partitioning,
             )
             weighted_torques[bead_index] = self.get_weighted_torques(
-                data_container, bead, rot_axes, force_partitioning
+                bead,
+                rot_axes,
+                center,
+                force_partitioning,
+                moment_of_inertia,
+                axes_manager,
             )
 
         # Create covariance submatrices
@@ -167,6 +194,7 @@ def get_matrices(
         # Enforce consistent shape before accumulation
         if force_matrix is None:
             force_matrix = np.zeros_like(force_block)
+            force_matrix = force_block  # add first set of forces
         elif force_matrix.shape != force_block.shape:
             raise ValueError(
                 f"Inconsistent force matrix shape: existing "
@@ -177,6 +205,7 @@ def get_matrices(
 
         if torque_matrix is None:
             torque_matrix = np.zeros_like(torque_block)
+            torque_matrix = torque_block  # add first set of torques
         elif torque_matrix.shape != torque_block.shape:
             raise ValueError(
                 f"Inconsistent torque matrix shape: existing "
@@ -187,6 +216,118 @@ def get_matrices(
 
         return force_matrix, torque_matrix
 
+    def get_combined_forcetorque_matrices(
+        self,
+        data_container,
+        level,
+        highest_level,
+        forcetorque_matrix,
+        force_partitioning,
+        customised_axes,
+    ):
+        """
+        Compute and accumulate combined force/torque covariance matrices for
+        a given level.
+
+        Parameters:
+          data_container (MDAnalysis.Universe): Data for a molecule or residue.
+          level (str): 'polymer', 'residue', or 'united_atom'.
+          highest_level (bool): Whether this is the top (largest bead size) level.
+          forcetorque_matrix (np.ndarray or None): Accumulated matrices to add
+          to.
+          force_partitioning (float): Factor to adjust force contributions,
+          default is 0.5.
+          customised_axes (bool): Whether to use customised axes for rotating
+          forces.
+
+        Returns:
+          forcetorque_matrix (np.ndarray): Accumulated torque covariance matrix.
+        """
+
+        # Make beads
+        list_of_beads = self.get_beads(data_container, level)
+
+        # number of beads and frames in trajectory
+        number_beads = len(list_of_beads)
+
+        # initialize force and torque arrays
+        weighted_forces = [None for _ in range(number_beads)]
+        weighted_torques = [None for _ in range(number_beads)]
+
+        # Create axes manager for custom axes
+        axes_manager = AxesManager()
+
+        # Calculate forces/torques for each bead
+        for bead_index in range(number_beads):
+            bead = list_of_beads[bead_index]
+            # Set up axes
+            # translation and rotation use different axes
+            # how the axes are defined depends on the level
+            if level == "residue" and customised_axes:
+                trans_axes, rot_axes, center, moment_of_inertia = (
+                    axes_manager.get_residue_axes(data_container, bead_index)
+                )
+            else:
+                # ensure molecule is whole for PA calcs
+                make_whole(data_container.atoms)
+                make_whole(bead)
+                trans_axes = data_container.atoms.principal_axes()
+                rot_axes, moment_of_inertia = axes_manager.get_vanilla_axes(bead)
+                center = bead.center_of_mass(unwrap=True)
+
+            # Sort out coordinates, forces, and torques for each atom in the bead
+            weighted_forces[bead_index] = self.get_weighted_forces(
+                data_container,
+                bead,
+                trans_axes,
+                highest_level,
+                force_partitioning,
+            )
+            weighted_torques[bead_index] = self.get_weighted_torques(
+                bead,
+                rot_axes,
+                center,
+                force_partitioning,
+                moment_of_inertia,
+                axes_manager,
+            )
+
+        # Create covariance submatrices
+        forcetorque_submatrix = [
+            [0 for _ in range(number_beads)] for _ in range(number_beads)
+        ]
+
+        for i in range(number_beads):
+            for j in range(i, number_beads):
+                ft_sub = self.create_FTsubmatrix(
+                    np.concatenate((weighted_forces[i], weighted_torques[i])),
+                    np.concatenate((weighted_forces[j], weighted_torques[j])),
+                )
+                forcetorque_submatrix[i][j] = ft_sub
+                forcetorque_submatrix[j][i] = ft_sub.T
+
+        # Convert block matrices to full matrix
+        forcetorque_block = np.block(
+            [
+                [forcetorque_submatrix[i][j] for j in range(number_beads)]
+                for i in range(number_beads)
+            ]
+        )
+
+        # Enforce consistent shape before accumulation
+        if forcetorque_matrix is None:
+            forcetorque_matrix = np.zeros_like(forcetorque_block)
+            forcetorque_matrix = forcetorque_block  # add first set of torques
+        elif forcetorque_matrix.shape != forcetorque_block.shape:
+            raise ValueError(
+                f"Inconsistent forcetorque matrix shape: existing "
+                f"{forcetorque_matrix.shape}, new {forcetorque_block.shape}"
+            )
+        else:
+            forcetorque_matrix = forcetorque_block
+
+        return forcetorque_matrix
+
     def get_beads(self, data_container, level):
         """
         Function to define beads depending on the level in the hierarchy.
@@ -296,7 +437,15 @@ def get_weighted_forces(
 
         return weighted_force
 
-    def get_weighted_torques(self, data_container, bead, rot_axes, force_partitioning):
+    def get_weighted_torques(
+        self,
+        bead,
+        rot_axes,
+        center,
+        force_partitioning,
+        moment_of_inertia,
+        axes_manager,
+    ):
         """
         Compute moment-of-inertia weighted torques for a bead.
 
@@ -311,9 +460,7 @@ def get_weighted_torques(self, data_container, bead, rot_axes, force_partitionin
         the sorted eigenvalues of the moment of inertia tensor.
 
         To ensure numerical stability:
-        - Torque components that are effectively zero are skipped.
-        - Zero moments of inertia result in zero weighted torque with a warning.
-        - Negative moments of inertia raise an error.
+        - Torque components that are effectively zero, zero or negative are skipped.
 
         Parameters
         ----------
@@ -326,55 +473,51 @@ def get_weighted_torques(self, data_container, bead, rot_axes, force_partitionin
         force_partitioning : float
             Scaling factor applied to forces to avoid over-counting correlated
             contributions (typically 0.5).
+        moment_of_inertia : np.ndarray
+            Moment of inertia (3,)
+        customised_axes: bool
+            Whether to use customised axes for rotating forces.
 
         Returns
         -------
         weighted_torque : np.ndarray
             Moment-of-inertia weighted torque acting on the bead.
-
-        Raises
-        ------
-        ValueError
-            If a negative principal moment of inertia is encountered.
         """
-        torques = np.zeros((3,))
-        weighted_torque = np.zeros((3,))
-        moment_of_inertia = np.zeros(3)
 
-        for atom in bead.atoms:
-            coords_rot = (
-                data_container.atoms[atom.index].position - bead.center_of_mass()
-            )
-            coords_rot = np.matmul(rot_axes, coords_rot)
-            forces_rot = np.matmul(rot_axes, data_container.atoms[atom.index].force)
-
-            forces_rot = force_partitioning * forces_rot
-
-            torques_local = np.cross(coords_rot, forces_rot)
-            torques += torques_local
-
-        eigenvalues, _ = np.linalg.eig(bead.moment_of_inertia())
-        moments_of_inertia = sorted(eigenvalues, reverse=True)
+        # translate and rotate positions and forces
+        # translated_coords = bead.positions - center
+        translated_coords = axes_manager.get_vector(
+            center, bead.positions, bead.dimensions[:3]
+        )
+        rotated_coords = np.tensordot(translated_coords, rot_axes.T, axes=1)
+        rotated_forces = np.tensordot(bead.forces, rot_axes.T, axes=1)
+        # scale forces
+        rotated_forces *= force_partitioning
+        # get torques here
+        torques = np.cross(rotated_coords, rotated_forces)
+        torques = np.sum(torques, axis=0)
 
+        weighted_torque = np.zeros((3,))
         for dimension in range(3):
             if np.isclose(torques[dimension], 0):
                 weighted_torque[dimension] = 0
                 continue
 
-            if np.isclose(moments_of_inertia[dimension], 0):
+            if np.iscomplex(moment_of_inertia[dimension]):
                 weighted_torque[dimension] = 0
-                logger.warning("Zero moment of inertia. Setting torque to 0")
                 continue
 
-            if moments_of_inertia[dimension] < 0:
-                raise ValueError(
-                    f"Negative value encountered for moment of inertia: "
-                    f"{moment_of_inertia[dimension]} "
-                    f"Cannot compute weighted torque."
-                )
+            if moment_of_inertia[dimension] == 0:
+                moment_of_inertia[dimension] = 0
+                continue
 
+            if moment_of_inertia[dimension] < 0:
+                moment_of_inertia[dimension] = 0
+                continue
+
+            # Compute weighted torque
             weighted_torque[dimension] = torques[dimension] / np.sqrt(
-                moments_of_inertia[dimension]
+                moment_of_inertia[dimension]
             )
 
         logger.debug(f"Weighted Torque: {weighted_torque}")
@@ -407,6 +550,30 @@ def create_submatrix(self, data_i, data_j):
 
         return submatrix
 
+    def create_FTsubmatrix(self, data_i, data_j):
+        """
+        Function for making covariance matrices.
+
+        Args
+        -----
+        data_i : values for bead i
+        data_j : values for bead j
+
+        Returns
+        ------
+        submatrix : 6x6 matrix for the covariance between i and j
+        """
+
+        # Start with 6 by 6 matrix of zeros
+        submatrix = np.zeros((6, 6))
+
+        # For each frame calculate the outer product (cross product) of the data from
+        # the two beads and add the result to the submatrix
+        outer_product_matrix = np.outer(data_i, data_j)
+        submatrix = np.add(submatrix, outer_product_matrix)
+
+        return submatrix
+
     def build_covariance_matrices(
         self,
         entropy_manager,
@@ -418,6 +585,8 @@ def build_covariance_matrices(
         step,
         number_frames,
         force_partitioning,
+        combined_forcetorque,
+        customised_axes,
     ):
         """
         Construct average force and torque covariance matrices for all molecules and
@@ -443,7 +612,8 @@ def build_covariance_matrices(
             Total number of frames to process.
         force_partitioning : float
             Factor to adjust force contributions, default is 0.5.
-
+        combined_forcetorque : bool
+         Whether to use combined forcetorque covariance matrix.
 
         Returns
         -------
@@ -466,6 +636,12 @@ def build_covariance_matrices(
             "poly": [None] * number_groups,
         }
 
+        forcetorque_avg = {
+            "ua": {},
+            "res": [None] * number_groups,
+            "poly": [None] * number_groups,
+        }
+
         total_steps = len(reduced_atom.trajectory[start:end:step])
         total_items = (
             sum(len(levels[mol_id]) for mols in groups.values() for mol_id in mols)
@@ -524,13 +700,16 @@ def build_covariance_matrices(
                                 number_frames,
                                 force_avg,
                                 torque_avg,
+                                forcetorque_avg,
                                 frame_counts,
                                 force_partitioning,
+                                combined_forcetorque,
+                                customised_axes,
                             )
 
                             progress.advance(task)
 
-        return force_avg, torque_avg, frame_counts
+        return force_avg, torque_avg, forcetorque_avg, frame_counts
 
     def update_force_torque_matrices(
         self,
@@ -543,8 +722,11 @@ def update_force_torque_matrices(
         num_frames,
         force_avg,
         torque_avg,
+        forcetorque_avg,
         frame_counts,
         force_partitioning,
+        combined_forcetorque,
+        customised_axes,
     ):
         """
         Update the running averages of force and torque covariance matrices
@@ -586,7 +768,11 @@ def update_force_torque_matrices(
             Dictionary holding the count of frames processed for each molecule/level
             combination.
         force_partitioning : float
-         Factor to adjust force contributions, default is 0.5.
+            Factor to adjust force contributions, default is 0.5.
+        combined_forcetorque : bool
+            Whether to use combined forcetorque covariance matrix.
+        customised_axes: bool
+            Whether to use bonded axes for UA rovib calculations
         Returns
         -------
         None
@@ -619,6 +805,7 @@ def update_force_torque_matrices(
                     None if key not in force_avg["ua"] else force_avg["ua"][key],
                     None if key not in torque_avg["ua"] else torque_avg["ua"][key],
                     force_partitioning,
+                    customised_axes,
                 )
 
                 if key not in force_avg["ua"]:
@@ -641,28 +828,58 @@ def update_force_torque_matrices(
             # Build the matrices, adding data from each timestep
             # Being careful for the first timestep when data has not yet
             # been added to the matrices
-            f_mat, t_mat = self.get_matrices(
-                mol,
-                level,
-                highest,
-                None if force_avg[key][group_id] is None else force_avg[key][group_id],
-                (
-                    None
-                    if torque_avg[key][group_id] is None
-                    else torque_avg[key][group_id]
-                ),
-                force_partitioning,
-            )
+            if highest and combined_forcetorque:
+                # use combined forcetorque covariance matrix for the highest level only
+                ft_mat = self.get_combined_forcetorque_matrices(
+                    mol,
+                    level,
+                    highest,
+                    (
+                        None
+                        if forcetorque_avg[key][group_id] is None
+                        else forcetorque_avg[key][group_id]
+                    ),
+                    force_partitioning,
+                    customised_axes,
+                )
 
-            if force_avg[key][group_id] is None:
-                force_avg[key][group_id] = f_mat.copy()
-                torque_avg[key][group_id] = t_mat.copy()
-                frame_counts[key][group_id] = 1
+                if forcetorque_avg[key][group_id] is None:
+                    forcetorque_avg[key][group_id] = ft_mat.copy()
+                    frame_counts[key][group_id] = 1
+                else:
+                    frame_counts[key][group_id] += 1
+                    n = frame_counts[key][group_id]
+                    forcetorque_avg[key][group_id] += (
+                        ft_mat - forcetorque_avg[key][group_id]
+                    ) / n
             else:
-                frame_counts[key][group_id] += 1
-                n = frame_counts[key][group_id]
-                force_avg[key][group_id] += (f_mat - force_avg[key][group_id]) / n
-                torque_avg[key][group_id] += (t_mat - torque_avg[key][group_id]) / n
+                f_mat, t_mat = self.get_matrices(
+                    mol,
+                    level,
+                    highest,
+                    (
+                        None
+                        if force_avg[key][group_id] is None
+                        else force_avg[key][group_id]
+                    ),
+                    (
+                        None
+                        if torque_avg[key][group_id] is None
+                        else torque_avg[key][group_id]
+                    ),
+                    force_partitioning,
+                    customised_axes,
+                )
+
+                if force_avg[key][group_id] is None:
+                    force_avg[key][group_id] = f_mat.copy()
+                    torque_avg[key][group_id] = t_mat.copy()
+                    frame_counts[key][group_id] = 1
+                else:
+                    frame_counts[key][group_id] += 1
+                    n = frame_counts[key][group_id]
+                    force_avg[key][group_id] += (f_mat - force_avg[key][group_id]) / n
+                    torque_avg[key][group_id] += (t_mat - torque_avg[key][group_id]) / n
 
         return frame_counts
 
diff --git a/CodeEntropy/mda_universe_operations.py b/CodeEntropy/mda_universe_operations.py
index d7f739e..306dfd3 100644
--- a/CodeEntropy/mda_universe_operations.py
+++ b/CodeEntropy/mda_universe_operations.py
@@ -54,8 +54,15 @@ def new_U_select_frame(self, u, start=None, end=None, step=1):
             .run()
             .results["timeseries"][start:end:step]
         )
+        dimensions = (
+            AnalysisFromFunction(lambda ag: ag.dimensions.copy(), select_atom)
+            .run()
+            .results["timeseries"][start:end:step]
+        )
         u2 = mda.Merge(select_atom)
-        u2.load_new(coordinates, format=MemoryReader, forces=forces)
+        u2.load_new(
+            coordinates, format=MemoryReader, forces=forces, dimensions=dimensions
+        )
         logger.debug(f"MDAnalysis.Universe - reduced universe: {u2}")
 
         return u2
@@ -89,8 +96,15 @@ def new_U_select_atom(self, u, select_string="all"):
             .run()
             .results["timeseries"]
         )
+        dimensions = (
+            AnalysisFromFunction(lambda ag: ag.dimensions.copy(), select_atom)
+            .run()
+            .results["timeseries"]
+        )
         u2 = mda.Merge(select_atom)
-        u2.load_new(coordinates, format=MemoryReader, forces=forces)
+        u2.load_new(
+            coordinates, format=MemoryReader, forces=forces, dimensions=dimensions
+        )
         logger.debug(f"MDAnalysis.Universe - reduced universe: {u2}")
 
         return u2
@@ -148,12 +162,18 @@ def merge_forces(self, tprfile, trrfile, forcefile, fileformat=None, kcal=False)
             .results["timeseries"]
         )
 
+        dimensions = (
+            AnalysisFromFunction(lambda ag: ag.dimensions.copy(), select_atom)
+            .run()
+            .results["timeseries"]
+        )
+
         if kcal:
             # Convert from kcal to kJ
             forces *= 4.184
 
         logger.debug("Merging forces with coordinates universe.")
         new_universe = mda.Merge(select_atom)
-        new_universe.load_new(coordinates, forces=forces)
+        new_universe.load_new(coordinates, forces=forces, dimensions=dimensions)
 
         return new_universe
diff --git a/docs/getting_started.rst b/docs/getting_started.rst
index 3d103e4..c7f8c4d 100644
--- a/docs/getting_started.rst
+++ b/docs/getting_started.rst
@@ -137,6 +137,18 @@ The top_traj_file argument is necessary to identify your simulation data, the ot
      - How to group molecules for averaging
      - ``molecules``
      - ``str``
+   * - ``--kcal_force_units``
+     - Set input units as kcal/mol 
+     - ``bool``
+     - ``False``
+   * - ``--combined_forcetorque``
+     - Use the combined force-torque covariance matrix for the highest level to match the 2019 paper
+     - ``bool``
+     - ``True``
+   * - ``--customised_axes``
+     - Use custom bonded axes to get COM, MOI and PA that match the 2019 paper
+     - ``bool``
+     - ``True``
 
 Averaging
 ^^^^^^^^^
diff --git a/tests/test_CodeEntropy/test_axes.py b/tests/test_CodeEntropy/test_axes.py
new file mode 100644
index 0000000..c5e2d38
--- /dev/null
+++ b/tests/test_CodeEntropy/test_axes.py
@@ -0,0 +1,479 @@
+from unittest.mock import MagicMock, patch
+
+import numpy as np
+
+from CodeEntropy.axes import AxesManager
+from tests.test_CodeEntropy.test_base import BaseTestCase
+
+
+class TestAxesManager(BaseTestCase):
+    def setUp(self):
+        super().setUp()
+
+    def test_get_residue_axes_no_bonds_custom_axes_branch(self):
+        """
+        Tests that: atom_set empty (len == 0) -> custom axes branch
+        """
+        axes_manager = AxesManager()
+        data_container = MagicMock()
+
+        atom_set = MagicMock()
+        atom_set.__len__.return_value = 0
+
+        residue = MagicMock()
+
+        data_container.select_atoms.side_effect = [atom_set, residue]
+
+        center = np.array([1.0, 2.0, 3.0])
+        residue.atoms.center_of_mass.return_value = center
+
+        UAs = MagicMock()
+        UAs.positions = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
+        residue.select_atoms.return_value = UAs
+
+        UA_masses = [12.0, 14.0]
+        axes_manager.get_UA_masses = MagicMock(return_value=UA_masses)
+
+        moi_tensor = np.eye(3) * 5.0
+        axes_manager.get_moment_of_inertia_tensor = MagicMock(return_value=moi_tensor)
+
+        rot_axes = np.eye(3)
+        moi = np.array([10.0, 9.0, 8.0])
+        axes_manager.get_custom_principal_axes = MagicMock(return_value=(rot_axes, moi))
+
+        trans_axes_out, rot_axes_out, center_out, moi_out = (
+            axes_manager.get_residue_axes(
+                data_container=data_container,
+                index=5,
+            )
+        )
+
+        calls = data_container.select_atoms.call_args_list
+        assert len(calls) == 2
+        assert calls[0].args[0] == "(resindex 4 or resindex 6) and bonded resid 5"
+        assert calls[1].args[0] == "resindex 5"
+
+        residue.select_atoms.assert_called_once_with("mass 2 to 999")
+
+        axes_manager.get_UA_masses.assert_called_once_with(residue)
+
+        axes_manager.get_moment_of_inertia_tensor.assert_called_once()
+        tensor_args, tensor_kwargs = axes_manager.get_moment_of_inertia_tensor.call_args
+        np.testing.assert_array_equal(tensor_args[0], center)
+        np.testing.assert_array_equal(tensor_args[1], UAs.positions)
+        assert tensor_args[2] == UA_masses
+        assert tensor_kwargs == {}
+
+        axes_manager.get_custom_principal_axes.assert_called_once_with(moi_tensor)
+
+        np.testing.assert_array_equal(trans_axes_out, rot_axes)
+        np.testing.assert_array_equal(rot_axes_out, rot_axes)
+        np.testing.assert_array_equal(center_out, center)
+        np.testing.assert_array_equal(moi_out, moi)
+
+    def test_get_residue_axes_bonded_default_axes_branch(self):
+        """
+        Tests that: atom_set non-empty (len != 0) -> default/bonded branch
+        """
+        axes_manager = AxesManager()
+        data_container = MagicMock()
+        data_container.atoms = MagicMock()
+
+        atom_set = MagicMock()
+        atom_set.__len__.return_value = 2
+
+        residue = MagicMock()
+        data_container.select_atoms.side_effect = [atom_set, residue]
+
+        trans_axes_expected = np.eye(3) * 2
+        data_container.atoms.principal_axes.return_value = trans_axes_expected
+
+        rot_axes_expected = np.eye(3) * 3
+        residue.principal_axes.return_value = rot_axes_expected
+
+        moi_tensor = np.eye(3)
+        residue.moment_of_inertia.return_value = moi_tensor
+
+        center_expected = np.array([9.0, 8.0, 7.0])
+        residue.atoms.center_of_mass.return_value = center_expected
+        residue.center_of_mass.return_value = center_expected
+
+        with (
+            patch("CodeEntropy.axes.make_whole", autospec=True),
+            patch.object(
+                AxesManager,
+                "get_vanilla_axes",
+                return_value=(rot_axes_expected, np.array([3.0, 2.0, 1.0])),
+            ),
+        ):
+            trans_axes_out, rot_axes_out, center_out, moi_out = (
+                axes_manager.get_residue_axes(
+                    data_container=data_container,
+                    index=5,
+                )
+            )
+
+        np.testing.assert_allclose(trans_axes_out, trans_axes_expected)
+        np.testing.assert_allclose(rot_axes_out, rot_axes_expected)
+        np.testing.assert_allclose(center_out, center_expected)
+        np.testing.assert_allclose(moi_out, np.array([3.0, 2.0, 1.0]))
+
+    def test_get_UA_axes_returns_expected_outputs(self):
+        """
+        Tests that: `get_UA_axes` returns expected UA axes.
+        """
+        axes = AxesManager()
+
+        dc = MagicMock()
+        dc.atoms = MagicMock()
+        dc.dimensions = np.array([1.0, 2.0, 3.0, 90.0, 90.0, 90.0])
+        dc.atoms.center_of_mass.return_value = np.array([0.0, 0.0, 0.0])
+
+        uas = MagicMock()
+        uas.positions = np.zeros((2, 3))
+
+        a0 = MagicMock()
+        a0.index = 7
+        a1 = MagicMock()
+        a1.index = 9
+        heavy_atoms = [a0, a1]
+
+        heavy_ag = MagicMock()
+        heavy_ag.positions = np.array([[9.9, 8.8, 7.7]])
+        heavy_ag.__getitem__.return_value = MagicMock()
+
+        dc.select_atoms.side_effect = [uas, heavy_atoms, heavy_ag]
+
+        axes.get_UA_masses = MagicMock(return_value=[1.0, 1.0])
+        axes.get_moment_of_inertia_tensor = MagicMock(return_value=np.eye(3))
+
+        trans_axes_expected = np.eye(3)
+        axes.get_custom_principal_axes = MagicMock(
+            return_value=(trans_axes_expected, np.array([1.0, 1.0, 1.0]))
+        )
+
+        rot_axes_expected = np.eye(3) * 2
+        moi_expected = np.array([3.0, 2.0, 1.0])
+        axes.get_bonded_axes = MagicMock(return_value=(rot_axes_expected, moi_expected))
+
+        trans_axes, rot_axes, center, moi = axes.get_UA_axes(dc, index=1)
+
+        np.testing.assert_array_equal(trans_axes, trans_axes_expected)
+        np.testing.assert_array_equal(rot_axes, rot_axes_expected)
+        np.testing.assert_array_equal(center, heavy_ag.positions[0])
+        np.testing.assert_array_equal(moi, moi_expected)
+
+        calls = [c.args[0] for c in dc.select_atoms.call_args_list]
+        assert calls[0] == "mass 2 to 999"
+        assert calls[1] == "prop mass > 1.1"
+        assert calls[2] == "index 9"
+
+    def test_get_bonded_axes_returns_none_for_light_atom(self):
+        """
+        Tests that: bonded axes return none for light atoms
+        """
+        axes = AxesManager()
+
+        atom = MagicMock()
+        atom.mass = 1.0
+        system = MagicMock()
+
+        out = axes.get_bonded_axes(
+            system=system, atom=atom, dimensions=np.array([1.0, 2.0, 3.0])
+        )
+        assert out is None
+
+    def test_get_bonded_axes_case2_one_heavy_zero_light(self):
+        """
+        Tests that: bonded return one heavy and zero light atoms
+        """
+        axes = AxesManager()
+
+        system = MagicMock()
+        atom = MagicMock()
+        atom.mass = 12.0
+        atom.index = 0
+        atom.position = np.array([0.0, 0.0, 0.0])
+
+        heavy0 = MagicMock()
+        heavy0.position = np.array([1.0, 0.0, 0.0])
+
+        heavy_bonded = [heavy0]
+        light_bonded = []
+
+        axes.find_bonded_atoms = MagicMock(return_value=(heavy_bonded, light_bonded))
+
+        custom_axes = np.eye(3)
+        axes.get_custom_axes = MagicMock(return_value=custom_axes)
+
+        moi = np.array([1.0, 2.0, 3.0])
+        axes.get_custom_moment_of_inertia = MagicMock(return_value=moi)
+
+        flipped_axes = np.eye(3) * 2
+        axes.get_flipped_axes = MagicMock(return_value=flipped_axes)
+
+        out_axes, out_moi = axes.get_bonded_axes(
+            system, atom, np.array([10.0, 10.0, 10.0])
+        )
+
+        np.testing.assert_array_equal(out_axes, flipped_axes)
+        np.testing.assert_array_equal(out_moi, moi)
+
+        axes.get_custom_axes.assert_called_once()
+        args, _ = axes.get_custom_axes.call_args
+        np.testing.assert_array_equal(args[0], atom.position)
+        assert len(args[1]) == 1
+        np.testing.assert_array_equal(args[1][0], heavy0.position)
+        np.testing.assert_array_equal(args[2], np.zeros(3))
+        np.testing.assert_array_equal(args[3], np.array([10.0, 10.0, 10.0]))
+
+    def test_get_bonded_axes_case3_one_heavy_with_light(self):
+        """
+        Tests that: bonded axes return one heavy with one light atom
+        """
+        axes = AxesManager()
+
+        system = MagicMock()
+        atom = MagicMock()
+        atom.mass = 12.0
+        atom.index = 0
+        atom.position = np.array([0.0, 0.0, 0.0])
+
+        heavy0 = MagicMock()
+        heavy0.position = np.array([1.0, 0.0, 0.0])
+
+        light0 = MagicMock()
+        light0.position = np.array([0.0, 1.0, 0.0])
+
+        heavy_bonded = [heavy0]
+        light_bonded = [light0]
+
+        axes.find_bonded_atoms = MagicMock(return_value=(heavy_bonded, light_bonded))
+
+        custom_axes = np.eye(3)
+        axes.get_custom_axes = MagicMock(return_value=custom_axes)
+        axes.get_custom_moment_of_inertia = MagicMock(
+            return_value=np.array([1.0, 1.0, 1.0])
+        )
+        axes.get_flipped_axes = MagicMock(return_value=custom_axes)
+
+        axes.get_bonded_axes(system, atom, np.array([10.0, 10.0, 10.0]))
+
+        axes.get_custom_axes.assert_called_once()
+        args, _ = axes.get_custom_axes.call_args
+
+        np.testing.assert_array_equal(args[2], light0.position)
+
+    def test_get_bonded_axes_case5_two_or_more_heavy(self):
+        """
+        Tests that: bonded axes return two or more heavy atoms
+        """
+        axes = AxesManager()
+
+        system = MagicMock()
+        atom = MagicMock()
+        atom.mass = 12.0
+        atom.index = 0
+        atom.position = np.array([0.0, 0.0, 0.0])
+
+        heavy0 = MagicMock()
+        heavy0.position = np.array([1.0, 0.0, 0.0])
+        heavy1 = MagicMock()
+        heavy1.position = np.array([0.0, 1.0, 0.0])
+
+        heavy_bonded = MagicMock()
+        heavy_bonded.__len__.return_value = 2
+        heavy_bonded.positions = np.array([heavy0.position, heavy1.position])
+
+        heavy_bonded.__getitem__.side_effect = lambda i: [heavy0, heavy1][i]
+
+        light_bonded = []
+
+        axes.find_bonded_atoms = MagicMock(return_value=(heavy_bonded, light_bonded))
+
+        custom_axes = np.eye(3)
+        axes.get_custom_axes = MagicMock(return_value=custom_axes)
+        axes.get_custom_moment_of_inertia = MagicMock(
+            return_value=np.array([9.0, 9.0, 9.0])
+        )
+        axes.get_flipped_axes = MagicMock(return_value=custom_axes)
+
+        axes.get_bonded_axes(system, atom, np.array([10.0, 10.0, 10.0]))
+
+        axes.get_custom_axes.assert_called_once()
+        args, _ = axes.get_custom_axes.call_args
+
+        np.testing.assert_array_equal(args[1], heavy_bonded.positions)
+        np.testing.assert_array_equal(args[2], heavy1.position)
+
+    def test_find_bonded_atoms_splits_heavy_and_h(self):
+        """
+        Tests that: Bonded atoms split into heavy and hydrogen.
+        """
+        axes = AxesManager()
+
+        system = MagicMock()
+        bonded = MagicMock()
+        heavy = MagicMock()
+        hydrogens = MagicMock()
+
+        system.select_atoms.return_value = bonded
+        bonded.select_atoms.side_effect = [heavy, hydrogens]
+
+        out_heavy, out_h = axes.find_bonded_atoms(5, system)
+
+        system.select_atoms.assert_called_once_with("bonded index 5")
+        assert bonded.select_atoms.call_args_list[0].args[0] == "mass 2 to 999"
+        assert bonded.select_atoms.call_args_list[1].args[0] == "mass 1 to 1.1"
+        assert out_heavy is heavy
+        assert out_h is hydrogens
+
+    def test_get_vector_wraps_pbc(self):
+        """
+        Tests that: The vector wraps across periodic boundary.
+        """
+        axes = AxesManager()
+
+        a = np.array([9.0, 0.0, 0.0])
+        b = np.array([1.0, 0.0, 0.0])
+        dims = np.array([10.0, 10.0, 10.0])
+
+        out = axes.get_vector(a, b, dims)
+        np.testing.assert_array_equal(out, np.array([2.0, 0.0, 0.0]))
+
+    def test_get_custom_axes_returns_unit_axes(self):
+        """
+        Tests that: `get_axes` returns normalized 3x3 axes.
+        """
+        axes = AxesManager()
+
+        a = np.zeros(3)
+        b_list = [np.array([1.0, 0.0, 0.0])]
+        c = np.array([0.0, 1.0, 0.0])
+        dims = np.array([100.0, 100.0, 100.0])
+
+        out = axes.get_custom_axes(a, b_list, c, dims)
+
+        assert out.shape == (3, 3)
+        norms = np.linalg.norm(out, axis=1)
+        np.testing.assert_allclose(norms, np.ones(3))
+
+    def test_get_custom_axes_uses_bc_vector_when_multiple_heavy_atoms(self):
+        """
+        Tests that: `get_custom_axes` uses c → b_list[0] vector when b_list has
+        ≥ 2 atoms.
+        """
+        axes = AxesManager()
+
+        a = np.zeros(3)
+        b0 = np.array([1.0, 0.0, 0.0])
+        b1 = np.array([0.0, 1.0, 0.0])
+        b_list = [b0, b1]
+        c = np.array([0.0, 0.0, 1.0])
+        dimensions = np.array([10.0, 10.0, 10.0])
+
+        # Track calls to get_vector
+        axes.get_vector = MagicMock(return_value=np.array([1.0, 0.0, 0.0]))
+
+        axes.get_custom_axes(a, b_list, c, dimensions)
+
+        # get_vector should be called
+        calls = axes.get_vector.call_args_list
+
+        # Last call must be (c, b_list[0], dimensions)
+        last_args = calls[-1].args
+        np.testing.assert_array_equal(last_args[0], c)
+        np.testing.assert_array_equal(last_args[1], b0)
+        np.testing.assert_array_equal(last_args[2], dimensions)
+
+    def test_get_custom_moment_of_inertia_len2_zeroed(self):
+        """
+        Tests that: `get_custom_moment_of_inertia` zeroes one MOI component for
+        two-atom UA.
+        """
+        axes = AxesManager()
+
+        UA = MagicMock()
+        UA.positions = np.array([[1, 0, 0], [0, 1, 0]])
+        UA.masses = np.array([12.0, 1.0])
+        UA.__len__.return_value = 2
+
+        dimensions = np.array([100.0, 100.0, 100.0])
+
+        moi = axes.get_custom_moment_of_inertia(UA, np.eye(3), np.zeros(3), dimensions)
+
+        assert moi.shape == (3,)
+        assert np.any(np.isclose(moi, 0.0))
+
+    def test_get_flipped_axes_flips_negative_dot(self):
+        """
+        Tests that: `get_flipped_axes` flips axis when dot product is negative.
+        """
+        axes = AxesManager()
+
+        UA = MagicMock()
+        atom0 = MagicMock()
+        atom0.position = np.zeros(3)
+        UA.__getitem__.return_value = atom0
+
+        axes.get_vector = MagicMock(return_value=np.array([-1.0, 0.0, 0.0]))
+
+        custom_axes = np.eye(3)
+        out = axes.get_flipped_axes(
+            UA, custom_axes, np.zeros(3), np.array([10, 10, 10])
+        )
+
+        np.testing.assert_array_equal(out[0], np.array([-1.0, 0.0, 0.0]))
+
+    def test_get_moment_of_inertia_tensor_simple(self):
+        """
+        Tests that: `get_moment_of_inertia` Computes inertia tensor correctly.
+        """
+        axes = AxesManager()
+
+        center = np.zeros(3)
+        pos = np.array([[1, 0, 0], [0, 1, 0]])
+        masses = np.array([1.0, 1.0])
+        dimensions = np.array([100.0, 100.0, 100.0])
+
+        tensor = axes.get_moment_of_inertia_tensor(center, pos, masses, dimensions)
+
+        expected = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 2]])
+        np.testing.assert_array_equal(tensor, expected)
+
+    def test_get_custom_principal_axes_flips_z(self):
+        """
+        Tests that: `get_custom_principle_axes` ensures right-handed axes.
+        """
+        axes = AxesManager()
+
+        with patch("CodeEntropy.axes.np.linalg.eig") as eig:
+            eig.return_value = (
+                np.array([3, 2, 1]),
+                np.array([[1, 0, 0], [0, 1, 0], [0, 0, -1]]),
+            )
+
+            axes_out, moi = axes.get_custom_principal_axes(np.eye(3))
+
+        np.testing.assert_array_equal(axes_out[2], np.array([0, 0, 1]))
+
+    def test_get_UA_masses_sums_hydrogens(self):
+        """
+        Tests that: `get_UA_masses` sums heavy atom with bonded hydrogens.
+        """
+        axes = AxesManager()
+
+        heavy = MagicMock(mass=12.0, index=0)
+        light = MagicMock(mass=1.0, index=1)
+
+        mol = MagicMock()
+        mol.__iter__.return_value = iter([heavy, light])
+
+        bonded = MagicMock()
+        H = MagicMock(mass=1.0)
+        mol.select_atoms.return_value = bonded
+        bonded.select_atoms.return_value = [H]
+
+        out = axes.get_UA_masses(mol)
+
+        assert out == [13.0]
diff --git a/tests/test_CodeEntropy/test_entropy.py b/tests/test_CodeEntropy/test_entropy.py
index 8fce98f..72b5466 100644
--- a/tests/test_CodeEntropy/test_entropy.py
+++ b/tests/test_CodeEntropy/test_entropy.py
@@ -8,6 +8,7 @@
 
 import MDAnalysis as mda
 import numpy as np
+import numpy.linalg as la
 import pytest
 
 import tests.data as data
@@ -83,7 +84,12 @@ def test_execute_full_workflow(self):
             return_value=(mock_reduced_atom, 3, mock_levels, mock_groups)
         )
         entropy_manager._level_manager.build_covariance_matrices = MagicMock(
-            return_value=("force_matrices", "torque_matrices", "frame_counts")
+            return_value=(
+                "force_matrices",
+                "torque_matrices",
+                "forcetorque_avg",
+                "frame_counts",
+            )
         )
         entropy_manager._dihedral_analysis.build_conformational_states = MagicMock(
             return_value=(["state_ua"], ["state_res"])
@@ -128,6 +134,7 @@ def test_execute_full_workflow(self):
             mock_groups,
             "force_matrices",
             "torque_matrices",
+            "forcetorque_avg",
             ["state_ua"],
             ["state_res"],
             "frame_counts",
@@ -173,7 +180,12 @@ def test_execute_triggers_handle_water_entropy_minimal(self):
             return_value=(MagicMock(), 3, {}, {0: [0]})
         )
         entropy_manager._level_manager.build_covariance_matrices = MagicMock(
-            return_value=("force_matrices", "torque_matrices", "frame_counts")
+            return_value=(
+                "force_matrices",
+                "torque_matrices",
+                "forcetorque_avg",
+                "frame_counts",
+            )
         )
         entropy_manager._dihedral_analysis.build_conformational_states = MagicMock(
             return_value=(["state_ua"], ["state_res"])
@@ -715,6 +727,8 @@ def test_process_vibrational_only_levels(self):
         ve = MagicMock()
         ve.vibrational_entropy_calculation.side_effect = [1.11, 2.22]
 
+        forcetorque_matrix = np.eye(6)
+
         # Run the method
         manager._process_vibrational_entropy(
             group_id=0,
@@ -724,6 +738,7 @@ def test_process_vibrational_only_levels(self):
             level="Vibrational",
             force_matrix=force_matrix,
             torque_matrix=torque_matrix,
+            forcetorque_matrix=forcetorque_matrix,
             highest=True,
         )
 
@@ -731,7 +746,7 @@ def test_process_vibrational_only_levels(self):
         df = data_logger.molecule_data
         self.assertEqual(len(df), 2)  # Transvibrational and Rovibrational
 
-        expected_types = {"Transvibrational", "Rovibrational"}
+        expected_types = {"FTmat-Transvibrational", "FTmat-Rovibrational"}
         actual_types = set(entry[2] for entry in df)
         self.assertSetEqual(actual_types, expected_types)
 
@@ -739,6 +754,77 @@ def test_process_vibrational_only_levels(self):
         self.assertIn(1.11, results)
         self.assertIn(2.22, results)
 
+    def test_process_vibrational_entropy_else_branch(self):
+        """
+        Atomic unit test for EntropyManager._process_vibrational_entropy else-branch:
+        - forcetorque_matrix is None
+        - force/torque matrices are filtered
+        - ve.vibrational_entropy_calculation called for force & torque
+        - results logged as Transvibrational/Rovibrational
+        - group label added from mol_container residues/atoms
+        """
+        manager = MagicMock()
+        manager._args = MagicMock(temperature=300)
+
+        manager._level_manager = MagicMock()
+        manager._data_logger = MagicMock()
+
+        force_matrix = np.eye(3)
+        torque_matrix = np.eye(3) * 2
+
+        filtered_force = np.eye(3) * 7
+        filtered_torque = np.eye(3) * 9
+        manager._level_manager.filter_zero_rows_columns.side_effect = [
+            filtered_force,
+            filtered_torque,
+        ]
+
+        ve = MagicMock()
+        ve.vibrational_entropy_calculation.side_effect = [1.11, 2.22]
+
+        res1 = MagicMock(resname="ALA")
+        res2 = MagicMock(resname="GLY")
+        res3 = MagicMock(resname="ALA")
+        mol_container = MagicMock()
+        mol_container.residues = [res1, res2, res3]
+        mol_container.atoms = [MagicMock(), MagicMock(), MagicMock(), MagicMock()]
+
+        EntropyManager._process_vibrational_entropy(
+            manager,
+            group_id=0,
+            mol_container=mol_container,
+            number_frames=10,
+            ve=ve,
+            level="Vibrational",
+            force_matrix=force_matrix,
+            torque_matrix=torque_matrix,
+            forcetorque_matrix=None,
+            highest=True,
+        )
+
+        filter_calls = manager._level_manager.filter_zero_rows_columns.call_args_list
+        assert len(filter_calls) == 2
+
+        np.testing.assert_array_equal(filter_calls[0].args[0], force_matrix)
+        np.testing.assert_array_equal(filter_calls[1].args[0], torque_matrix)
+
+        ve_calls = ve.vibrational_entropy_calculation.call_args_list
+        assert len(ve_calls) == 2
+
+        np.testing.assert_array_equal(ve_calls[0].args[0], filtered_force)
+        assert ve_calls[0].args[1:] == ("force", 300, True)
+
+        np.testing.assert_array_equal(ve_calls[1].args[0], filtered_torque)
+        assert ve_calls[1].args[1:] == ("torque", 300, True)
+
+        manager._data_logger.add_results_data.assert_any_call(
+            0, "Vibrational", "Transvibrational", 1.11
+        )
+        manager._data_logger.add_results_data.assert_any_call(
+            0, "Vibrational", "Rovibrational", 2.22
+        )
+        manager._data_logger.add_group_label.assert_called_once_with(0, "ALA_GLY", 3, 4)
+
     def test_compute_entropies_polymer_branch(self):
         """
         Test _compute_entropies triggers _process_vibrational_entropy for 'polymer'
@@ -789,6 +875,7 @@ def test_compute_entropies_polymer_branch(self):
             groups,
             force_matrices,
             torque_matrices,
+            force_matrices,
             states_ua,
             states_res,
             frame_counts,
@@ -942,6 +1029,8 @@ def test_compute_entropies_united_atom(self):
         universe_operations.get_molecule_container = MagicMock(return_value=mol_mock)
         manager._process_united_atom_entropy = MagicMock()
 
+        force_torque_matrices = MagicMock()
+
         ve = MagicMock()
         ce = MagicMock()
 
@@ -951,6 +1040,7 @@ def test_compute_entropies_united_atom(self):
             groups,
             force_matrices,
             torque_matrices,
+            force_torque_matrices,
             states_ua,
             states_res,
             frame_counts,
@@ -1017,6 +1107,8 @@ def test_compute_entropies_residue(self):
         manager._process_vibrational_entropy = MagicMock()
         manager._process_conformational_entropy = MagicMock()
 
+        force_torque_matrices = MagicMock()
+
         # Mock entropy calculators
         ve = MagicMock()
         ce = MagicMock()
@@ -1028,6 +1120,7 @@ def test_compute_entropies_residue(self):
             groups,
             force_matrices,
             torque_matrices,
+            force_torque_matrices,
             states_ua,
             states_res,
             frame_counts,
@@ -1048,6 +1141,7 @@ def test_compute_entropies_polymer(self):
         data_logger = DataLogger()
         group_molecules = MagicMock()
         dihedral_analysis = MagicMock()
+
         manager = EntropyManager(
             run_manager,
             args,
@@ -1066,9 +1160,10 @@ def test_compute_entropies_polymer(self):
 
         force_matrices = {"poly": {0: "force_poly"}}
         torque_matrices = {"poly": {0: "torque_poly"}}
+        force_torque_matrices = {"poly": {0: "ft_poly"}}
+
         states_ua = {}
         states_res = []
-
         frame_counts = {"poly": {(0, 0): 10}}
 
         mol_mock = MagicMock()
@@ -1085,6 +1180,7 @@ def test_compute_entropies_polymer(self):
             groups,
             force_matrices,
             torque_matrices,
+            force_torque_matrices,
             states_ua,
             states_res,
             frame_counts,
@@ -1101,6 +1197,7 @@ def test_compute_entropies_polymer(self):
             "polymer",
             force_matrices["poly"][0],
             torque_matrices["poly"][0],
+            force_torque_matrices["poly"][0],
             True,
         )
 
@@ -1521,6 +1618,108 @@ def test_vibrational_entropy_polymer_torque(self):
 
         assert S_vib == pytest.approx(48.45003266069881)
 
+    def test_vibrational_entropy_calculation_forcetorqueTRANS(self):
+        """
+        Test for matrix_type='forcetorqueTRANS':
+        - verifies S_vib_total = sum(S_components[:3])
+        """
+        run_manager = MagicMock()
+        run_manager.change_lambda_units.side_effect = lambda x: x
+        kT = 2.47e-21
+        run_manager.get_KT2J.return_value = kT
+
+        ve = VibrationalEntropy(
+            run_manager,
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+        )
+
+        orig_eigvals = la.eigvals
+        la.eigvals = lambda m: np.array(
+            [1.0] * 6
+        )  # length 6 -> 6 frequencies/components
+
+        try:
+            freqs = np.array([6.0, 5.0, 4.0, 3.0, 2.0, 1.0])
+            ve.frequency_calculation = MagicMock(return_value=freqs)
+
+            matrix = np.identity(6)
+
+            result = ve.vibrational_entropy_calculation(
+                matrix=matrix,
+                matrix_type="forcetorqueTRANS",
+                temp=298,
+                highest_level=True,
+            )
+
+            sorted_freqs = np.sort(freqs)
+            exponent = ve._PLANCK_CONST * sorted_freqs / kT
+            power_positive = np.exp(exponent)
+            power_negative = np.exp(-exponent)
+            S_components = exponent / (power_positive - 1) - np.log(1 - power_negative)
+            S_components *= ve._GAS_CONST
+
+            expected = float(np.sum(S_components[:3]))
+            self.assertAlmostEqual(result, expected, places=6)
+
+        finally:
+            la.eigvals = orig_eigvals
+
+    def test_vibrational_entropy_calculation_forcetorqueROT(self):
+        """
+        Test for matrix_type='forcetorqueROT':
+        - verifies S_vib_total = sum(S_components[3:])
+        """
+        run_manager = MagicMock()
+        run_manager.change_lambda_units.side_effect = lambda x: x
+        kT = 2.47e-21
+        run_manager.get_KT2J.return_value = kT
+
+        ve = VibrationalEntropy(
+            run_manager,
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+            MagicMock(),
+        )
+
+        orig_eigvals = la.eigvals
+        la.eigvals = lambda m: np.array([1.0] * 6)
+
+        try:
+            freqs = np.array([6.0, 5.0, 4.0, 3.0, 2.0, 1.0])
+            ve.frequency_calculation = MagicMock(return_value=freqs)
+
+            matrix = np.identity(6)
+
+            result = ve.vibrational_entropy_calculation(
+                matrix=matrix,
+                matrix_type="forcetorqueROT",
+                temp=298,
+                highest_level=True,
+            )
+
+            sorted_freqs = np.sort(freqs)
+            exponent = ve._PLANCK_CONST * sorted_freqs / kT
+            power_positive = np.exp(exponent)
+            power_negative = np.exp(-exponent)
+            S_components = exponent / (power_positive - 1) - np.log(1 - power_negative)
+            S_components *= ve._GAS_CONST
+
+            expected = float(np.sum(S_components[3:]))
+            self.assertAlmostEqual(result, expected, places=3)
+
+        finally:
+            la.eigvals = orig_eigvals
+
     def test_calculate_water_orientational_entropy(self):
         """
         Test that orientational entropy values are correctly extracted from Sorient_dict
diff --git a/tests/test_CodeEntropy/test_levels.py b/tests/test_CodeEntropy/test_levels.py
index 98a2684..eb853df 100644
--- a/tests/test_CodeEntropy/test_levels.py
+++ b/tests/test_CodeEntropy/test_levels.py
@@ -1,7 +1,8 @@
-from unittest.mock import MagicMock
+from unittest.mock import MagicMock, patch
 
 import numpy as np
 
+from CodeEntropy.axes import AxesManager
 from CodeEntropy.levels import LevelManager
 from CodeEntropy.mda_universe_operations import UniverseOperations
 from tests.test_CodeEntropy.test_base import BaseTestCase
@@ -60,12 +61,15 @@ def test_select_levels(self):
 
     def test_get_matrices(self):
         """
-        Test `get_matrices` with mocked internal methods and a simple setup.
-        Ensures that the method returns correctly shaped force and torque matrices.
+        Atomic unit test for LevelManager.get_matrices:
+        - AxesManager is mocked
+        - No inertia / MDAnalysis math
+        - Verifies block matrix construction and shape only
         """
         universe_operations = UniverseOperations()
         level_manager = LevelManager(universe_operations)
 
+        # Two beads
         bead1 = MagicMock()
         bead1.principal_axes.return_value = np.ones(3)
 
@@ -73,39 +77,62 @@ def test_get_matrices(self):
         bead2.principal_axes.return_value = np.ones(3)
 
         level_manager.get_beads = MagicMock(return_value=[bead1, bead2])
+
         level_manager.get_weighted_forces = MagicMock(
             return_value=np.array([1.0, 2.0, 3.0])
         )
         level_manager.get_weighted_torques = MagicMock(
             return_value=np.array([0.5, 1.5, 2.5])
         )
-        level_manager.create_submatrix = MagicMock(return_value=np.identity(3))
+
+        # Deterministic 3x3 submatrix for every (i,j) call
+        I3 = np.identity(3)
+        level_manager.create_submatrix = MagicMock(return_value=I3)
 
         data_container = MagicMock()
         data_container.atoms = MagicMock()
         data_container.atoms.principal_axes.return_value = np.ones(3)
 
-        force_matrix, torque_matrix = level_manager.get_matrices(
-            data_container=data_container,
-            level="residue",
-            highest_level=True,
-            force_matrix=None,
-            torque_matrix=None,
-            force_partitioning=0.5,
-        )
+        dummy_trans_axes = np.eye(3)
+        dummy_rot_axes = np.eye(3)
+        dummy_center = np.zeros(3)
+        dummy_moi = np.eye(3)
+
+        with patch("CodeEntropy.levels.AxesManager") as AxesManagerMock:
+            axes = AxesManagerMock.return_value
+            axes.get_residue_axes.return_value = (
+                dummy_trans_axes,
+                dummy_rot_axes,
+                dummy_center,
+                dummy_moi,
+            )
 
-        self.assertIsInstance(force_matrix, np.ndarray)
-        self.assertIsInstance(torque_matrix, np.ndarray)
+            force_matrix, torque_matrix = level_manager.get_matrices(
+                data_container=data_container,
+                level="residue",
+                highest_level=True,
+                force_matrix=None,
+                torque_matrix=None,
+                force_partitioning=0.5,
+                customised_axes=True,
+            )
 
-        self.assertEqual(force_matrix.shape, (6, 6))
-        self.assertEqual(torque_matrix.shape, (6, 6))
+        # Shape: 2 beads × 3 dof => 6×6
+        assert force_matrix.shape == (6, 6)
+        assert torque_matrix.shape == (6, 6)
 
-        level_manager.get_beads.assert_called_once_with(data_container, "residue")
+        # Expected block structure when every block is I3:
+        expected = np.block([[I3, I3], [I3, I3]])
+        np.testing.assert_array_equal(force_matrix, expected)
+        np.testing.assert_array_equal(torque_matrix, expected)
 
-        self.assertEqual(level_manager.get_weighted_forces.call_count, 2)
-        self.assertEqual(level_manager.get_weighted_torques.call_count, 2)
+        # Lightweight behavioral assertions
+        level_manager.get_beads.assert_called_once_with(data_container, "residue")
+        assert axes.get_residue_axes.call_count == 2
 
-        self.assertEqual(level_manager.create_submatrix.call_count, 6)
+        # For 2 beads: (0,0), (0,1), (1,1) => 3 pairs;
+        # each pair calls create_submatrix twice (force+torque)
+        assert level_manager.create_submatrix.call_count == 6
 
     def test_get_matrices_force_shape_mismatch(self):
         """
@@ -115,6 +142,7 @@ def test_get_matrices_force_shape_mismatch(self):
         universe_operations = UniverseOperations()
         level_manager = LevelManager(universe_operations)
 
+        # Two beads -> force_block will be 6x6
         bead1 = MagicMock()
         bead1.principal_axes.return_value = np.ones(3)
 
@@ -129,6 +157,7 @@ def test_get_matrices_force_shape_mismatch(self):
         level_manager.get_weighted_torques = MagicMock(
             return_value=np.array([0.5, 1.5, 2.5])
         )
+
         level_manager.create_submatrix = MagicMock(return_value=np.identity(3))
 
         data_container = MagicMock()
@@ -138,17 +167,32 @@ def test_get_matrices_force_shape_mismatch(self):
         bad_force_matrix = np.zeros((3, 3))
         correct_torque_matrix = np.zeros((6, 6))
 
-        with self.assertRaises(ValueError) as context:
-            level_manager.get_matrices(
-                data_container=data_container,
-                level="residue",
-                highest_level=True,
-                force_matrix=bad_force_matrix,
-                torque_matrix=correct_torque_matrix,
-                force_partitioning=0.5,
+        dummy_trans_axes = np.eye(3)
+        dummy_rot_axes = np.eye(3)
+        dummy_center = np.zeros(3)
+        dummy_moi = np.eye(3)
+
+        with patch("CodeEntropy.levels.AxesManager") as AxesManagerMock:
+            axes = AxesManagerMock.return_value
+            axes.get_residue_axes.return_value = (
+                dummy_trans_axes,
+                dummy_rot_axes,
+                dummy_center,
+                dummy_moi,
             )
 
-        self.assertIn("Inconsistent force matrix shape", str(context.exception))
+            with self.assertRaises(ValueError) as context:
+                level_manager.get_matrices(
+                    data_container=data_container,
+                    level="residue",
+                    highest_level=True,
+                    force_matrix=bad_force_matrix,
+                    torque_matrix=correct_torque_matrix,
+                    force_partitioning=0.5,
+                    customised_axes=True,
+                )
+
+        assert "force matrix shape" in str(context.exception)
 
     def test_get_matrices_torque_shape_mismatch(self):
         """
@@ -179,19 +223,35 @@ def test_get_matrices_torque_shape_mismatch(self):
         data_container.atoms.principal_axes.return_value = np.ones(3)
 
         correct_force_matrix = np.zeros((6, 6))
-        bad_torque_matrix = np.zeros((3, 3))  # Incorrect shape
-
-        with self.assertRaises(ValueError) as context:
-            level_manager.get_matrices(
-                data_container=data_container,
-                level="residue",
-                highest_level=True,
-                force_matrix=correct_force_matrix,
-                torque_matrix=bad_torque_matrix,
-                force_partitioning=0.5,
+        bad_torque_matrix = np.zeros((3, 3))  # Incorrect shape (should be 6x6)
+
+        # Mock AxesManager return tuple to satisfy unpacking
+        dummy_trans_axes = np.eye(3)
+        dummy_rot_axes = np.eye(3)
+        dummy_center = np.zeros(3)
+        dummy_moi = np.eye(3)
+
+        with patch("CodeEntropy.levels.AxesManager") as AxesManagerMock:
+            axes = AxesManagerMock.return_value
+            axes.get_residue_axes.return_value = (
+                dummy_trans_axes,
+                dummy_rot_axes,
+                dummy_center,
+                dummy_moi,
             )
 
-        self.assertIn("Inconsistent torque matrix shape", str(context.exception))
+            with self.assertRaises(ValueError) as context:
+                level_manager.get_matrices(
+                    data_container=data_container,
+                    level="residue",
+                    highest_level=True,
+                    force_matrix=correct_force_matrix,
+                    torque_matrix=bad_torque_matrix,
+                    force_partitioning=0.5,
+                    customised_axes=True,
+                )
+
+        assert "torque matrix shape" in str(context.exception)
 
     def test_get_matrices_torque_consistency(self):
         """
@@ -224,27 +284,458 @@ def test_get_matrices_torque_consistency(self):
         initial_force_matrix = np.zeros((6, 6))
         initial_torque_matrix = np.zeros((6, 6))
 
-        force_matrix_1, torque_matrix_1 = level_manager.get_matrices(
-            data_container=data_container,
-            level="residue",
-            highest_level=True,
-            force_matrix=initial_force_matrix.copy(),
-            torque_matrix=initial_torque_matrix.copy(),
-            force_partitioning=0.5,
-        )
+        # Mock AxesManager return tuple (unpacked by get_matrices)
+        dummy_trans_axes = np.eye(3)
+        dummy_rot_axes = np.eye(3)
+        dummy_center = np.zeros(3)
+        dummy_moi = np.eye(3)
+
+        with patch("CodeEntropy.levels.AxesManager") as AxesManagerMock:
+            axes = AxesManagerMock.return_value
+            axes.get_residue_axes.return_value = (
+                dummy_trans_axes,
+                dummy_rot_axes,
+                dummy_center,
+                dummy_moi,
+            )
 
-        force_matrix_2, torque_matrix_2 = level_manager.get_matrices(
-            data_container=data_container,
-            level="residue",
-            highest_level=True,
-            force_matrix=initial_force_matrix.copy(),
-            torque_matrix=initial_torque_matrix.copy(),
-            force_partitioning=0.5,
-        )
+            force_matrix_1, torque_matrix_1 = level_manager.get_matrices(
+                data_container=data_container,
+                level="residue",
+                highest_level=True,
+                force_matrix=initial_force_matrix.copy(),
+                torque_matrix=initial_torque_matrix.copy(),
+                force_partitioning=0.5,
+                customised_axes=True,
+            )
+
+            force_matrix_2, torque_matrix_2 = level_manager.get_matrices(
+                data_container=data_container,
+                level="residue",
+                highest_level=True,
+                force_matrix=initial_force_matrix.copy(),
+                torque_matrix=initial_torque_matrix.copy(),
+                force_partitioning=0.5,
+                customised_axes=True,
+            )
 
         np.testing.assert_array_equal(force_matrix_1, force_matrix_2)
         np.testing.assert_array_equal(torque_matrix_1, torque_matrix_2)
 
+        assert force_matrix_1.shape == (6, 6)
+        assert torque_matrix_1.shape == (6, 6)
+
+    def test_get_matrices_united_atom_customised_axes(self):
+        """
+        Test that: level='united_atom' with customised_axes=True
+        Verifies:
+        - UA axes path is taken
+        - block matrix shape is correct for 1 bead (3x3)
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        bead = MagicMock()
+        level_manager.get_beads = MagicMock(return_value=[bead])
+
+        level_manager.get_weighted_forces = MagicMock(
+            return_value=np.array([1.0, 2.0, 3.0])
+        )
+        level_manager.get_weighted_torques = MagicMock(
+            return_value=np.array([0.5, 1.5, 2.5])
+        )
+
+        I3 = np.identity(3)
+        level_manager.create_submatrix = MagicMock(return_value=I3)
+
+        data_container = MagicMock()
+        data_container.atoms = MagicMock()
+        data_container.atoms.principal_axes.return_value = np.ones(3)
+
+        dummy_trans_axes = np.eye(3)
+        dummy_rot_axes = np.eye(3)
+        dummy_center = np.zeros(3)
+        dummy_moi = np.array([1.0, 1.0, 1.0])
+
+        with patch("CodeEntropy.levels.AxesManager") as AxesManagerMock:
+            axes = AxesManagerMock.return_value
+            axes.get_UA_axes.return_value = (
+                dummy_trans_axes,
+                dummy_rot_axes,
+                dummy_center,
+                dummy_moi,
+            )
+
+            force_matrix, torque_matrix = level_manager.get_matrices(
+                data_container=data_container,
+                level="united_atom",
+                highest_level=True,
+                force_matrix=None,
+                torque_matrix=None,
+                force_partitioning=0.5,
+                customised_axes=True,
+            )
+
+        assert force_matrix.shape == (3, 3)
+        assert torque_matrix.shape == (3, 3)
+        np.testing.assert_array_equal(force_matrix, I3)
+        np.testing.assert_array_equal(torque_matrix, I3)
+
+        axes.get_UA_axes.assert_called_once()
+        assert axes.get_residue_axes.call_count == 0
+
+    def test_get_matrices_non_customised_axes_path_atomic(self):
+        """
+        Tests that `customised_axes=False` triggers the non-customised axes path.
+
+        Verifies that:
+        - translational axes are taken from `data_container.atoms.principal_axes()`
+        - rotational axes are taken from `bead.principal_axes()` (real-valued)
+        - bead moment of inertia and center of mass are queried
+        - force and torque matrices are assembled with size (3N, 3N) for N beads
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        bead1, bead2 = MagicMock(), MagicMock()
+        bead1.principal_axes.return_value = np.eye(3) * (1 + 2j)
+        bead2.principal_axes.return_value = np.eye(3) * (1 + 2j)
+        bead1.center_of_mass.return_value = np.zeros(3)
+        bead2.center_of_mass.return_value = np.zeros(3)
+        bead1.moment_of_inertia.return_value = np.eye(3)
+        bead2.moment_of_inertia.return_value = np.eye(3)
+
+        level_manager.get_beads = MagicMock(return_value=[bead1, bead2])
+        level_manager.get_weighted_forces = MagicMock(
+            return_value=np.array([1.0, 2.0, 3.0])
+        )
+        level_manager.get_weighted_torques = MagicMock(
+            return_value=np.array([0.5, 1.5, 2.5])
+        )
+        level_manager.create_submatrix = MagicMock(return_value=np.eye(3))
+
+        data_container = MagicMock()
+        data_container.atoms = MagicMock()
+        data_container.atoms.principal_axes.return_value = np.eye(3)
+
+        with (
+            patch("CodeEntropy.levels.make_whole", autospec=True),
+            patch(
+                "CodeEntropy.levels.np.linalg.eig",
+                return_value=(np.array([1.0, 3.0, 2.0]), None),
+            ),
+        ):
+            force_matrix, torque_matrix = level_manager.get_matrices(
+                data_container=data_container,
+                level="polymer",
+                highest_level=True,
+                force_matrix=None,
+                torque_matrix=None,
+                force_partitioning=0.5,
+                customised_axes=False,
+            )
+
+        data_container.atoms.principal_axes.assert_called()
+        bead1.principal_axes.assert_called()
+        bead2.principal_axes.assert_called()
+        bead1.center_of_mass.assert_called()
+        bead2.center_of_mass.assert_called()
+        bead1.moment_of_inertia.assert_called()
+        bead2.moment_of_inertia.assert_called()
+
+        assert force_matrix.shape == (6, 6)
+        assert torque_matrix.shape == (6, 6)
+
+    def test_get_matrices_accepts_existing_same_shape(self):
+        """
+        Test that: if force_matrix and torque_matrix are provided with correct shape,
+        no error is raised and returned matrices match the newly computed blocks.
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        bead1 = MagicMock()
+        bead2 = MagicMock()
+        level_manager.get_beads = MagicMock(return_value=[bead1, bead2])
+
+        level_manager.get_weighted_forces = MagicMock(
+            return_value=np.array([1.0, 2.0, 3.0])
+        )
+        level_manager.get_weighted_torques = MagicMock(
+            return_value=np.array([0.5, 1.5, 2.5])
+        )
+
+        I3 = np.identity(3)
+        level_manager.create_submatrix = MagicMock(return_value=I3)
+
+        data_container = MagicMock()
+        data_container.atoms = MagicMock()
+        data_container.atoms.principal_axes.return_value = np.ones(3)
+
+        dummy_trans_axes = np.eye(3)
+        dummy_rot_axes = np.eye(3)
+        dummy_center = np.zeros(3)
+        dummy_moi = np.array([1.0, 1.0, 1.0])
+
+        existing_force = np.zeros((6, 6))
+        existing_torque = np.zeros((6, 6))
+
+        with patch("CodeEntropy.levels.AxesManager") as AxesManagerMock:
+            axes = AxesManagerMock.return_value
+            axes.get_residue_axes.return_value = (
+                dummy_trans_axes,
+                dummy_rot_axes,
+                dummy_center,
+                dummy_moi,
+            )
+
+            force_matrix, torque_matrix = level_manager.get_matrices(
+                data_container=data_container,
+                level="residue",
+                highest_level=True,
+                force_matrix=existing_force,
+                torque_matrix=existing_torque,
+                force_partitioning=0.5,
+                customised_axes=True,
+            )
+
+        expected = np.block([[I3, I3], [I3, I3]])
+        np.testing.assert_array_equal(force_matrix, expected)
+        np.testing.assert_array_equal(torque_matrix, expected)
+
+    def test_get_combined_forcetorque_matrices_residue_customised_init(self):
+        """
+        Test: level='residue', customised_axes=True uses AxesManager.get_residue_axes
+        and returns a (6N x 6N) block matrix for N beads.
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        bead1 = MagicMock()
+        bead2 = MagicMock()
+        level_manager.get_beads = MagicMock(return_value=[bead1, bead2])
+
+        wf = np.array([1.0, 2.0, 3.0])
+        wt = np.array([4.0, 5.0, 6.0])
+        level_manager.get_weighted_forces = MagicMock(return_value=wf)
+        level_manager.get_weighted_torques = MagicMock(return_value=wt)
+
+        I6 = np.identity(6)
+        level_manager.create_FTsubmatrix = MagicMock(return_value=I6)
+
+        data_container = MagicMock()
+        data_container.atoms = MagicMock()
+
+        dummy_trans_axes = np.eye(3)
+        dummy_rot_axes = np.eye(3)
+        dummy_center = np.zeros(3)
+        dummy_moi = np.array([1.0, 1.0, 1.0])
+
+        with patch("CodeEntropy.levels.AxesManager") as AxesManagerMock:
+            axes = AxesManagerMock.return_value
+            axes.get_residue_axes.return_value = (
+                dummy_trans_axes,
+                dummy_rot_axes,
+                dummy_center,
+                dummy_moi,
+            )
+
+            ft_matrix = level_manager.get_combined_forcetorque_matrices(
+                data_container=data_container,
+                level="residue",
+                highest_level=True,
+                forcetorque_matrix=None,
+                force_partitioning=0.5,
+                customised_axes=True,
+            )
+
+        assert ft_matrix.shape == (12, 12)
+
+        expected = np.block([[I6, I6], [I6, I6]])
+        np.testing.assert_array_equal(ft_matrix, expected)
+
+        assert axes.get_residue_axes.call_count == 2
+        assert level_manager.create_FTsubmatrix.call_count == 3
+
+    def test_get_combined_forcetorque_matrices_noncustomised_axes_path(self):
+        """
+        Test that: customised_axes=False forces else-path:
+        - make_whole(data_container.atoms) and make_whole(bead) called
+        - trans_axes = data_container.atoms.principal_axes()
+        - rot_axes, moment_of_inertia = AxesManager.get_vanilla_axes(bead)
+        - center = bead.center_of_mass(unwrap=True)
+        - FT block matrix assembled via create_FTsubmatrix and np.block
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        bead1 = MagicMock(name="bead1")
+        bead2 = MagicMock(name="bead2")
+        beads = [bead1, bead2]
+
+        level_manager.get_beads = MagicMock(return_value=beads)
+
+        data_container = MagicMock(name="data_container")
+        data_container.atoms = MagicMock(name="atoms")
+        data_container.atoms.principal_axes.return_value = np.eye(3)
+
+        # Forces/torques are 3-vectors -> concatenated to length 6
+        level_manager.get_weighted_forces = MagicMock(
+            side_effect=[
+                np.array([1.0, 2.0, 3.0]),
+                np.array([1.1, 2.1, 3.1]),
+            ]
+        )
+        level_manager.get_weighted_torques = MagicMock(
+            side_effect=[
+                np.array([4.0, 5.0, 6.0]),
+                np.array([4.1, 5.1, 6.1]),
+            ]
+        )
+
+        level_manager.create_FTsubmatrix = MagicMock(return_value=np.identity(6))
+
+        rot_axes_expected = np.eye(3)
+        moi_expected = np.array([3.0, 2.0, 1.0])
+
+        with (
+            patch("CodeEntropy.levels.make_whole", autospec=True) as mw_mock,
+            patch(
+                "CodeEntropy.axes.AxesManager.get_vanilla_axes",
+                autospec=True,
+                return_value=(rot_axes_expected, moi_expected),
+            ) as vanilla_mock,
+        ):
+            bead1.center_of_mass.return_value = np.zeros(3)
+            bead2.center_of_mass.return_value = np.zeros(3)
+
+            ft_matrix = level_manager.get_combined_forcetorque_matrices(
+                data_container=data_container,
+                level="polymer",
+                highest_level=True,
+                forcetorque_matrix=None,
+                force_partitioning=0.5,
+                customised_axes=False,
+            )
+
+        data_container.atoms.principal_axes.assert_called()
+        bead1.center_of_mass.assert_called_with(unwrap=True)
+        bead2.center_of_mass.assert_called_with(unwrap=True)
+
+        assert vanilla_mock.call_count == 2  # once per bead
+
+        # make_whole is called twice per bead: on data_container.atoms and on bead
+        assert mw_mock.call_count == 4
+        mw_mock.assert_any_call(data_container.atoms)
+        mw_mock.assert_any_call(bead1)
+        mw_mock.assert_any_call(bead2)
+
+        # result shape: (6N, 6N) with N=2
+        assert ft_matrix.shape == (12, 12)
+
+    def test_get_combined_forcetorque_matrices_shape_mismatch_raises(self):
+        """
+        Test that: raises ValueError when existing forcetorque_matrix has wrong shape.
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        bead1 = MagicMock()
+        bead2 = MagicMock()
+        level_manager.get_beads = MagicMock(return_value=[bead1, bead2])
+
+        level_manager.get_weighted_forces = MagicMock(
+            return_value=np.array([1.0, 2.0, 3.0])
+        )
+        level_manager.get_weighted_torques = MagicMock(
+            return_value=np.array([4.0, 5.0, 6.0])
+        )
+        level_manager.create_FTsubmatrix = MagicMock(return_value=np.identity(6))
+
+        data_container = MagicMock()
+        data_container.atoms = MagicMock()
+
+        dummy_trans_axes = np.eye(3)
+        dummy_rot_axes = np.eye(3)
+        dummy_center = np.zeros(3)
+        dummy_moi = np.array([1.0, 1.0, 1.0])
+
+        bad_existing = np.zeros((6, 6))
+
+        with patch("CodeEntropy.levels.AxesManager") as AxesManagerMock:
+            axes = AxesManagerMock.return_value
+            axes.get_residue_axes.return_value = (
+                dummy_trans_axes,
+                dummy_rot_axes,
+                dummy_center,
+                dummy_moi,
+            )
+
+            with self.assertRaises(ValueError) as ctx:
+                level_manager.get_combined_forcetorque_matrices(
+                    data_container=data_container,
+                    level="residue",
+                    highest_level=True,
+                    forcetorque_matrix=bad_existing,
+                    force_partitioning=0.5,
+                    customised_axes=True,
+                )
+
+        assert "forcetorque matrix shape" in str(ctx.exception)
+
+    def test_get_combined_forcetorque_matrices_existing_same_shape(self):
+        """
+        Test that: if existing forcetorque_matrix has correct shape, function returns
+        the newly computed block (no ValueError).
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        bead1 = MagicMock()
+        bead2 = MagicMock()
+        level_manager.get_beads = MagicMock(return_value=[bead1, bead2])
+
+        level_manager.get_weighted_forces = MagicMock(
+            return_value=np.array([1.0, 2.0, 3.0])
+        )
+        level_manager.get_weighted_torques = MagicMock(
+            return_value=np.array([4.0, 5.0, 6.0])
+        )
+
+        I6 = np.identity(6)
+        level_manager.create_FTsubmatrix = MagicMock(return_value=I6)
+
+        data_container = MagicMock()
+        data_container.atoms = MagicMock()
+
+        dummy_trans_axes = np.eye(3)
+        dummy_rot_axes = np.eye(3)
+        dummy_center = np.zeros(3)
+        dummy_moi = np.array([1.0, 1.0, 1.0])
+
+        existing_ok = np.zeros((12, 12))
+
+        with patch("CodeEntropy.levels.AxesManager") as AxesManagerMock:
+            axes = AxesManagerMock.return_value
+            axes.get_residue_axes.return_value = (
+                dummy_trans_axes,
+                dummy_rot_axes,
+                dummy_center,
+                dummy_moi,
+            )
+
+            ft_matrix = level_manager.get_combined_forcetorque_matrices(
+                data_container=data_container,
+                level="residue",
+                highest_level=True,
+                forcetorque_matrix=existing_ok,
+                force_partitioning=0.5,
+                customised_axes=True,
+            )
+
+        expected = np.block([[I6, I6], [I6, I6]])
+        np.testing.assert_array_equal(ft_matrix, expected)
+
     def test_get_beads_polymer_level(self):
         """
         Test `get_beads` for 'polymer' level.
@@ -453,70 +944,75 @@ def test_get_weighted_forces_negative_mass_raises_value_error(self):
 
     def test_get_weighted_torques_weighted_torque_basic(self):
         """
-        Test basic torque calculation with non-zero moment of inertia and torques.
+        Test basic weighted torque calculation for a single-atom bead.
+
+        Setup:
+        r = [1, 0, 0], F = [0, 1, 0]  => r x F = [0, 0, 1]
+        With force_partitioning=0.5, rot_axes=I, MOI=[1,1,1],
+        expected weighted torque is [0, 0, 0.5].
         """
         universe_operations = UniverseOperations()
-
         level_manager = LevelManager(universe_operations)
+        axes_manager = AxesManager()
 
-        # Mock atom
-        atom = MagicMock()
-        atom.index = 0
-
-        # Mock bead
         bead = MagicMock()
-        bead.atoms = [atom]
-        bead.center_of_mass.return_value = np.array([0.0, 0.0, 0.0])
-        bead.moment_of_inertia.return_value = np.identity(3)
-
-        # Mock data_container
-        data_container = MagicMock()
-        data_container.atoms.__getitem__.return_value.position = np.array(
-            [1.0, 0.0, 0.0]
-        )
-        data_container.atoms.__getitem__.return_value.force = np.array([0.0, 1.0, 0.0])
-
-        # Rotation axes (identity matrix)
-        rot_axes = np.identity(3)
+        bead.positions = np.array([[1.0, 0.0, 0.0]])
+        bead.forces = np.array([[0.0, 1.0, 0.0]])
+        bead.dimensions = np.array([10.0, 10.0, 10.0])
 
+        rot_axes = np.eye(3)
+        center = np.zeros(3)
         force_partitioning = 0.5
+        moment_of_inertia = np.array([1.0, 1.0, 1.0])
+
+        with patch.object(
+            AxesManager, "get_vector", return_value=bead.positions - center
+        ) as gv_mock:
+            result = level_manager.get_weighted_torques(
+                bead=bead,
+                rot_axes=rot_axes,
+                center=center,
+                force_partitioning=force_partitioning,
+                moment_of_inertia=moment_of_inertia,
+                axes_manager=axes_manager,
+            )
 
-        result = level_manager.get_weighted_torques(
-            data_container, bead, rot_axes, force_partitioning
-        )
+        gv_mock.assert_called()
 
-        np.testing.assert_array_almost_equal(result, np.array([0.0, 0.0, 0.5]))
+        expected = np.array([0.0, 0.0, 0.5])
+        np.testing.assert_allclose(result, expected)
 
     def test_get_weighted_torques_zero_torque_skips_division(self):
         """
         Test that zero torque components skip division and remain zero.
         """
         universe_operations = UniverseOperations()
-
         level_manager = LevelManager(universe_operations)
-
-        atom = MagicMock()
-        atom.index = 0
+        axes_manager = AxesManager()
 
         bead = MagicMock()
-        bead.atoms = [atom]
-        bead.center_of_mass.return_value = np.array([0.0, 0.0, 0.0])
-        bead.moment_of_inertia.return_value = np.identity(3)
-
-        data_container = MagicMock()
-        data_container.atoms.__getitem__.return_value.position = np.array(
-            [0.0, 0.0, 0.0]
-        )
-        data_container.atoms.__getitem__.return_value.force = np.array([0.0, 0.0, 0.0])
+        bead.positions = np.array([[0.0, 0.0, 0.0]])
+        bead.forces = np.array([[0.0, 0.0, 0.0]])
+        bead.dimensions = np.array([10.0, 10.0, 10.0])
 
         rot_axes = np.identity(3)
-
+        center = np.array([0.0, 0.0, 0.0])
         force_partitioning = 0.5
+        moment_of_inertia = np.array([1.0, 2.0, 3.0])
+
+        with patch.object(
+            AxesManager, "get_vector", return_value=bead.positions - center
+        ):
+            result = level_manager.get_weighted_torques(
+                bead=bead,
+                rot_axes=rot_axes,
+                center=center,
+                force_partitioning=force_partitioning,
+                moment_of_inertia=moment_of_inertia,
+                axes_manager=axes_manager,
+            )
 
-        result = level_manager.get_weighted_torques(
-            data_container, bead, rot_axes, force_partitioning
-        )
-        np.testing.assert_array_almost_equal(result, np.zeros(3))
+        np.testing.assert_array_equal(result, np.zeros(3))
 
     def test_get_weighted_torques_zero_moi(self):
         """
@@ -524,78 +1020,65 @@ def test_get_weighted_torques_zero_moi(self):
         and torque in that dimension is non-zero.
         """
         universe_operations = UniverseOperations()
-
         level_manager = LevelManager(universe_operations)
-
-        atom = MagicMock()
-        atom.index = 0
+        axes_manager = AxesManager()
 
         bead = MagicMock()
-        bead.atoms = [atom]
-        bead.center_of_mass.return_value = np.array([0.0, 0.0, 0.0])
-
-        # Set moment of inertia with zero in dimension 2
-        moi = np.zeros((3, 3))
-        moi[2, 2] = 0.0
-        bead.moment_of_inertia.return_value = moi
-
-        data_container = MagicMock()
-        # Position and force that will produce a non-zero torque in z (dimension 2)
-        data_container.atoms.__getitem__.return_value.position = np.array(
-            [1.0, 0.0, 0.0]
-        )
-        data_container.atoms.__getitem__.return_value.force = np.array([0.0, 1.0, 0.0])
+        bead.positions = np.array([[1.0, 0.0, 0.0]])
+        bead.forces = np.array([[0.0, 1.0, 0.0]])
+        bead.dimensions = np.array([10.0, 10.0, 10.0])
 
         rot_axes = np.identity(3)
-
+        center = np.array([0.0, 0.0, 0.0])
         force_partitioning = 0.5
+        moment_of_inertia = np.array([1.0, 1.0, 0.0])
+
+        with patch.object(
+            AxesManager, "get_vector", return_value=bead.positions - center
+        ):
+            torque = level_manager.get_weighted_torques(
+                bead=bead,
+                rot_axes=rot_axes,
+                center=center,
+                force_partitioning=force_partitioning,
+                moment_of_inertia=moment_of_inertia,
+                axes_manager=axes_manager,
+            )
 
-        torque = level_manager.get_weighted_torques(
-            data_container, bead, rot_axes, force_partitioning
-        )
-
-        self.assertEqual(torque[2], 0)
+        np.testing.assert_array_equal(torque, np.zeros(3))
 
-    def test_get_weighted_torques_negative_moi_raises(self):
+    def test_get_weighted_torques_negative_moi_sets_zero(self):
         """
-        Should raise ValueError when moment of inertia is negative in a dimension
-        and torque in that dimension is non-zero.
+        Negative moment of inertia components should be skipped and set to 0
+        even if the corresponding torque component is non-zero.
         """
         universe_operations = UniverseOperations()
-
         level_manager = LevelManager(universe_operations)
-
-        atom = MagicMock()
-        atom.index = 0
+        axes_manager = AxesManager()
 
         bead = MagicMock()
-        bead.atoms = [atom]
-        bead.center_of_mass.return_value = np.array([0.0, 0.0, 0.0])
-
-        # Set moment of inertia with negative value in dimension 2
-        moi = np.identity(3)
-        moi *= -1.0
-        bead.moment_of_inertia.return_value = moi
-
-        data_container = MagicMock()
-        # Position and force that will produce a non-zero torque in z (dimension 2)
-        data_container.atoms.__getitem__.return_value.position = np.array(
-            [1.0, 0.0, 0.0]
-        )
-        data_container.atoms.__getitem__.return_value.force = np.array([0.0, 1.0, 0.0])
+        bead.positions = np.array([[1.0, 0.0, 0.0]])
+        bead.forces = np.array([[0.0, 1.0, 0.0]])
+        bead.dimensions = np.array([10.0, 10.0, 10.0])
 
         rot_axes = np.identity(3)
-
-        foce_partitioning = 0.5
-
-        with self.assertRaises(ValueError) as context:
-            level_manager.get_weighted_torques(
-                data_container, bead, rot_axes, foce_partitioning
+        center = np.array([0.0, 0.0, 0.0])
+        force_partitioning = 0.5
+        moment_of_inertia = np.array([1.0, 1.0, -1.0])
+
+        with patch.object(
+            AxesManager, "get_vector", return_value=bead.positions - center
+        ):
+            result = level_manager.get_weighted_torques(
+                bead=bead,
+                rot_axes=rot_axes,
+                center=center,
+                force_partitioning=force_partitioning,
+                moment_of_inertia=moment_of_inertia,
+                axes_manager=axes_manager,
             )
 
-        self.assertIn(
-            "Negative value encountered for moment of inertia", str(context.exception)
-        )
+        np.testing.assert_array_equal(result, np.zeros(3))
 
     def test_create_submatrix_basic_outer_product(self):
         """
@@ -656,22 +1139,80 @@ def test_create_submatrix_symmetric_result_when_data_equal(self):
 
         self.assertTrue(np.allclose(result, result.T))  # Check symmetry
 
-    def test_build_covariance_matrices_atomic(self):
+    def test_create_FTsubmatrix_basic_outer_product(self):
         """
-        Test `build_covariance_matrices` to ensure it correctly orchestrates
-        calls and returns dictionaries with the expected structure.
+        Test that:
+        - create_FTsubmatrix returns the outer product of two 6D vectors
+        - shape is (6, 6)
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        data_i = np.array([1, 2, 3, 4, 5, 6], dtype=float)
+        data_j = np.array([6, 5, 4, 3, 2, 1], dtype=float)
+
+        result = level_manager.create_FTsubmatrix(data_i, data_j)
+
+        expected = np.outer(data_i, data_j)
+
+        assert result.shape == (6, 6)
+        np.testing.assert_array_equal(result, expected)
 
-        This test mocks dependencies including the entropy_manager, reduced_atom
-        trajectory, levels, groups, and internal method
-        `update_force_torque_matrices`.
+    def test_create_FTsubmatrix_zero_input(self):
         """
+        Test that:
+        - if either input is zero, the result is a zero matrix
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        data_i = np.zeros(6)
+        data_j = np.array([1, 2, 3, 4, 5, 6], dtype=float)
+
+        result = level_manager.create_FTsubmatrix(data_i, data_j)
+
+        np.testing.assert_array_equal(result, np.zeros((6, 6)))
+
+    def test_create_FTsubmatrix_transpose_property(self):
+        """
+        Test that:
+        - outer(i, j).T == outer(j, i)
+        - required by block-matrix symmetry logic
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        data_i = np.arange(1, 7, dtype=float)
+        data_j = np.arange(7, 13, dtype=float)
 
-        # Instantiate your class (replace YourClass with actual class name)
+        sub_ij = level_manager.create_FTsubmatrix(data_i, data_j)
+        sub_ji = level_manager.create_FTsubmatrix(data_j, data_i)
+
+        np.testing.assert_array_equal(sub_ij.T, sub_ji)
+
+    def test_create_FTsubmatrix_dtype(self):
+        """
+        Test that:
+        - output dtype follows NumPy outer-product rules
+        """
         universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        data_i = np.ones(6, dtype=np.float64)
+        data_j = np.ones(6, dtype=np.float64)
+
+        result = level_manager.create_FTsubmatrix(data_i, data_j)
 
+        assert result.dtype == np.float64
+
+    def test_build_covariance_matrices_atomic(self):
+        """
+        Test `build_covariance_matrices` to ensure it correctly orchestrates
+        calls and returns dictionaries with the expected structure.
+        """
+        universe_operations = UniverseOperations()
         level_manager = LevelManager(universe_operations)
 
-        # Mock entropy_manager and get_molecule_container
         entropy_manager = MagicMock()
 
         # Fake atom with minimal attributes
@@ -680,30 +1221,25 @@ def test_build_covariance_matrices_atomic(self):
         atom.resid = 1
         atom.segid = "A"
 
-        # Fake molecule with atoms list
         fake_mol = MagicMock()
         fake_mol.atoms = [atom]
 
-        # Always return fake_mol from get_molecule_container
         universe_operations.get_molecule_container = MagicMock(return_value=fake_mol)
 
-        # Mock reduced_atom with trajectory yielding two timesteps
         timestep1 = MagicMock()
         timestep1.frame = 0
         timestep2 = MagicMock()
         timestep2.frame = 1
+
         reduced_atom = MagicMock()
         reduced_atom.trajectory.__getitem__.return_value = [timestep1, timestep2]
 
-        # Setup groups and levels dictionaries
         groups = {"ua": ["mol1", "mol2"]}
         levels = {"mol1": ["level1", "level2"], "mol2": ["level1"]}
 
-        # Mock update_force_torque_matrices to just track calls
         level_manager.update_force_torque_matrices = MagicMock()
 
-        # Call the method under test
-        force_matrices, torque_matrices, _ = level_manager.build_covariance_matrices(
+        force_matrices, torque_matrices, *_ = level_manager.build_covariance_matrices(
             entropy_manager=entropy_manager,
             reduced_atom=reduced_atom,
             levels=levels,
@@ -713,24 +1249,21 @@ def test_build_covariance_matrices_atomic(self):
             step=1,
             number_frames=2,
             force_partitioning=0.5,
+            combined_forcetorque=False,
+            customised_axes=True,
         )
 
-        # Assert returned matrices are dictionaries with correct keys
         self.assertIsInstance(force_matrices, dict)
         self.assertIsInstance(torque_matrices, dict)
         self.assertSetEqual(set(force_matrices.keys()), {"ua", "res", "poly"})
         self.assertSetEqual(set(torque_matrices.keys()), {"ua", "res", "poly"})
 
-        # Assert 'res' and 'poly' entries are lists of correct length
         self.assertIsInstance(force_matrices["res"], list)
         self.assertIsInstance(force_matrices["poly"], list)
         self.assertEqual(len(force_matrices["res"]), len(groups))
         self.assertEqual(len(force_matrices["poly"]), len(groups))
 
-        # Check get_molecule_container call count: 2 timesteps * 2 molecules = 4 calls
         self.assertEqual(universe_operations.get_molecule_container.call_count, 4)
-
-        # Check update_force_torque_matrices call count:
         self.assertEqual(level_manager.update_force_torque_matrices.call_count, 6)
 
     def test_update_force_torque_matrices_united_atom(self):
@@ -771,6 +1304,7 @@ def test_update_force_torque_matrices_united_atom(self):
 
         force_avg = {"ua": {}, "res": [None], "poly": [None]}
         torque_avg = {"ua": {}, "res": [None], "poly": [None]}
+        forcetorque_avg = {"ua": {}, "res": [None], "poly": [None]}
         frame_counts = {"ua": {}, "res": [None], "poly": [None]}
 
         level_manager.update_force_torque_matrices(
@@ -783,8 +1317,11 @@ def test_update_force_torque_matrices_united_atom(self):
             num_frames=10,
             force_avg=force_avg,
             torque_avg=torque_avg,
+            forcetorque_avg=forcetorque_avg,
             frame_counts=frame_counts,
             force_partitioning=0.5,
+            combined_forcetorque=False,
+            customised_axes=True,
         )
 
         assert (0, 0) in force_avg["ua"]
@@ -800,6 +1337,8 @@ def test_update_force_torque_matrices_united_atom(self):
         assert frame_counts["ua"][(0, 0)] == 1
         assert frame_counts["ua"][(0, 1)] == 1
 
+        assert forcetorque_avg["ua"] == {}
+
     def test_update_force_torque_matrices_united_atom_increment(self):
         """
         Test that update_force_torque_matrices() correctly updates (increments)
@@ -824,7 +1363,6 @@ def test_update_force_torque_matrices_united_atom_increment(self):
         selected_atoms = MagicMock()
         selected_atoms.trajectory = MagicMock()
         selected_atoms.trajectory.__getitem__.return_value = None
-
         universe_operations.new_U_select_atom.return_value = selected_atoms
 
         f_mat_1 = np.array([[1.0]])
@@ -833,10 +1371,13 @@ def test_update_force_torque_matrices_united_atom_increment(self):
         f_mat_2 = np.array([[3.0]])
         t_mat_2 = np.array([[4.0]])
 
-        level_manager.get_matrices = MagicMock(return_value=(f_mat_1, t_mat_1))
+        level_manager.get_matrices = MagicMock(
+            side_effect=[(f_mat_1, t_mat_1), (f_mat_2, t_mat_2)]
+        )
 
         force_avg = {"ua": {}, "res": [None], "poly": [None]}
         torque_avg = {"ua": {}, "res": [None], "poly": [None]}
+        forcetorque_avg = {"ua": {}, "res": [None], "poly": [None]}
         frame_counts = {"ua": {}, "res": [None], "poly": [None]}
 
         level_manager.update_force_torque_matrices(
@@ -849,12 +1390,14 @@ def test_update_force_torque_matrices_united_atom_increment(self):
             num_frames=10,
             force_avg=force_avg,
             torque_avg=torque_avg,
+            forcetorque_avg=forcetorque_avg,
             frame_counts=frame_counts,
             force_partitioning=0.5,
+            combined_forcetorque=False,
+            customised_axes=True,
         )
 
-        level_manager.get_matrices = MagicMock(return_value=(f_mat_2, t_mat_2))
-
+        # Second update
         level_manager.update_force_torque_matrices(
             entropy_manager=entropy_manager,
             mol=mol,
@@ -865,8 +1408,11 @@ def test_update_force_torque_matrices_united_atom_increment(self):
             num_frames=10,
             force_avg=force_avg,
             torque_avg=torque_avg,
+            forcetorque_avg=forcetorque_avg,
             frame_counts=frame_counts,
             force_partitioning=0.5,
+            combined_forcetorque=False,
+            customised_axes=True,
         )
 
         expected_force = f_mat_1 + (f_mat_2 - f_mat_1) / 2
@@ -874,7 +1420,9 @@ def test_update_force_torque_matrices_united_atom_increment(self):
 
         np.testing.assert_array_almost_equal(force_avg["ua"][(0, 0)], expected_force)
         np.testing.assert_array_almost_equal(torque_avg["ua"][(0, 0)], expected_torque)
-        self.assertEqual(frame_counts["ua"][(0, 0)], 2)
+        assert frame_counts["ua"][(0, 0)] == 2
+
+        assert forcetorque_avg["ua"] == {}
 
     def test_update_force_torque_matrices_residue(self):
         """
@@ -883,8 +1431,8 @@ def test_update_force_torque_matrices_residue(self):
         incrementing frame counts.
         """
         universe_operations = UniverseOperations()
-
         level_manager = LevelManager(universe_operations)
+
         entropy_manager = MagicMock()
         mol = MagicMock()
         mol.trajectory.__getitem__.return_value = None
@@ -895,6 +1443,7 @@ def test_update_force_torque_matrices_residue(self):
 
         force_avg = {"ua": {}, "res": [None], "poly": [None]}
         torque_avg = {"ua": {}, "res": [None], "poly": [None]}
+        forcetorque_avg = {"ua": {}, "res": [None], "poly": [None]}
         frame_counts = {"ua": {}, "res": [None], "poly": [None]}
 
         level_manager.update_force_torque_matrices(
@@ -907,13 +1456,18 @@ def test_update_force_torque_matrices_residue(self):
             num_frames=10,
             force_avg=force_avg,
             torque_avg=torque_avg,
+            forcetorque_avg=forcetorque_avg,
             frame_counts=frame_counts,
             force_partitioning=0.5,
+            combined_forcetorque=False,
+            customised_axes=True,
         )
 
         np.testing.assert_array_equal(force_avg["res"][0], f_mat_mock)
         np.testing.assert_array_equal(torque_avg["res"][0], t_mat_mock)
-        self.assertEqual(frame_counts["res"][0], 1)
+        assert frame_counts["res"][0] == 1
+
+        assert forcetorque_avg["res"][0] is None
 
     def test_update_force_torque_matrices_incremental_average(self):
         """
@@ -923,8 +1477,8 @@ def test_update_force_torque_matrices_incremental_average(self):
         Ensures that float precision is maintained and no casting errors occur.
         """
         universe_operations = UniverseOperations()
-
         level_manager = LevelManager(universe_operations)
+
         entropy_manager = MagicMock()
         mol = MagicMock()
         mol.trajectory.__getitem__.return_value = None
@@ -941,6 +1495,7 @@ def test_update_force_torque_matrices_incremental_average(self):
 
         force_avg = {"ua": {}, "res": [None], "poly": [None]}
         torque_avg = {"ua": {}, "res": [None], "poly": [None]}
+        forcetorque_avg = {"ua": {}, "res": [None], "poly": [None]}
         frame_counts = {"ua": {}, "res": [None], "poly": [None]}
 
         # First update
@@ -954,8 +1509,11 @@ def test_update_force_torque_matrices_incremental_average(self):
             num_frames=10,
             force_avg=force_avg,
             torque_avg=torque_avg,
+            forcetorque_avg=forcetorque_avg,
             frame_counts=frame_counts,
             force_partitioning=0.5,
+            combined_forcetorque=False,
+            customised_axes=True,
         )
 
         # Second update
@@ -969,8 +1527,11 @@ def test_update_force_torque_matrices_incremental_average(self):
             num_frames=10,
             force_avg=force_avg,
             torque_avg=torque_avg,
+            forcetorque_avg=forcetorque_avg,
             frame_counts=frame_counts,
             force_partitioning=0.5,
+            combined_forcetorque=False,
+            customised_axes=True,
         )
 
         expected_force = f_mat_1 + (f_mat_2 - f_mat_1) / 2
@@ -978,7 +1539,140 @@ def test_update_force_torque_matrices_incremental_average(self):
 
         np.testing.assert_array_almost_equal(force_avg["res"][0], expected_force)
         np.testing.assert_array_almost_equal(torque_avg["res"][0], expected_torque)
-        self.assertEqual(frame_counts["res"][0], 2)
+
+        assert frame_counts["res"][0] == 2
+        assert forcetorque_avg["res"][0] is None
+
+    def test_update_force_torque_matrices_residue_combined_ft_init(self):
+        """
+        Test that: When highest=True and combined_forcetorque=True at residue level,
+        update_force_torque_matrices should:
+        - call get_combined_forcetorque_matrices
+        - store ft_mat into forcetorque_avg['res'][group_id]
+        - set frame_counts['res'][group_id] = 1
+        - NOT touch force_avg/torque_avg
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        mol = MagicMock()
+        mol.trajectory.__getitem__.return_value = None
+
+        ft_mat = np.array([[1.0, 2.0], [3.0, 4.0]], dtype=float)
+
+        level_manager.get_combined_forcetorque_matrices = MagicMock(return_value=ft_mat)
+        level_manager.get_matrices = MagicMock()
+
+        force_avg = {"ua": {}, "res": [None], "poly": [None]}
+        torque_avg = {"ua": {}, "res": [None], "poly": [None]}
+        forcetorque_avg = {"ua": {}, "res": [None], "poly": [None]}
+        frame_counts = {"ua": {}, "res": [None], "poly": [None]}
+
+        level_manager.update_force_torque_matrices(
+            entropy_manager=MagicMock(),
+            mol=mol,
+            group_id=0,
+            level="residue",
+            level_list=["united_atom", "residue"],
+            time_index=0,
+            num_frames=10,
+            force_avg=force_avg,
+            torque_avg=torque_avg,
+            forcetorque_avg=forcetorque_avg,
+            frame_counts=frame_counts,
+            force_partitioning=0.5,
+            combined_forcetorque=True,
+            customised_axes=True,
+        )
+
+        level_manager.get_combined_forcetorque_matrices.assert_called_once()
+        args = level_manager.get_combined_forcetorque_matrices.call_args.args
+        assert args[0] is mol
+        assert args[1] == "residue"
+        assert args[2] is True
+        assert args[3] is None
+        assert args[4] == 0.5
+        assert args[5] is True
+
+        np.testing.assert_array_equal(forcetorque_avg["res"][0], ft_mat)
+        assert frame_counts["res"][0] == 1
+
+        level_manager.get_matrices.assert_not_called()
+
+        assert force_avg["res"][0] is None
+        assert torque_avg["res"][0] is None
+
+    def test_update_force_torque_matrices_residue_combined_ft_incremental_avg_no_helper(
+        self,
+    ):
+        """
+        Test that: highest=True and combined_forcetorque=True
+        - initializes forcetorque_avg on first call
+        - updates it via incremental mean on second call
+        - avoids asserting the mutable 'existing' arg passed into the mock
+        """
+        universe_operations = UniverseOperations()
+        level_manager = LevelManager(universe_operations)
+
+        mol = MagicMock()
+        mol.trajectory.__getitem__.return_value = None
+
+        ft1 = np.array([[1.0, 1.0], [1.0, 1.0]], dtype=float)
+        ft2 = np.array([[3.0, 3.0], [3.0, 3.0]], dtype=float)
+
+        level_manager.get_combined_forcetorque_matrices = MagicMock(
+            side_effect=[ft1, ft2]
+        )
+        level_manager.get_matrices = MagicMock()
+
+        force_avg = {"ua": {}, "res": [None], "poly": [None]}
+        torque_avg = {"ua": {}, "res": [None], "poly": [None]}
+        forcetorque_avg = {"ua": {}, "res": [None], "poly": [None]}
+        frame_counts = {"ua": {}, "res": [None], "poly": [None]}
+
+        level_manager.update_force_torque_matrices(
+            entropy_manager=MagicMock(),
+            mol=mol,
+            group_id=0,
+            level="residue",
+            level_list=["united_atom", "residue"],
+            time_index=0,
+            num_frames=10,
+            force_avg=force_avg,
+            torque_avg=torque_avg,
+            forcetorque_avg=forcetorque_avg,
+            frame_counts=frame_counts,
+            force_partitioning=0.5,
+            combined_forcetorque=True,
+            customised_axes=True,
+        )
+
+        np.testing.assert_array_equal(forcetorque_avg["res"][0], ft1)
+        assert frame_counts["res"][0] == 1
+
+        level_manager.update_force_torque_matrices(
+            entropy_manager=MagicMock(),
+            mol=mol,
+            group_id=0,
+            level="residue",
+            level_list=["united_atom", "residue"],
+            time_index=1,
+            num_frames=10,
+            force_avg=force_avg,
+            torque_avg=torque_avg,
+            forcetorque_avg=forcetorque_avg,
+            frame_counts=frame_counts,
+            force_partitioning=0.5,
+            combined_forcetorque=True,
+            customised_axes=True,
+        )
+
+        expected = ft1 + (ft2 - ft1) / 2.0
+        np.testing.assert_array_almost_equal(forcetorque_avg["res"][0], expected)
+        assert frame_counts["res"][0] == 2
+
+        level_manager.get_matrices.assert_not_called()
+        assert level_manager.get_combined_forcetorque_matrices.call_count == 2
 
     def test_filter_zero_rows_columns_no_zeros(self):
         """
diff --git a/tests/test_CodeEntropy/test_mda_universe_operations.py b/tests/test_CodeEntropy/test_mda_universe_operations.py
index fb4c891..46e68a7 100644
--- a/tests/test_CodeEntropy/test_mda_universe_operations.py
+++ b/tests/test_CodeEntropy/test_mda_universe_operations.py
@@ -28,13 +28,6 @@ def setUp(self):
     def test_new_U_select_frame(self, MockMerge, MockAnalysisFromFunction):
         """
         Unit test for UniverseOperations.new_U_select_frame().
-
-        Verifies that:
-        - The Universe is queried with select_atoms("all", updating=True)
-        - AnalysisFromFunction is called to obtain coordinates and forces
-        - mda.Merge is called with the selected AtomGroup
-        - The new universe returned by Merge.load_new receives the correct arrays
-        - The method returns the merged universe
         """
         # Mock Universe and its components
         mock_universe = MagicMock()
@@ -48,6 +41,7 @@ def test_new_U_select_frame(self, MockMerge, MockAnalysisFromFunction):
         # Mock AnalysisFromFunction results for coordinates, forces, and dimensions
         coords = np.random.rand(10, 100, 3)
         forces = np.random.rand(10, 100, 3)
+        dims = np.random.rand(10, 6)
 
         mock_coords_analysis = MagicMock()
         mock_coords_analysis.run.return_value.results = {"timeseries": coords}
@@ -55,35 +49,35 @@ def test_new_U_select_frame(self, MockMerge, MockAnalysisFromFunction):
         mock_forces_analysis = MagicMock()
         mock_forces_analysis.run.return_value.results = {"timeseries": forces}
 
-        # Set the side effects for the three AnalysisFromFunction calls
+        mock_dims_analysis = MagicMock()
+        mock_dims_analysis.run.return_value.results = {"timeseries": dims}
+
         MockAnalysisFromFunction.side_effect = [
             mock_coords_analysis,
             mock_forces_analysis,
+            mock_dims_analysis,
         ]
 
-        # Mock the merge operation
+        # Mock merge operation
         mock_merged_universe = MagicMock()
         MockMerge.return_value = mock_merged_universe
 
         ops = UniverseOperations()
         result = ops.new_U_select_frame(mock_universe)
 
+        # Basic behavior checks
         mock_universe.select_atoms.assert_called_once_with("all", updating=True)
-        MockMerge.assert_called_once_with(mock_select_atoms)
 
-        # Ensure the 'load_new' method was called with the correct arguments
-        mock_merged_universe.load_new.assert_called_once()
-        args, kwargs = mock_merged_universe.load_new.call_args
+        # AnalysisFromFunction called 3 times (coords, forces, dimensions)
+        assert MockAnalysisFromFunction.call_count == 3
+        mock_coords_analysis.run.assert_called_once()
+        mock_forces_analysis.run.assert_called_once()
+        mock_dims_analysis.run.assert_called_once()
 
-        # Assert that the arrays are passed correctly
-        np.testing.assert_array_equal(args[0], coords)
-        np.testing.assert_array_equal(kwargs["forces"], forces)
-
-        # Check if format was included in the kwargs
-        self.assertIn("format", kwargs)
+        # Merge called with selected AtomGroup
+        MockMerge.assert_called_once_with(mock_select_atoms)
 
-        # Ensure the result is the mock merged universe
-        self.assertEqual(result, mock_merged_universe)
+        assert result == mock_merged_universe
 
     @patch("CodeEntropy.mda_universe_operations.AnalysisFromFunction")
     @patch("CodeEntropy.mda_universe_operations.mda.Merge")
@@ -93,7 +87,7 @@ def test_new_U_select_atom(self, MockMerge, MockAnalysisFromFunction):
 
         Ensures that:
         - The Universe is queried with the correct selection string
-        - Coordinates and forces are extracted via AnalysisFromFunction
+        - Coordinates, forces, and dimensions are extracted via AnalysisFromFunction
         - mda.Merge receives the AtomGroup from select_atoms
         - The new universe is populated with the expected data via load_new()
         - The returned universe is the object created by Merge
@@ -106,6 +100,7 @@ def test_new_U_select_atom(self, MockMerge, MockAnalysisFromFunction):
         # Mock AnalysisFromFunction results for coordinates, forces, and dimensions
         coords = np.random.rand(10, 100, 3)
         forces = np.random.rand(10, 100, 3)
+        dims = np.random.rand(10, 6)
 
         mock_coords_analysis = MagicMock()
         mock_coords_analysis.run.return_value.results = {"timeseries": coords}
@@ -113,10 +108,13 @@ def test_new_U_select_atom(self, MockMerge, MockAnalysisFromFunction):
         mock_forces_analysis = MagicMock()
         mock_forces_analysis.run.return_value.results = {"timeseries": forces}
 
-        # Set the side effects for the three AnalysisFromFunction calls
+        mock_dims_analysis = MagicMock()
+        mock_dims_analysis.run.return_value.results = {"timeseries": dims}
+
         MockAnalysisFromFunction.side_effect = [
             mock_coords_analysis,
             mock_forces_analysis,
+            mock_dims_analysis,
         ]
 
         # Mock the merge operation
@@ -126,23 +124,19 @@ def test_new_U_select_atom(self, MockMerge, MockAnalysisFromFunction):
         ops = UniverseOperations()
 
         result = ops.new_U_select_atom(mock_universe, select_string="resid 1-10")
-
         mock_universe.select_atoms.assert_called_once_with("resid 1-10", updating=True)
-        MockMerge.assert_called_once_with(mock_select_atoms)
-
-        # Ensure the 'load_new' method was called with the correct arguments
-        mock_merged_universe.load_new.assert_called_once()
-        args, kwargs = mock_merged_universe.load_new.call_args
 
-        # Assert that the arrays are passed correctly
-        np.testing.assert_array_equal(args[0], coords)
-        np.testing.assert_array_equal(kwargs["forces"], forces)
+        # AnalysisFromFunction called for coords, forces, dimensions
+        assert MockAnalysisFromFunction.call_count == 3
+        mock_coords_analysis.run.assert_called_once()
+        mock_forces_analysis.run.assert_called_once()
+        mock_dims_analysis.run.assert_called_once()
 
-        # Check if format was included in the kwargs
-        self.assertIn("format", kwargs)
+        # Merge called with the selected AtomGroup
+        MockMerge.assert_called_once_with(mock_select_atoms)
 
-        # Ensure the result is the mock merged universe
-        self.assertEqual(result, mock_merged_universe)
+        # Returned universe should be the merged universe
+        assert result == mock_merged_universe
 
     def test_get_molecule_container(self):
         """
@@ -179,25 +173,13 @@ def test_get_molecule_container(self):
     def test_merge_forces(self, MockUniverse, MockMerge, MockAnalysisFromFunction):
         """
         Unit test for UniverseOperations.merge_forces().
-
-        This test ensures that:
-        - Two MDAnalysis Universes are created: one for coordinates
-        (tprfile + trrfile) and one for forces (tprfile + forcefile).
-        - Both Universes correctly return AtomGroups via select_atoms("all").
-        - Coordinates and forces are extracted using AnalysisFromFunction.
-        - mda.Merge is called with the coordinate AtomGroup.
-        - The merged Universe receives the correct coordinate and force arrays
-        through load_new().
-        - When kcal=False, force values are passed through unchanged
-        (no kcal→kJ conversion).
-        - The returned universe is the same object returned by mda.Merge().
         """
-
+        # Two Universes created: coords and forces
         mock_u_coords = MagicMock()
         mock_u_force = MagicMock()
         MockUniverse.side_effect = [mock_u_coords, mock_u_force]
 
-        # Each universe returns a mock AtomGroup from select_atoms("all")
+        # Each universe returns an AtomGroup from select_atoms("all")
         mock_ag_coords = MagicMock()
         mock_ag_force = MagicMock()
         mock_u_coords.select_atoms.return_value = mock_ag_coords
@@ -205,6 +187,7 @@ def test_merge_forces(self, MockUniverse, MockMerge, MockAnalysisFromFunction):
 
         coords = np.random.rand(5, 10, 3)
         forces = np.random.rand(5, 10, 3)
+        dims = np.random.rand(5, 6)
 
         mock_coords_analysis = MagicMock()
         mock_coords_analysis.run.return_value.results = {"timeseries": coords}
@@ -212,10 +195,13 @@ def test_merge_forces(self, MockUniverse, MockMerge, MockAnalysisFromFunction):
         mock_forces_analysis = MagicMock()
         mock_forces_analysis.run.return_value.results = {"timeseries": forces}
 
-        # Two calls: first for coordinates, second for forces
+        mock_dims_analysis = MagicMock()
+        mock_dims_analysis.run.return_value.results = {"timeseries": dims}
+
         MockAnalysisFromFunction.side_effect = [
             mock_coords_analysis,
             mock_forces_analysis,
+            mock_dims_analysis,
         ]
 
         mock_merged = MagicMock()
@@ -230,28 +216,24 @@ def test_merge_forces(self, MockUniverse, MockMerge, MockAnalysisFromFunction):
             kcal=False,
         )
 
-        self.assertEqual(MockUniverse.call_count, 2)
-        MockUniverse.assert_any_call("topol.tpr", "traj.trr", format=None)
-        MockUniverse.assert_any_call("topol.tpr", "forces.trr", format=None)
+        # Universe construction
+        assert MockUniverse.call_count == 2
 
+        # Selection
         mock_u_coords.select_atoms.assert_called_once_with("all")
         mock_u_force.select_atoms.assert_called_once_with("all")
 
-        self.assertEqual(MockAnalysisFromFunction.call_count, 2)
+        # AnalysisFromFunction usage
+        assert MockAnalysisFromFunction.call_count == 3
+        mock_coords_analysis.run.assert_called_once()
+        mock_forces_analysis.run.assert_called_once()
+        mock_dims_analysis.run.assert_called_once()
 
+        # Merge called with coordinate AtomGroup
         MockMerge.assert_called_once_with(mock_ag_coords)
 
-        mock_merged.load_new.assert_called_once()
-        args, kwargs = mock_merged.load_new.call_args
-
-        # Coordinates passed positionally
-        np.testing.assert_array_equal(args[0], coords)
-
-        # Forces passed via kwargs
-        np.testing.assert_array_equal(kwargs["forces"], forces)
-
-        # Finally the function returns the merged universe
-        self.assertEqual(result, mock_merged)
+        # Returned object is merged universe
+        assert result == mock_merged
 
     @patch("CodeEntropy.mda_universe_operations.AnalysisFromFunction")
     @patch("CodeEntropy.mda_universe_operations.mda.Merge")
@@ -262,15 +244,6 @@ def test_merge_forces_kcal_conversion(
         """
         Unit test for UniverseOperations.merge_forces() covering the kcal→kJ
         conversion branch.
-
-        Verifies that:
-        - Two Universe objects are constructed for coords and forces.
-        - Each Universe returns an AtomGroup via select_atoms("all").
-        - AnalysisFromFunction is called twice.
-        - Forces are multiplied EXACTLY once by 4.184 when kcal=True.
-        - Merge() is called with the coordinate AtomGroup.
-        - load_new() receives the correct coordinates and converted forces.
-        - The returned Universe is the Merge() result.
         """
         mock_u_coords = MagicMock()
         mock_u_force = MagicMock()
@@ -286,6 +259,8 @@ def test_merge_forces_kcal_conversion(
         original_forces = np.ones((2, 3, 3))
         mock_forces_array = original_forces.copy()
 
+        dims = np.ones((2, 6))
+
         # Mock AnalysisFromFunction return values
         mock_coords_analysis = MagicMock()
         mock_coords_analysis.run.return_value.results = {"timeseries": coords}
@@ -295,9 +270,13 @@ def test_merge_forces_kcal_conversion(
             "timeseries": mock_forces_array
         }
 
+        mock_dims_analysis = MagicMock()
+        mock_dims_analysis.run.return_value.results = {"timeseries": dims}
+
         MockAnalysisFromFunction.side_effect = [
             mock_coords_analysis,
             mock_forces_analysis,
+            mock_dims_analysis,
         ]
 
         mock_merged = MagicMock()
@@ -306,10 +285,20 @@ def test_merge_forces_kcal_conversion(
         ops = UniverseOperations()
         result = ops.merge_forces("t.tpr", "c.trr", "f.trr", kcal=True)
 
-        _, kwargs = mock_merged.load_new.call_args
+        # select_atoms("all") (your code uses no updating=True)
+        mock_u_coords.select_atoms.assert_called_once_with("all")
+        mock_u_force.select_atoms.assert_called_once_with("all")
+
+        # AnalysisFromFunction called three times
+        assert MockAnalysisFromFunction.call_count == 3
 
-        expected_forces = original_forces * 4.184
-        np.testing.assert_array_equal(kwargs["forces"], expected_forces)
-        np.testing.assert_array_equal(mock_merged.load_new.call_args[0][0], coords)
+        # Forces are multiplied exactly once by 4.184 when kcal=True
+        np.testing.assert_allclose(
+            mock_forces_array, original_forces * 4.184, rtol=0, atol=0
+        )
+
+        # Merge called with coordinate AtomGroup
+        MockMerge.assert_called_once_with(mock_ag_coords)
 
-        self.assertEqual(result, mock_merged)
+        # Returned universe is the merged universe
+        assert result == mock_merged