14059657073

Committed 25 Mar 2025 12:31PM UTC coverage: 38.633% (+8.9%) from 29.725%

Build # 14059657073

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #73 from CCPBioSim/34-implement-python-logging

Implement Python Logging

Run Details

137 of 211 new or added lines in 9 files covered. (64.93%)

5 existing lines in 1 file now uncovered.

243 of 629 relevant lines covered (38.63%)

1.16 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

52.33

/CodeEntropy/calculations/EntropyFunctions.py

import logging
import math

# 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.calculations import ConformationFunctions as CONF
from CodeEntropy.calculations import UnitsAndConversions as UAC

logger = logging.getLogger(__name__)

# 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.

    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π

    Input
    -----
       lambdas : array of floats - eigenvalues of the covariance matrix
       temp: float - temperature

    Returns
    -------
       frequencies : array of floats - corresponding vibrational frequencies
    """
    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


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,
    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 (whole molecule)

    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
    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
            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
        else:
            S_vib_total = sum(S_components[6:])

    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


def conformational_entropy(
    data_container, dihedrals, bin_width, start, end, step, number_frames
):
    """
    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.

    Uses the adaptive enumeration method (AEM).

    Input
    -----
    dihedrals : array - array of dihedrals in the molecule
    Returns
    -------
       S_conf_total : float - conformational entropy
    """

    S_conf_total = 0

    # For each dihedral, identify the conformation in each frame
    num_dihedrals = len(dihedrals)
    conformation = np.zeros((num_dihedrals, number_frames))
    index = 0
    for dihedral in dihedrals:
        conformation[index] = CONF.assign_conformation(
            data_container, dihedral, number_frames, bin_width, start, end, step
        )
        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)
        S_conf_total += 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


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.

    σ 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?

    Input
    -----
      neighbours_dict :  dictionary - dictionary of neighbours for the molecule -
          should contain the type of neighbour molecule and the number of neighbour
          molecules of that species

    Returns
    -------
       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 neighbour in []:  # water molecules - call POSEIDON functions
            pass  # TODO temporary until function is written
        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

1	import logging	3✔
2	import math	3✔
3
4	# import matplotlib.pyplot as plt
5	# import MDAnalysis as mda
6	import numpy as np	3✔
7
8	# import pandas as pd
9	from numpy import linalg as la	3✔
10
11	from CodeEntropy.calculations import ConformationFunctions as CONF	3✔
12	from CodeEntropy.calculations import UnitsAndConversions as UAC	3✔
13
14	logger = logging.getLogger(__name__)	3✔
15
16	# import sys
17	# from ast import arg
18
19
20	# from CodeEntropy import NeighbourFunctions as NF
21
22
23	def frequency_calculation(lambdas, temp):	3✔
24	"""
25	Function to calculate an array of vibrational frequencies from the eigenvalues of
26	the covariance matrix.
27
28	Calculated from eq. (3) in Higham, S.-Y. Chou, F. Gräter and R. H. Henchman,
29	Molecular Physics, 2018, 116, 1965–1976//eq. (3) in A. Chakravorty, J. Higham
30	and R. H. Henchman, J. Chem. Inf. Model., 2020, 60, 5540–5551
31
32	frequency=sqrt(λ/kT)/2π
33
34	Input
35	-----
36	lambdas : array of floats - eigenvalues of the covariance matrix
37	temp: float - temperature
38
39	Returns
40	-------
41	frequencies : array of floats - corresponding vibrational frequencies
42	"""
43	pi = np.pi	3✔
44	# get kT in Joules from given temperature
45	kT = UAC.get_KT2J(temp)	3✔
46	logger.debug(f"Temperature: {temp}, kT: {kT}")	3✔
47
48	lambdas = np.array(lambdas) # Ensure input is a NumPy array	3✔
49	logger.debug(f"Eigenvalues (lambdas): {lambdas}")	3✔
50
51	# Check for negatives and raise an error if any are found
52	if np.any(lambdas < 0):	3✔
NEW 53	logger.error(f"Negative eigenvalues encountered: {lambdas[lambdas < 0]}")	×
54	raise ValueError(f"Negative eigenvalues encountered: {lambdas[lambdas < 0]}")	×
55
56	# Compute frequencies safely
57	frequencies = 1 / (2 * pi) * np.sqrt(lambdas / kT)	3✔
58	logger.debug(f"Calculated frequencies: {frequencies}")	3✔
59
60	return frequencies	3✔
61
62
63	def vibrational_entropy(matrix, matrix_type, temp, highest_level):	3✔
64	"""
65	Function to calculate the vibrational entropy for each level calculated from eq. (4)
66	in J. Higham, S.-Y. Chou, F. Gräter and R. H. Henchman, Molecular Physics, 2018,
67	116, 1965–1976 / eq. (2) in A. Chakravorty, J. Higham and R. H. Henchman,
68	J. Chem. Inf. Model., 2020, 60, 5540–5551.
69
70	Input
71	-----
72	matrix : matrix - force/torque covariance matrix
73	matrix_type: string
74	temp: float - temperature
75	highest_level: bool - is this the highest level of the heirarchy (whole molecule)
76
77	Returns
78	-------
79	S_vib_total : float - transvibrational/rovibrational entropy
80	"""
81	# N beads at a level => 3N x 3N covariance matrix => 3N eigenvalues
82	# Get eigenvalues of the given matrix and change units to SI units
83	lambdas = la.eigvals(matrix)	3✔
84	logger.debug(f"Eigenvalues (lambdas) before unit change: {lambdas}")	3✔
85	lambdas = UAC.change_lambda_units(lambdas)	3✔
86	logger.debug(f"Eigenvalues (lambdas) after unit change: {lambdas}")	3✔
87
88	# Calculate frequencies from the eigenvalues
89	frequencies = frequency_calculation(lambdas, temp)	3✔
90	logger.debug(f"Calculated frequencies: {frequencies}")	3✔
91
92	# Sort frequencies lowest to highest
93	frequencies = np.sort(frequencies)	3✔
94	logger.debug(f"Sorted frequencies: {frequencies}")	3✔
95
96	kT = UAC.get_KT2J(temp)	3✔
97	logger.debug(f"Temperature: {temp}, kT: {kT}")	3✔
98	exponent = UAC.PLANCK_CONST * frequencies / kT	3✔
99	logger.debug(f"Exponent values: {exponent}")	3✔
100	power_positive = np.power(np.e, exponent)	3✔
101	power_negative = np.power(np.e, -exponent)	3✔
102	logger.debug(f"Power positive values: {power_positive}")	3✔
103	logger.debug(f"Power negative values: {power_negative}")	3✔
104	S_components = exponent / (power_positive - 1) - np.log(1 - power_negative)	3✔
105	S_components = (	3✔
106	S_components * UAC.GAS_CONST
107	) # multiply by R - get entropy in J mol^{-1} K^{-1}
108	logger.debug(f"Entropy components: {S_components}")	3✔
109	# N beads at a level => 3N x 3N covariance matrix => 3N eigenvalues
110	if matrix_type == "force": # force covariance matrix	3✔
111	if highest_level: # whole molecule level - we take all frequencies into account	3✔
112	S_vib_total = sum(S_components)	3✔
113
114	# discard the 6 lowest frequencies to discard translation and rotation of the
115	# whole unit the overall translation and rotation of a unit is an internal
116	# motion of the level above
117	else:
118	S_vib_total = sum(S_components[6:])	×
119
120	else: # torque covariance matrix - we always take all values into account
121	S_vib_total = sum(S_components)	3✔
122
123	logger.debug(f"Total vibrational entropy: {S_vib_total}")	3✔
124
125	return S_vib_total	3✔
126
127
128	def conformational_entropy(	3✔
129	data_container, dihedrals, bin_width, start, end, step, number_frames
130	):
131	"""
132	Function to calculate conformational entropies using eq. (7) in Higham, S.-Y. Chou,
133	F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116, 1965–1976 / eq. (4) in
134	A. Chakravorty, J. Higham and R. H. Henchman, J. Chem. Inf. Model., 2020, 60,
135	5540–5551.
136
137	Uses the adaptive enumeration method (AEM).
138
139	Input
140	-----
141	dihedrals : array - array of dihedrals in the molecule
142	Returns
143	-------
144	S_conf_total : float - conformational entropy
145	"""
146
147	S_conf_total = 0	×
148
149	# For each dihedral, identify the conformation in each frame
150	num_dihedrals = len(dihedrals)	×
151	conformation = np.zeros((num_dihedrals, number_frames))	×
152	index = 0	×
153	for dihedral in dihedrals:	×
154	conformation[index] = CONF.assign_conformation(	×
155	data_container, dihedral, number_frames, bin_width, start, end, step
156	)
157	index += 1	×
158
NEW 159	logger.debug(f"Conformation matrix: {conformation}")	×
160
161	# For each frame, convert the conformation of all dihedrals into a state string
162	states = ["" for x in range(number_frames)]	×
163	for frame_index in range(number_frames):	×
164	for index in range(num_dihedrals):	×
165	states[frame_index] += str(conformation[index][frame_index])	×
166
NEW 167	logger.debug(f"States: {states}")	×
168
169	# Count how many times each state occurs, then use the probability to get the
170	# entropy
171	# entropy = sum over states p*ln(p)
172	values, counts = np.unique(states, return_counts=True)	×
173	for state in range(len(values)):	×
NEW 174	logger.debug(f"Unique states: {values}")	×
NEW 175	logger.debug(f"Counts: {counts}")	×
176	count = counts[state]	×
177	probability = count / number_frames	×
178	entropy = probability * np.log(probability)	×
179	S_conf_total += entropy	×
180
181	# multiply by gas constant to get the units J/mol/K
182	S_conf_total = -1 UAC.GAS_CONST	×
183
NEW 184	logger.debug(f"Total conformational entropy: {S_conf_total}")	×
185
186	return S_conf_total	×
187
188
189	def orientational_entropy(neighbours_dict):	3✔
190	"""
191	Function to calculate orientational entropies from eq. (10) in J. Higham, S.-Y.
192	Chou, F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116,3 1965–1976.
193	Number of orientations, Ω, is calculated using eq. (8) in J. Higham,
194	S.-Y. Chou, F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116,
195	3 1965–1976.
196
197	σ is assumed to be 1 for the molecules we're concerned with and hence,
198	max {1, (Nc^3π)^(1/2)} will always be (Nc^3π)^(1/2).
199
200	TODO future release - function for determing symmetry and symmetry numbers
201	maybe?
202
203	Input
204	-----
205	neighbours_dict : dictionary - dictionary of neighbours for the molecule -
206	should contain the type of neighbour molecule and the number of neighbour
207	molecules of that species
208
209	Returns
210	-------
211	S_or_total : float - orientational entropy
212	"""
213
214	# Replaced molecule with neighbour as this is what the for loop uses
215	S_or_total = 0	×
216	for neighbour in neighbours_dict: # we are going through neighbours	×
217	if neighbour in []: # water molecules - call POSEIDON functions	×
218	pass # TODO temporary until function is written	×
219	else:
220	# the bound ligand is always going to be a neighbour
221	omega = np.sqrt((neighbours_dict[neighbour] ** 3) * math.pi)	×
NEW 222	logger.debug(f"Omega for neighbour {neighbour}: {omega}")	×
223	# orientational entropy arising from each neighbouring species
224	# - we know the species is going to be a neighbour
225	S_or_component = math.log(omega)	×
NEW 226	logger.debug(	×
227	f"S_or_component (log(omega)) for neighbour {neighbour}: "
228	f"{S_or_component}"
229	)
230	S_or_component *= UAC.GAS_CONST	×
NEW 231	logger.debug(	×
232	f"S_or_component after multiplying by GAS_CONST for neighbour "
233	f"{neighbour}: {S_or_component}"
234	)
235	S_or_total += S_or_component	×
NEW 236	logger.debug(	×
237	f"S_or_total after adding component for neighbour {neighbour}: "
238	f"{S_or_total}"
239	)
240	# TODO for future releases
241	# implement a case for molecules with hydrogen bonds but to a lesser
242	# extent than water
243
NEW 244	logger.debug(f"Final total orientational entropy: {S_or_total}")	×
245
246	return S_or_total	×

CCPBioSim / CodeEntropy / 14059657073

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous