Refactor and Reformat entire codebase:

- Introduced a .pre-commi-config.yaml to enforce standards when commiting to this repository.
 - Includes Black, pre-commit-hooks, flake8 and isort.
- Applied 'pre-commit run --all-files', this automatically formatted a majority of the code, some manual intervention was needed.
- Addressed issue #41 with ensuring any instances of numpy is refered to as np. and not .nmp

Author: harryswift01

.pre-commit-config.yaml

@@ -15,20 +15,14 @@ repos:
       - id: check-toml
       - id: check-yaml
 
-  # - repo: https://github.com/pycqa/flake8
-  #   rev: 7.0.0
-  #   hooks:
-  #     - id: flake8
-  #       additional_dependencies: [flake8-pyproject]
+  - repo: https://github.com/pycqa/flake8
+    rev: 7.0.0
+    hooks:
+      - id: flake8
+        additional_dependencies: [flake8-pyproject]
 
   - repo: https://github.com/pycqa/isort
     rev: 5.12.0
     hooks:
       - id: isort
-        args: ["--profile=black"]
-
-  # - repo: https://github.com/PyCQA/pylint
-  #   rev: v3.3.4
-  #   hooks:
-  #     - id: pylint
-  #       additional_dependencies: ['pylint>=2.16.0']
+        args: ["--profile=black"]

CodeEntropy/ConformationFunctions.py

@@ -1,5 +1,6 @@
-import numpy as nmp
-from MDAnalysis.analysis.dihedrals import Dihedral
+import numpy as np
+
+# from MDAnalysis.analysis.dihedrals import Dihedral
 
 
 def assign_conformation(
@@ -10,7 +11,8 @@ def assign_conformation(
     as a function of time. The function creates a histogram from the timeseries of the
     dihedral angle values and identifies points of dominant occupancy
     (called CONVEX TURNING POINTS).
-    Based on the identified TPs, states are assigned to each configuration of the dihedral.
+    Based on the identified TPs, states are assigned to each configuration of the
+    dihedral.
 
     Input
     -----
@@ -26,22 +28,23 @@ def assign_conformation(
     A timeseries with integer labels describing the state at each point in time.
 
     """
-    conformations = nmp.zeros(number_frames)
-    phi = nmp.zeros(number_frames)
+    conformations = np.zeros(number_frames)
+    phi = np.zeros(number_frames)
 
     # get the values of the angle for the dihedral
     # dihedral angle values have a range from -180 to 180
     for timestep in data_container.trajectory[start:end:step]:
         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 (i.e. -15 degrees = +345 degrees)
+        # using the fact that dihedrals have circular symetry
+        # (i.e. -15 degrees = +345 degrees)
         if value < 0:
             value += 360
         phi[timestep.frame] = value
 
     # create a histogram using numpy
     number_bins = int(360 / bin_width)
-    popul, bin_edges = nmp.histogram(a=phi, bins=number_bins, range=(0, 360))
+    popul, bin_edges = np.histogram(a=phi, bins=number_bins, range=(0, 360))
     bin_value = [0.5 * (bin_edges[i] + bin_edges[i + 1]) for i in range(0, len(popul))]
 
     # identify "convex turning-points" and populate a list of peaks
@@ -53,7 +56,8 @@ def assign_conformation(
         # if there is no dihedrals in a bin then it cannot be a peak
         if popul[bin_index] == 0:
             pass
-        # being careful of the last bin (dihedrals have circular symmetry, the histogram does not)
+        # being careful of the last bin
+        # (dihedrals have circular symmetry, the histogram does not)
         elif (
             bin_index == number_bins - 1
         ):  # the -1 is because the index starts with 0 not 1
@@ -73,7 +77,7 @@ def assign_conformation(
     for frame in range(number_frames):
         # find the TP that the snapshot is least distant from
         distances = [abs(phi[frame] - peak) for peak in peak_values]
-        conformations[frame] = nmp.argmin(distances)
+        conformations[frame] = np.argmin(distances)
 
     return conformations
 

CodeEntropy/EntropyFunctions.py

@@ -1,27 +1,30 @@
 import math
-import sys
-from ast import arg
 
-import matplotlib.pyplot as plt
-import MDAnalysis as mda
-import numpy as nmp
-import pandas as pd
+# import matplotlib.pyplot as plt
+# import MDAnalysis as mda
+import numpy as np
+
+# import pandas as pd
 from numpy import linalg as la
 
 from CodeEntropy import ConformationFunctions as CONF
 from CodeEntropy import UnitsAndConversions as UAC
 
+# import sys
+# from ast import arg
+
+
 # from CodeEntropy import NeighbourFunctions as NF
 
 
 def frequency_calculation(lambdas, temp):
     """
-    Function to calculate an array of vibrational frequencies from the eigenvalues of the
-    covariance matrix.
+    Function to calculate an array of vibrational frequencies from the eigenvalues of
+    the covariance matrix.
 
-    Calculated from eq. (3) in Higham, S.-Y. Chou, F. Gräter and  R. H. Henchman, Molecular
-    Physics, 2018, 116, 1965–1976//eq. (3) in A. Chakravorty, J. Higham and R. H. Henchman,
-    J. Chem. Inf. Model., 2020, 60, 5540–5551
+    Calculated from eq. (3) in Higham, S.-Y. Chou, F. Gräter and  R. H. Henchman,
+    Molecular Physics, 2018, 116, 1965–1976//eq. (3) in A. Chakravorty, J. Higham
+    and R. H. Henchman, J. Chem. Inf. Model., 2020, 60, 5540–5551
 
     frequency=sqrt(λ/kT)/2π
 
@@ -34,10 +37,10 @@ def frequency_calculation(lambdas, temp):
     -------
        frequencies : array of floats - corresponding vibrational frequencies
     """
-    pi = nmp.pi
+    pi = np.pi
     # get kT in Joules from given temperature
     kT = UAC.get_KT2J(temp)
-    frequencies = 1 / (2 * pi) * nmp.sqrt(lambdas / kT)
+    frequencies = 1 / (2 * pi) * np.sqrt(lambdas / kT)
 
     return frequencies
 
@@ -47,9 +50,9 @@ def frequency_calculation(lambdas, temp):
 
 def vibrational_entropy(matrix, matrix_type, temp, highest_level):
     """
-    Function to calculate the vibrational entropy for each level calculated from eq. (4) in
-    J. Higham, S.-Y. Chou, F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116,
-    1965–1976 / eq. (2) in A. Chakravorty, J. Higham and R. H. Henchman,
+    Function to calculate the vibrational entropy for each level calculated from eq. (4)
+    in J. Higham, S.-Y. Chou, F. Gräter and R. H. Henchman, Molecular 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
@@ -72,13 +75,13 @@ def vibrational_entropy(matrix, matrix_type, temp, highest_level):
     frequencies = frequency_calculation(lambdas, temp)
 
     # Sort frequencies lowest to highest
-    frequencies = nmp.sort(frequencies)
+    frequencies = np.sort(frequencies)
 
     kT = UAC.get_KT2J(temp)
     exponent = UAC.PLANCK_CONST * frequencies / kT
-    power_positive = nmp.power(nmp.e, exponent)
-    power_negative = nmp.power(nmp.e, -exponent)
-    S_components = exponent / (power_positive - 1) - nmp.log(1 - power_negative)
+    power_positive = np.power(np.e, exponent)
+    power_negative = np.power(np.e, -exponent)
+    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}
@@ -87,8 +90,9 @@ def vibrational_entropy(matrix, matrix_type, temp, highest_level):
         if highest_level:  # whole molecule level - we take all frequencies into account
             S_vib_total = sum(S_components)
 
-        # discard the 6 lowest frequencies to discard translation and rotation of the whole unit
-        # the overall translation and rotation of a unit is an internal motion of the level above
+        # discard the 6 lowest frequencies to discard translation and rotation of the
+        # whole unit the overall translation and rotation of a unit is an internal
+        # motion of the level above
         else:
             S_vib_total = sum(S_components[6:])
 
@@ -107,7 +111,8 @@ def conformational_entropy(
     """
     Function to calculate conformational entropies using eq. (7) in Higham, S.-Y. Chou,
     F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116, 1965–1976 / eq. (4) in
-    A. Chakravorty, J. Higham and R. H. Henchman, J. Chem. Inf. Model., 2020, 60, 5540–5551.
+    A. Chakravorty, J. Higham and R. H. Henchman, J. Chem. Inf. Model., 2020, 60,
+    5540–5551.
 
     Uses the adaptive enumeration method (AEM).
 
@@ -123,7 +128,7 @@ def conformational_entropy(
 
     # For each dihedral, identify the conformation in each frame
     num_dihedrals = len(dihedrals)
-    conformation = nmp.zeros((num_dihedrals, number_frames))
+    conformation = np.zeros((num_dihedrals, number_frames))
     index = 0
     for dihedral in dihedrals:
         conformation[index] = CONF.assign_conformation(
@@ -137,13 +142,14 @@ def conformational_entropy(
         for index in range(num_dihedrals):
             states[frame_index] += str(conformation[index][frame_index])
 
-    # Count how many times each state occurs, then use the probability to get the entropy
+    # Count how many times each state occurs, then use the probability to get the
+    # entropy
     # entropy = sum over states p*ln(p)
-    values, counts = nmp.unique(states, return_counts=True)
+    values, counts = np.unique(states, return_counts=True)
     for state in range(len(values)):
         count = counts[state]
         probability = count / number_frames
-        entropy = probability * nmp.log(probability)
+        entropy = probability * np.log(probability)
         S_conf_total += entropy
 
     # multiply by gas constant to get the units J/mol/K
@@ -157,15 +163,17 @@ def conformational_entropy(
 
 def orientational_entropy(neighbours_dict):
     """
-    Function to calculate orientational entropies from eq. (10) in J. Higham, S.-Y. Chou,
-    F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116,3 1965–1976. Number of
-    orientations, Ω, is calculated using eq. (8) in  J. Higham, S.-Y. Chou, F. Gräter and
-    R. H. Henchman,  Molecular Physics, 2018, 116,3 1965–1976.
+    Function to calculate orientational entropies from eq. (10) in J. Higham, S.-Y.
+    Chou, F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116,3 1965–1976.
+    Number of orientations, Ω, is calculated using eq. (8) in  J. Higham,
+    S.-Y. Chou, F. Gräter and R. H. Henchman,  Molecular Physics, 2018, 116,
+    3 1965–1976.
 
     σ is assumed to be 1 for the molecules we're concerned with and hence,
     max {1, (Nc^3*π)^(1/2)} will always be (Nc^3*π)^(1/2).
 
-    TODO future release - function for determing symmetry and symmetry numbers maybe?
+    TODO future release - function for determing symmetry and symmetry numbers
+    maybe?
 
     Input
     -----
@@ -177,20 +185,23 @@ def orientational_entropy(neighbours_dict):
     -------
        S_or_total : float - orientational entropy
     """
+
+    # Replaced molecule with neighbour as this is what the for loop uses
     S_or_total = 0
     for neighbour in neighbours_dict:  # we are going through neighbours
-        if molecule in []:  # water molecules - call POSEIDON functions
+        if neighbour in []:  # water molecules - call POSEIDON functions
             pass  # TODO temporary until function is written
         else:
-            omega = nmp.sqrt(
-                (neighbours_dict[molecule] ** 3) * math.pi
-            )  # the bound ligand is always going to be a neighbour
-            S_or_component = math.log(
-                omega
-            )  # orientational entropy arising from each neighbouring species - we know the species is going to be a neighbour
+            # the bound ligand is always going to be a neighbour
+            omega = np.sqrt((neighbours_dict[neighbour] ** 3) * math.pi)
+            # orientational entropy arising from each neighbouring species
+            # - we know the species is going to be a neighbour
+            S_or_component = math.log(omega)
             S_or_component *= UAC.GAS_CONST
         S_or_total += S_or_component
-    # TODO for future releases -  implement a case for molecules with hydrogen bonds but to a lesser extent than water
+    # TODO for future releases
+    # implement a case for molecules with hydrogen bonds but to a lesser
+    # extent than water
     return S_or_total