13769110208

Committed 10 Mar 2025 03:57PM UTC coverage: 1.774% (+0.05%) from 1.726%

Build # 13769110208

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #53 from CCPBioSim/7-nan-error-handling

Add error handling for NaN and invalid values in calculations:

Run Details

2 of 33 new or added lines in 3 files covered. (6.06%)

3 existing lines in 2 files now uncovered.

55 of 3101 relevant lines covered (1.77%)

0.05 hits per line

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

50.0

/CodeEntropy/EntropyFunctions.py

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

    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)

    lambdas = np.array(lambdas)  # Ensure input is a NumPy array

    # Check for negatives and raise an error if any are found
    if np.any(lambdas < 0):
        raise ValueError(f"Negative eigenvalues encountered: {lambdas[lambdas < 0]}")

    # Compute frequencies safely
    frequencies = 1 / (2 * pi) * np.sqrt(lambdas / kT)

    return frequencies


# END frequency_calculation


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)
    lambdas = UAC.change_lambda_units(lambdas)

    # Calculate frequencies from the eigenvalues
    frequencies = frequency_calculation(lambdas, temp)

    # Sort frequencies lowest to highest
    frequencies = np.sort(frequencies)

    kT = UAC.get_KT2J(temp)
    exponent = UAC.PLANCK_CONST * frequencies / kT
    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}
    # 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)

    return S_vib_total


# END vibrational_entropy


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

    # 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])

    # 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)):
        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

    return S_conf_total


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

    σ 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)
            # 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
    return S_or_total


# END orientational_entropy

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

CCPBioSim / CodeEntropy / 13769110208

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