1- from ast import arg
2- import sys
31import math
2+ import sys
3+ from ast import arg
4+
5+ import matplotlib .pyplot as plt
6+ import MDAnalysis as mda
47import numpy as nmp
5- from numpy import linalg as la
68import pandas as pd
7- import MDAnalysis as mda
8- import matplotlib . pyplot as plt
9+ from numpy import linalg as la
10+
911from CodeEntropy import ConformationFunctions as CONF
1012from CodeEntropy import UnitsAndConversions as UAC
11- #from CodeEntropy import NeighbourFunctions as NF
1213
13- def frequency_calculation (lambdas ,temp ):
14+ # from CodeEntropy import NeighbourFunctions as NF
15+
16+
17+ def frequency_calculation (lambdas , temp ):
1418 """
1519 Function to calculate an array of vibrational frequencies from the eigenvalues of the
1620 covariance matrix.
17-
18- Calculated from eq. (3) in Higham, S.-Y. Chou, F. Gräter and R. H. Henchman, Molecular
19- Physics, 2018, 116, 1965–1976//eq. (3) in A. Chakravorty, J. Higham and R. H. Henchman,
21+
22+ Calculated from eq. (3) in Higham, S.-Y. Chou, F. Gräter and R. H. Henchman, Molecular
23+ Physics, 2018, 116, 1965–1976//eq. (3) in A. Chakravorty, J. Higham and R. H. Henchman,
2024 J. Chem. Inf. Model., 2020, 60, 5540–5551
21-
25+
2226 frequency=sqrt(λ/kT)/2π
23-
27+
2428 Input
2529 -----
2630 lambdas : array of floats - eigenvalues of the covariance matrix
@@ -33,18 +37,21 @@ def frequency_calculation(lambdas,temp):
3337 pi = nmp .pi
3438 # get kT in Joules from given temperature
3539 kT = UAC .get_KT2J (temp )
36- frequencies = 1 / ( 2 * pi )* nmp .sqrt (lambdas / kT )
40+ frequencies = 1 / ( 2 * pi ) * nmp .sqrt (lambdas / kT )
3741
3842 return frequencies
39- #END frequency_calculation
43+
44+
45+ # END frequency_calculation
46+
4047
4148def vibrational_entropy (matrix , matrix_type , temp , highest_level ):
4249 """
4350 Function to calculate the vibrational entropy for each level calculated from eq. (4) in
44- J. Higham, S.-Y. Chou, F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116,
51+ J. Higham, S.-Y. Chou, F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116,
4552 1965–1976 / eq. (2) in A. Chakravorty, J. Higham and R. H. Henchman,
4653 J. Chem. Inf. Model., 2020, 60, 5540–5551.
47-
54+
4855 Input
4956 -----
5057 matrix : matrix - force/torque covariance matrix
@@ -60,36 +67,43 @@ def vibrational_entropy(matrix, matrix_type, temp, highest_level):
6067 # Get eigenvalues of the given matrix and change units to SI units
6168 lambdas = la .eigvals (matrix )
6269 lambdas = UAC .change_lambda_units (lambdas )
63-
70+
6471 # Calculate frequencies from the eigenvalues
65- frequencies = frequency_calculation (lambdas ,temp )
66-
72+ frequencies = frequency_calculation (lambdas , temp )
73+
6774 # Sort frequencies lowest to highest
6875 frequencies = nmp .sort (frequencies )
6976
7077 kT = UAC .get_KT2J (temp )
71- exponent = UAC .PLANCK_CONST * frequencies / kT
72- power_positive = nmp .power (nmp .e ,exponent )
73- power_negative = nmp .power (nmp .e ,- exponent )
74- S_components = exponent / (power_positive - 1 ) - nmp .log (1 - power_negative )
75- S_components = S_components * UAC .GAS_CONST #multiply by R - get entropy in J mol^{-1} K^{-1}
78+ exponent = UAC .PLANCK_CONST * frequencies / kT
79+ power_positive = nmp .power (nmp .e , exponent )
80+ power_negative = nmp .power (nmp .e , - exponent )
81+ S_components = exponent / (power_positive - 1 ) - nmp .log (1 - power_negative )
82+ S_components = (
83+ S_components * UAC .GAS_CONST
84+ ) # multiply by R - get entropy in J mol^{-1} K^{-1}
7685 # N beads at a level => 3N x 3N covariance matrix => 3N eigenvalues
77- if matrix_type == ' force' : # force covariance matrix
78- if highest_level : # whole molecule level - we take all frequencies into account
86+ if matrix_type == " force" : # force covariance matrix
87+ if highest_level : # whole molecule level - we take all frequencies into account
7988 S_vib_total = sum (S_components )
80-
89+
8190 # discard the 6 lowest frequencies to discard translation and rotation of the whole unit
8291 # the overall translation and rotation of a unit is an internal motion of the level above
8392 else :
8493 S_vib_total = sum (S_components [6 :])
85-
86- else : # torque covariance matrix - we always take all values into account
94+
95+ else : # torque covariance matrix - we always take all values into account
8796 S_vib_total = sum (S_components )
88-
97+
8998 return S_vib_total
90- #END vibrational_entropy
9199
92- def conformational_entropy (data_container , dihedrals , bin_width , start , end , step , number_frames ):
100+
101+ # END vibrational_entropy
102+
103+
104+ def conformational_entropy (
105+ data_container , dihedrals , bin_width , start , end , step , number_frames
106+ ):
93107 """
94108 Function to calculate conformational entropies using eq. (7) in Higham, S.-Y. Chou,
95109 F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116, 1965–1976 / eq. (4) in
@@ -105,19 +119,21 @@ def conformational_entropy(data_container, dihedrals, bin_width, start, end, ste
105119 S_conf_total : float - conformational entropy
106120 """
107121
108- S_conf_total = 0
122+ S_conf_total = 0
109123
110124 # For each dihedral, identify the conformation in each frame
111125 num_dihedrals = len (dihedrals )
112126 conformation = nmp .zeros ((num_dihedrals , number_frames ))
113127 index = 0
114128 for dihedral in dihedrals :
115- conformation [index ] = CONF .assign_conformation (data_container , dihedral , number_frames , bin_width , start , end , step )
129+ conformation [index ] = CONF .assign_conformation (
130+ data_container , dihedral , number_frames , bin_width , start , end , step
131+ )
116132 index += 1
117133
118134 # For each frame, convert the conformation of all dihedrals into a state string
119- states = ['' for x in range (number_frames )]
120- for frame_index in range (number_frames ):
135+ states = ["" for x in range (number_frames )]
136+ for frame_index in range (number_frames ):
121137 for index in range (num_dihedrals ):
122138 states [frame_index ] += str (conformation [index ][frame_index ])
123139
@@ -134,20 +150,23 @@ def conformational_entropy(data_container, dihedrals, bin_width, start, end, ste
134150 S_conf_total *= - 1 * UAC .GAS_CONST
135151
136152 return S_conf_total
137- #END conformational_entropy
153+
154+
155+ # END conformational_entropy
156+
138157
139158def orientational_entropy (neighbours_dict ):
140159 """
141160 Function to calculate orientational entropies from eq. (10) in J. Higham, S.-Y. Chou,
142161 F. Gräter and R. H. Henchman, Molecular Physics, 2018, 116,3 1965–1976. Number of
143162 orientations, Ω, is calculated using eq. (8) in J. Higham, S.-Y. Chou, F. Gräter and
144- R. H. Henchman, Molecular Physics, 2018, 116,3 1965–1976.
163+ R. H. Henchman, Molecular Physics, 2018, 116,3 1965–1976.
145164
146165 σ is assumed to be 1 for the molecules we're concerned with and hence,
147166 max {1, (Nc^3*π)^(1/2)} will always be (Nc^3*π)^(1/2).
148-
167+
149168 TODO future release - function for determing symmetry and symmetry numbers maybe?
150-
169+
151170 Input
152171 -----
153172 neighbours_dict : dictionary - dictionary of neighbours for the molecule -
@@ -158,15 +177,21 @@ def orientational_entropy(neighbours_dict):
158177 -------
159178 S_or_total : float - orientational entropy
160179 """
161- S_or_total = 0
162- for neighbour in neighbours_dict : # we are going through neighbours
163- if molecule in [] : # water molecules - call POSEIDON functions
164- pass # TODO temporary until function is written
180+ S_or_total = 0
181+ for neighbour in neighbours_dict : # we are going through neighbours
182+ if molecule in []: # water molecules - call POSEIDON functions
183+ pass # TODO temporary until function is written
165184 else :
166- omega = nmp .sqrt ((neighbours_dict [molecule ]** 3 )* math .pi ) #the bound ligand is always going to be a neighbour
167- S_or_component = math .log (omega ) #orientational entropy arising from each neighbouring species - we know the species is going to be a neighbour
168- S_or_component *= UAC .GAS_CONST
169- S_or_total += S_or_component
170- #TODO for future releases - implement a case for molecules with hydrogen bonds but to a lesser extent than water
185+ omega = nmp .sqrt (
186+ (neighbours_dict [molecule ] ** 3 ) * math .pi
187+ ) # the bound ligand is always going to be a neighbour
188+ S_or_component = math .log (
189+ omega
190+ ) # orientational entropy arising from each neighbouring species - we know the species is going to be a neighbour
191+ S_or_component *= UAC .GAS_CONST
192+ S_or_total += S_or_component
193+ # TODO for future releases - implement a case for molecules with hydrogen bonds but to a lesser extent than water
171194 return S_or_total
172- #END orientational_entropy
195+
196+
197+ # END orientational_entropy
0 commit comments