Update Documentation for v1.0 Release Candidate of CodeEntropy

.gitignore

@@ -64,6 +64,7 @@ instance/
 
 # Sphinx documentation
 docs/_build/
+docs/autosummary/
 
 # PyBuilder
 target/

CodeEntropy/entropy.py

@@ -38,6 +38,8 @@ def __init__(
             universe: MDAnalysis universe representing the simulation system.
             data_logger: Logger for storing and exporting entropy data.
             level_manager: Provides level-specific data such as matrices and dihedrals.
+            group_molecules: includes the grouping functions for averaging over
+            molecules.
         """
         self._run_manager = run_manager
         self._args = args
@@ -58,6 +60,7 @@ def execute(self):
         - Computing vibrational and conformational entropies.
         - Finalizing and logging results.
         """
+        # Set up initial information
         start, end, step = self._get_trajectory_bounds()
         number_frames = self._get_number_frames(start, end, step)
 
@@ -118,6 +121,8 @@ def execute(self):
             )
         )
 
+        # Identify the conformational states from dihedral angles for the
+        # conformational entropy calculations
         states_ua, states_res = self._level_manager.build_conformational_states(
             self,
             reduced_atom,
@@ -131,6 +136,7 @@ def execute(self):
             ce,
         )
 
+        # Complete the entropy calculations
         self._compute_entropies(
             reduced_atom,
             levels,
@@ -145,6 +151,7 @@ def execute(self):
             ce,
         )
 
+        # Print the results in a nicely formated way
         self._finalize_molecule_results()
         self._data_logger.log_tables()
 
@@ -190,9 +197,15 @@ def _initialize_molecules(self):
                 - reduced_atom (Universe): The reduced atom selection.
                 - number_molecules (int): Number of molecules in the system.
                 - levels (list): List of entropy levels per molecule.
+                - groups (dict): Groups for averaging over molecules.
         """
+        # Based on the selection string, create a new MDAnalysis universe
         reduced_atom = self._get_reduced_universe()
+
+        # Count the molecules and identify the length scale levels for each one
         number_molecules, levels = self._level_manager.select_levels(reduced_atom)
+
+        # Group the molecules for averaging
         grouping = self._args.grouping
         groups = self._group_molecules.grouping_molecules(reduced_atom, grouping)
 
@@ -226,14 +239,16 @@ def _compute_entropies(
 
         Parameters:
             reduced_atom (Universe): The reduced atom selection from the trajectory.
-            number_molecules (int): Number of molecules in the system.
             levels (list): List of entropy levels per molecule.
+            groups (dict): Groups for averaging over molecules.
             force_matrices (dict): Precomputed force covariance matrices.
             torque_matrices (dict): Precomputed torque covariance matrices.
             states_ua (dict): Dictionary to store united-atom conformational states.
             states_res (list): List to store residue-level conformational states.
             frames_count (dict): Dictionary to store the frame counts
             number_frames (int): Total number of trajectory frames to process.
+            ve: Vibrational Entropy object
+            ce: Conformational Entropy object
         """
         with Progress(
             SpinnerColumn(),
@@ -363,8 +378,11 @@ def _get_reduced_universe(self):
         Returns:
             MDAnalysis.Universe: Selected subset of the system.
         """
+        # If selection string is "all" the universe does not change
         if self._args.selection_string == "all":
             return self._universe
+
+        # Otherwise create a new (smaller) universe based on the selection
         reduced = self._run_manager.new_U_select_atom(
             self._universe, self._args.selection_string
         )
@@ -383,8 +401,11 @@ def _get_molecule_container(self, universe, molecule_id):
         Returns:
             MDAnalysis.Universe: Universe containing only the selected molecule.
         """
+        # Identify the atoms in the molecule
         frag = universe.atoms.fragments[molecule_id]
         selection_string = f"index {frag.indices[0]}:{frag.indices[-1]}"
+
+        # Build a new universe with only the one molecule
         return self._run_manager.new_U_select_atom(universe, selection_string)
 
     def _process_united_atom_entropy(
@@ -406,7 +427,7 @@ def _process_united_atom_entropy(
         united-atom level.
 
         Args:
-            mol_id (int): ID of the molecule.
+            group_id (int): ID of the group.
             mol_container (Universe): Universe for the selected molecule.
             ve: VibrationalEntropy object.
             ce: ConformationalEntropy object.
@@ -415,43 +436,56 @@ def _process_united_atom_entropy(
             n_frames (int): Number of trajectory frames.
             frame_counts: Number of frames counted
             highest (bool): Whether this is the highest level of resolution for
-            the molecule.
+             the molecule.
+            number_frames (int): The number of frames analysed.
         """
         S_trans, S_rot, S_conf = 0, 0, 0
 
+        # The united atom entropy is calculated separately for each residue
+        # This is to allow residue by residue information
+        # and prevents the matrices from becoming too large
         for residue_id, residue in enumerate(mol_container.residues):
 
             key = (group_id, residue_id)
 
+            # Find the relevant force and torque matrices and tidy them up
+            # by removing rows and columns that are all zeros
             f_matrix = force_matrix[key]
             f_matrix = self._level_manager.filter_zero_rows_columns(f_matrix)
 
             t_matrix = torque_matrix[key]
             t_matrix = self._level_manager.filter_zero_rows_columns(t_matrix)
 
+            # Calculate the vibrational entropy
             S_trans_res = ve.vibrational_entropy_calculation(
                 f_matrix, "force", self._args.temperature, highest
             )
             S_rot_res = ve.vibrational_entropy_calculation(
                 t_matrix, "torque", self._args.temperature, highest
             )
 
+            # Get the relevant conformational states
             values = states[key]
-
+            # Check if there is information in the states array
             contains_non_empty_states = (
                 np.any(values) if isinstance(values, np.ndarray) else any(values)
             )
 
+            # Calculate the conformational entropy
+            # If there are no conformational states (i.e. no dihedrals)
+            # then the conformational entropy is zero
             S_conf_res = (
                 ce.conformational_entropy_calculation(values, number_frames)
                 if contains_non_empty_states
                 else 0
             )
 
+            # Add the data to the united atom level entropy
             S_trans += S_trans_res
             S_rot += S_rot_res
             S_conf += S_conf_res
 
+            # Print out the data for each residue
             self._data_logger.add_residue_data(
                 group_id,
                 residue.resname,
@@ -477,6 +511,7 @@ def _process_united_atom_entropy(
                 S_conf_res,
             )
 
+        # Print the total united atom level data 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)
         self._data_logger.add_results_data(group_id, level, "Conformational", S_conf)
@@ -502,8 +537,7 @@ def _process_vibrational_entropy(
         Calculates vibrational entropy.
 
         Args:
-            mol_id (int): Molecule ID.
-            mol_container (Universe): Selected molecule's universe.
+            group_id (int): Group ID.
             ve: VibrationalEntropy object.
             level (str): Current granularity level.
             force_matrix : Force covariance matrix
@@ -512,17 +546,21 @@ def _process_vibrational_entropy(
             highest (bool): Flag indicating if this is the highest granularity
             level.
         """
+        # 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)
 
         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)
 
@@ -547,9 +585,11 @@ def _process_conformational_entropy(
             mol_container (Universe): Selected molecule's universe.
             ce: ConformationalEntropy object.
             level (str): Level name (should be 'residue').
-            start, end, step (int): Frame bounds.
-            n_frames (int): Number of frames used.
+            states (array): The conformational states.
+            number_frames (int): Number of frames used.
         """
+        # Get the relevant conformational states
+        # Check if there is information in the states array
         group_states = states[group_id] if group_id < len(states) else None
 
         if group_states is not None:
@@ -561,6 +601,9 @@ def _process_conformational_entropy(
         else:
             contains_state_data = False
 
+        # Calculate the conformational entropy
+        # If there are no conformational states (i.e. no dihedrals)
+        # then the conformational entropy is zero
         S_conf = (
             ce.conformational_entropy_calculation(group_states, number_frames)
             if contains_state_data
@@ -784,14 +827,12 @@ def frequency_calculation(self, lambdas, temp):
 
         frequency=sqrt(λ/kT)/2π
 
-        Input
-        -----
-        lambdas : array of floats - eigenvalues of the covariance matrix
-        temp: float - temperature
+        Args:
+            lambdas : array of floats - eigenvalues of the covariance matrix
+            temp: float - temperature
 
-        Returns
-        -------
-        frequencies : array of floats - corresponding vibrational frequencies
+        Returns:
+            frequencies : array of floats - corresponding vibrational frequencies
         """
         pi = np.pi
         # get kT in Joules from given temperature
@@ -801,11 +842,16 @@ def frequency_calculation(self, lambdas, temp):
         lambdas = np.array(lambdas)  # Ensure input is a NumPy array
         logger.debug(f"Eigenvalues (lambdas): {lambdas}")
 
+        # Filter out lambda values that are negative or imaginary numbers
+        # As these will produce supurious entropy results that can crash
+        # the calculation
         lambdas = np.real_if_close(lambdas, tol=1000)
         valid_mask = (
             np.isreal(lambdas) & (lambdas > 0) & (~np.isclose(lambdas, 0, atol=1e-07))
         )
 
+        # If any lambdas were removed by the filter, warn the user
+        # as this will suggest insufficient sampling in the simulation data
         if len(lambdas) > np.count_nonzero(valid_mask):
             logger.warning(
                 f"{len(lambdas) - np.count_nonzero(valid_mask)} "
@@ -827,16 +873,14 @@ def vibrational_entropy_calculation(self, matrix, matrix_type, temp, highest_lev
         Physics, 2018, 116, 1965–1976 / eq. (2) in A. Chakravorty, J. Higham and
         R. H. Henchman, J. Chem. Inf. Model., 2020, 60, 5540–5551.
 
-        Input
-        -----
-        matrix : matrix - force/torque covariance matrix
-        matrix_type: string
-        temp: float - temperature
-        highest_level: bool - is this the highest level of the heirarchy
+        Args:
+            matrix : matrix - force/torque covariance matrix
+            matrix_type: string
+            temp: float - temperature
+            highest_level: bool - is this the highest level of the heirarchy
 
-        Returns
-        -------
-        S_vib_total : float - transvibrational/rovibrational entropy
+        Returns:
+            S_vib_total : float - transvibrational/rovibrational entropy
         """
         # N beads at a level => 3N x 3N covariance matrix => 3N eigenvalues
         # Get eigenvalues of the given matrix and change units to SI units
@@ -917,18 +961,18 @@ def assign_conformation(
         Based on the identified TPs, states are assigned to each configuration of the
         dihedral.
 
-        Input
-        -----
-        dihedral_atom_group : the group of 4 atoms defining the dihedral
-        number_frames : number of frames in the trajectory
-        bin_width : the width of the histogram bit, default 30 degrees
-        start : int, starting frame, will default to 0
-        end : int, ending frame, will default to -1 (last frame in trajectory)
-        step : int, spacing between frames, will default to 1
+        Args:
+            data_container (MDAnalysis Universe): data for the molecule/residue unit
+            dihedral (array): The dihedral angles in the unit
+            number_frames (int): number of frames in the trajectory
+            bin_width (int): the width of the histogram bit, default 30 degrees
+            start (int): starting frame, will default to 0
+            end (int): ending frame, will default to -1 (last frame in trajectory)
+            step (int): spacing between frames, will default to 1
 
-        Return
-        ------
-        A timeseries with integer labels describing the state at each point in time.
+        Returns:
+            conformations (array): A timeseries with integer labels describing the
+            state at each point in time.
 
         """
         conformations = np.zeros(number_frames)
@@ -943,7 +987,7 @@ def assign_conformation(
             timestep_index = timestep_index
             value = dihedral.value()
             # we want postive values in range 0 to 360 to make the peak assignment
-            # work using the fact that dihedrals have circular symetry
+            # works using the fact that dihedrals have circular symetry
             # (i.e. -15 degrees = +345 degrees)
             if value < 0:
                 value += 360
@@ -958,7 +1002,8 @@ def assign_conformation(
 
         # identify "convex turning-points" and populate a list of peaks
         # peak : a bin whose neighboring bins have smaller population
-        # NOTE might have problems if the peak is wide with a flat or sawtooth top
+        # NOTE might have problems if the peak is wide with a flat or sawtooth
+        # top in which case check you have a sensible bin width
         peak_values = []
 
         for bin_index in range(number_bins):
@@ -1001,12 +1046,12 @@ def conformational_entropy_calculation(self, states, number_frames):
 
         Uses the adaptive enumeration method (AEM).
 
-        Input
-        -----
-        dihedrals : array - array of dihedrals in the molecule
-        Returns
-        -------
-        S_conf_total : float - conformational entropy
+        Args:
+           states (array): Conformational states in the molecule
+           number_frames (int): The number of frames analysed
+
+        Returns:
+            S_conf_total (float) : conformational entropy
         """
 
         S_conf_total = 0