Merge branch 'main' of https://github.com/CCPBioSim/CodeEntropy into 74-yaml-input-bug

Author: harryswift01

CodeEntropy/ConformationFunctions.py …py/calculations/ConformationFunctions.pyCodeEntropy/ConformationFunctions.py renamed to CodeEntropy/calculations/ConformationFunctions.py CodeEntropy/ConformationFunctions.py renamed to CodeEntropy/calculations/ConformationFunctions.py

@@ -1,5 +1,9 @@
+import logging
+
 import numpy as np
 
+logger = logging.getLogger(__name__)
+
 # from MDAnalysis.analysis.dihedrals import Dihedral
 
 
@@ -79,4 +83,6 @@ def assign_conformation(
         distances = [abs(phi[frame] - peak) for peak in peak_values]
         conformations[frame] = np.argmin(distances)
 
+    logger.debug(f"Final conformations: {conformations}")
+
     return conformations

CodeEntropy/EntropyFunctions.py …Entropy/calculations/EntropyFunctions.pyCodeEntropy/EntropyFunctions.py renamed to CodeEntropy/calculations/EntropyFunctions.py CodeEntropy/EntropyFunctions.py renamed to CodeEntropy/calculations/EntropyFunctions.py

@@ -1,3 +1,4 @@
+import logging
 import math
 
 # import matplotlib.pyplot as plt
@@ -7,8 +8,10 @@
 # import pandas as pd
 from numpy import linalg as la
 
-from CodeEntropy import ConformationFunctions as CONF
-from CodeEntropy import UnitsAndConversions as UAC
+from CodeEntropy.calculations import ConformationFunctions as CONF
+from CodeEntropy.calculations import UnitsAndConversions as UAC
+
+logger = logging.getLogger(__name__)
 
 # import sys
 # from ast import arg
@@ -40,15 +43,19 @@ def frequency_calculation(lambdas, temp):
     pi = np.pi
     # get kT in Joules from given temperature
     kT = UAC.get_KT2J(temp)
+    logger.debug(f"Temperature: {temp}, kT: {kT}")
 
     lambdas = np.array(lambdas)  # Ensure input is a NumPy array
+    logger.debug(f"Eigenvalues (lambdas): {lambdas}")
 
     # Check for negatives and raise an error if any are found
     if np.any(lambdas < 0):
+        logger.error(f"Negative eigenvalues encountered: {lambdas[lambdas < 0]}")
         raise ValueError(f"Negative eigenvalues encountered: {lambdas[lambdas < 0]}")
 
     # Compute frequencies safely
     frequencies = 1 / (2 * pi) * np.sqrt(lambdas / kT)
+    logger.debug(f"Calculated frequencies: {frequencies}")
 
     return frequencies
 
@@ -74,22 +81,31 @@ def vibrational_entropy(matrix, matrix_type, temp, highest_level):
     # N beads at a level => 3N x 3N covariance matrix => 3N eigenvalues
     # Get eigenvalues of the given matrix and change units to SI units
     lambdas = la.eigvals(matrix)
+    logger.debug(f"Eigenvalues (lambdas) before unit change: {lambdas}")
     lambdas = UAC.change_lambda_units(lambdas)
+    logger.debug(f"Eigenvalues (lambdas) after unit change: {lambdas}")
 
     # Calculate frequencies from the eigenvalues
     frequencies = frequency_calculation(lambdas, temp)
+    logger.debug(f"Calculated frequencies: {frequencies}")
 
     # Sort frequencies lowest to highest
     frequencies = np.sort(frequencies)
+    logger.debug(f"Sorted frequencies: {frequencies}")
 
     kT = UAC.get_KT2J(temp)
+    logger.debug(f"Temperature: {temp}, kT: {kT}")
     exponent = UAC.PLANCK_CONST * frequencies / kT
+    logger.debug(f"Exponent values: {exponent}")
     power_positive = np.power(np.e, exponent)
     power_negative = np.power(np.e, -exponent)
+    logger.debug(f"Power positive values: {power_positive}")
+    logger.debug(f"Power negative values: {power_negative}")
     S_components = exponent / (power_positive - 1) - np.log(1 - power_negative)
     S_components = (
         S_components * UAC.GAS_CONST
     )  # multiply by R - get entropy in J mol^{-1} K^{-1}
+    logger.debug(f"Entropy components: {S_components}")
     # N beads at a level => 3N x 3N covariance matrix => 3N eigenvalues
     if matrix_type == "force":  # force covariance matrix
         if highest_level:  # whole molecule level - we take all frequencies into account
@@ -104,6 +120,8 @@ def vibrational_entropy(matrix, matrix_type, temp, highest_level):
     else:  # torque covariance matrix - we always take all values into account
         S_vib_total = sum(S_components)
 
+    logger.debug(f"Total vibrational entropy: {S_vib_total}")
+
     return S_vib_total
 
 
@@ -138,17 +156,23 @@ def conformational_entropy(
         )
         index += 1
 
+    logger.debug(f"Conformation matrix: {conformation}")
+
     # For each frame, convert the conformation of all dihedrals into a state string
     states = ["" for x in range(number_frames)]
     for frame_index in range(number_frames):
         for index in range(num_dihedrals):
             states[frame_index] += str(conformation[index][frame_index])
 
+    logger.debug(f"States: {states}")
+
     # Count how many times each state occurs, then use the probability to get the
     # entropy
     # entropy = sum over states p*ln(p)
     values, counts = np.unique(states, return_counts=True)
     for state in range(len(values)):
+        logger.debug(f"Unique states: {values}")
+        logger.debug(f"Counts: {counts}")
         count = counts[state]
         probability = count / number_frames
         entropy = probability * np.log(probability)
@@ -157,6 +181,8 @@ def conformational_entropy(
     # multiply by gas constant to get the units J/mol/K
     S_conf_total *= -1 * UAC.GAS_CONST
 
+    logger.debug(f"Total conformational entropy: {S_conf_total}")
+
     return S_conf_total
 
 
@@ -193,12 +219,28 @@ def orientational_entropy(neighbours_dict):
         else:
             # the bound ligand is always going to be a neighbour
             omega = np.sqrt((neighbours_dict[neighbour] ** 3) * math.pi)
+            logger.debug(f"Omega for neighbour {neighbour}: {omega}")
             # orientational entropy arising from each neighbouring species
             # - we know the species is going to be a neighbour
             S_or_component = math.log(omega)
+            logger.debug(
+                f"S_or_component (log(omega)) for neighbour {neighbour}: "
+                f"{S_or_component}"
+            )
             S_or_component *= UAC.GAS_CONST
+            logger.debug(
+                f"S_or_component after multiplying by GAS_CONST for neighbour "
+                f"{neighbour}: {S_or_component}"
+            )
         S_or_total += S_or_component
+        logger.debug(
+            f"S_or_total after adding component for neighbour {neighbour}: "
+            f"{S_or_total}"
+        )
     # TODO for future releases
     # implement a case for molecules with hydrogen bonds but to a lesser
     # extent than water
+
+    logger.debug(f"Final total orientational entropy: {S_or_total}")
+
     return S_or_total

CodeEntropy/GeometricFunctions.py …tropy/calculations/GeometricFunctions.pyCodeEntropy/GeometricFunctions.py renamed to CodeEntropy/calculations/GeometricFunctions.py CodeEntropy/GeometricFunctions.py renamed to CodeEntropy/calculations/GeometricFunctions.py

@@ -1,5 +1,9 @@
+import logging
+
 import numpy as np
 
+logger = logging.getLogger(__name__)
+
 
 def get_beads(data_container, level):
     """
@@ -40,6 +44,8 @@ def get_beads(data_container, level):
             )
             list_of_beads.append(data_container.select_atoms(atom_group))
 
+    logger.debug(f"List of beads: {list_of_beads}")
+
     return list_of_beads
 
 
@@ -120,6 +126,9 @@ def get_axes(data_container, level, index=0):
         # use spherical coordinates function to get rotational axes
         rot_axes = get_sphCoord_axes(vector)
 
+    logger.debug(f"Translational Axes: {trans_axes}")
+    logger.debug(f"Rotational Axes: {rot_axes}")
+
     return trans_axes, rot_axes
 
 
@@ -159,6 +168,8 @@ def get_avg_pos(atom_set, center):
     # transform the average position to a coordinate system with the origin at center
     avg_position = avg_position - center
 
+    logger.debug(f"Average Position: {avg_position}")
+
     return avg_position
 
 
@@ -231,6 +242,8 @@ def get_sphCoord_axes(arg_r):
     # Phi^
     spherical_basis[2, :] = np.asarray([-sin_phi, cos_phi, 0.0])
 
+    logger.debug(f"Spherical Basis: {spherical_basis}")
+
     return spherical_basis
 
 
@@ -276,6 +289,8 @@ def get_weighted_forces(
 
     weighted_force = forces_trans / np.sqrt(mass)
 
+    logger.debug(f"Weighted Force: {weighted_force}")
+
     return weighted_force
 
 
@@ -361,6 +376,8 @@ def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5)
             moment_of_inertia[dimension, dimension]
         )
 
+    logger.debug(f"Weighted Torque: {weighted_torque}")
+
     return weighted_torque
 
 
@@ -390,6 +407,8 @@ def create_submatrix(data_i, data_j, number_frames):
     # Divide by the number of frames to get the average
     submatrix /= number_frames
 
+    logger.debug(f"Submatrix: {submatrix}")
+
     return submatrix
 
 
@@ -432,11 +451,13 @@ def filter_zero_rows_columns(arg_matrix, verbose):
     # get the final shape
     final_shape = np.shape(arg_matrix)
 
-    if verbose and init_shape != final_shape:
-        print(
-            "A shape change has occured ({},{}) -> ({}, {})".format(
+    if init_shape != final_shape:
+        logger.debug(
+            "A shape change has occurred ({},{}) -> ({}, {})".format(
                 *init_shape, *final_shape
             )
         )
 
+    logger.debug(f"arg_matrix: {arg_matrix}")
+
     return arg_matrix

CodeEntropy/LevelFunctions.py …deEntropy/calculations/LevelFunctions.pyCodeEntropy/LevelFunctions.py renamed to CodeEntropy/calculations/LevelFunctions.py CodeEntropy/LevelFunctions.py renamed to CodeEntropy/calculations/LevelFunctions.py

@@ -1,7 +1,11 @@
 # import MDAnalysis as mda
+import logging
+
 import numpy as np
 
-from CodeEntropy import GeometricFunctions as GF
+from CodeEntropy.calculations import GeometricFunctions as GF
+
+logger = logging.getLogger(__name__)
 
 
 def select_levels(data_container, verbose):
@@ -23,7 +27,7 @@ def select_levels(data_container, verbose):
 
     # fragments is MDAnalysis terminology for what chemists would call molecules
     number_molecules = len(data_container.atoms.fragments)
-    print("The number of molecules is {}.".format(number_molecules))
+    logger.debug("The number of molecules is {}.".format(number_molecules))
     fragments = data_container.atoms.fragments
     levels = [[] for _ in range(number_molecules)]
 
@@ -42,8 +46,7 @@ def select_levels(data_container, verbose):
             if number_residues > 1:
                 levels[molecule].append("polymer")
 
-    if verbose:
-        print(levels)
+    logger.debug(f"levels {levels}")
 
     return number_molecules, levels
 
@@ -142,12 +145,8 @@ def get_matrices(
     force_matrix = GF.filter_zero_rows_columns(force_matrix, verbose)
     torque_matrix = GF.filter_zero_rows_columns(torque_matrix, verbose)
 
-    if verbose:
-        with open("matrix.out", "a") as f:
-            print("force_matrix \n", file=f)
-            print(force_matrix, file=f)
-            print("torque_matrix \n", file=f)
-            print(torque_matrix, file=f)
+    logger.debug(f"Force Matrix: {force_matrix}")
+    logger.debug(f"Torque Matrix: {torque_matrix}")
 
     return force_matrix, torque_matrix
 
@@ -181,7 +180,7 @@ def get_dihedrals(data_container, level):
     if level == "residue":
         num_residues = len(data_container.residues)
         if num_residues < 4:
-            print("no residue level dihedrals")
+            logger.debug("no residue level dihedrals")
 
         else:
             # find bonds between residues N-3:N-2 and N-1:N
@@ -224,4 +223,6 @@ def get_dihedrals(data_container, level):
                 atom_group = atom1 + atom2 + atom3 + atom4
                 dihedrals.append(atom_group.dihedral)
 
+    logger.debug(f"Dihedrals: {dihedrals}")
+
     return dihedrals

CodeEntropy/MDAUniverseHelper.py …ntropy/calculations/MDAUniverseHelper.pyCodeEntropy/MDAUniverseHelper.py renamed to CodeEntropy/calculations/MDAUniverseHelper.py CodeEntropy/MDAUniverseHelper.py renamed to CodeEntropy/calculations/MDAUniverseHelper.py

@@ -1,9 +1,12 @@
+import logging
 import pickle
 
 import MDAnalysis as mda
 from MDAnalysis.analysis.base import AnalysisFromFunction
 from MDAnalysis.coordinates.memory import MemoryReader
 
+logger = logging.getLogger(__name__)
+
 
 def new_U_select_frame(u, start=None, end=None, step=1):
     """Create a reduced universe by dropping frames according to user selection
@@ -46,6 +49,7 @@ def new_U_select_frame(u, start=None, end=None, step=1):
     )
     u2 = mda.Merge(select_atom)
     u2.load_new(coordinates, format=MemoryReader, forces=forces, dimensions=dimensions)
+    logger.debug(f"MDAnalysis.Universe - reduced universe: {u2}")
     return u2
 
 
@@ -83,6 +87,7 @@ def new_U_select_atom(u, select_string="all"):
     )
     u2 = mda.Merge(select_atom)
     u2.load_new(coordinates, format=MemoryReader, forces=forces, dimensions=dimensions)
+    logger.debug(f"MDAnalysis.Universe - reduced universe: {u2}")
     return u2
 
 

CodeEntropy/NeighbourFunctions.py …tropy/calculations/NeighbourFunctions.pyCodeEntropy/NeighbourFunctions.py renamed to CodeEntropy/calculations/NeighbourFunctions.py CodeEntropy/NeighbourFunctions.py renamed to CodeEntropy/calculations/NeighbourFunctions.py

@@ -1,10 +1,14 @@
 # import os
+import logging
 import sys
 
 # import matplotlib.pyplot as plt
 import MDAnalysis as mda
 import numpy as np
 
+logger = logging.getLogger(__name__)
+
+
 # from ast import arg
 
 
@@ -85,4 +89,8 @@ def get_neighbours(molecule_i, reduced_atom):
                 # print(r_ij)
                 neighbours_array.append(molecule_j.atoms.resids[0])
                 neighbours_dict[molecule_j.atoms.resnames[0]] = 1
+
+    logger.debug(f"Neighbours dictionary: {neighbours_dict}")
+    logger.debug(f"Neighbours array: {neighbours_array}")
+
     return neighbours_dict, neighbours_array