13726590978

Committed 07 Mar 2025 06:07PM UTC coverage: 1.726%. First build

Build # 13726590978

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #45 from CCPBioSim/39-implement-ci-pipeline-pre-commit-hooks

Implement CI pipeline pre-commit hooks

Run Details

27 of 1135 new or added lines in 23 files covered. (2.38%)

53 of 3071 relevant lines covered (1.73%)

0.05 hits per line

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

49.12

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

CCPBioSim / CodeEntropy / 13726590978

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