FAQ | This is a LIVE service | Changelog

Skip to content
Snippets Groups Projects
Commit 1fab004c authored by M.D. Driver's avatar M.D. Driver
Browse files

update form of association constant in solvation equation.

parent b6e17090
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
solventmapcreator/solvationcalculation/solvationcalculator.py
+ 13
− 5
View file @ 1fab004c
......@@ -5,7 +5,7 @@ Script containing functions for calculating free energy of interaction for 2
SSIPs in solution
The energy for the interaction between 2 SSIPs is:
2*R * T * ln( [-1 + sqrt(1 + 4K*theta)]/[2K*theta]) - S1 - S2
2*R * T * ln( [-1 + sqrt(1 + 8K*theta)]/[4K*theta]) - S1 - S2
where R is gas constant,
T is temperature,
......@@ -21,7 +21,7 @@ First term in expression. Derivation:
and noting that [i_{free}] = [j_{free}]
probability that SSIP is unbound is:
P_{f} = [i_{free}]/[i] = -1 + sqrt(1 + 4K*theta)]/[2K*theta]
P_{f} = [i_{free}]/[i] = [-1 + sqrt(1 + 4K*theta)]/[2K*theta]
so the correction per SSIP is
R* T * ln( [-1 + sqrt(1 + 4K*theta)]/[2K*theta])
Since both SSIPs have a probility of being free then this becomes
......@@ -77,8 +77,8 @@ def calculate_binding_energy(epsilon_i, epsilon_j, temperature, theta):
This is in kJ/mol.
"""
association_constant = calculate_association_constant(epsilon_i, epsilon_j, temperature)
fraction_numerator = -1.0 + np.sqrt(1.0 + 4*association_constant*theta)
fraction_denominator = 2 * association_constant * theta
fraction_numerator = -1.0 + np.sqrt(1.0 + 8*association_constant*theta)
fraction_denominator = 4 * association_constant * theta
return 2 * (GAS_CONSTANT * CONVERSION_FROM_J_TO_KJ) * temperature * np.log(fraction_numerator/fraction_denominator)
def calculate_solvation_energy(epsilon, polynomial_coefficients):
......@@ -87,10 +87,18 @@ def calculate_solvation_energy(epsilon, polynomial_coefficients):
"""
return polyvalcalc.calculate_polynomial_values(epsilon, polynomial_coefficients)
def calculate_association_constant(epsilon_i, epsilon_j, temperature):
"""This returns the association constant
"""
k_ij = calculate_assoc_const(epsilon_i, epsilon_j, temperature)
k_ii = calculate_assoc_const(epsilon_i, epsilon_i, temperature)
k_jj = calculate_assoc_const(epsilon_j, epsilon_j, temperature)
total_k = 0.5* k_ij + 0.25 * k_ii + 0.25 *k_jj
return total_k
def calculate_assoc_const(epsilon_i, epsilon_j, temperature):
"""This calculates the association constant for a single interaction.
"""
interaction_energy = (epsilon_i * epsilon_j) + EVDW
exponent = -(interaction_energy * np.power(10,3))/(GAS_CONSTANT * temperature)
return 0.5 * np.exp(exponent)
solventmapcreator/test/solvationcalculationtest/solvationcalculatortest.py
+ 9
− 3
View file @ 1fab004c
......@@ -27,7 +27,7 @@ class SolvationCalculatorTestCase(unittest.TestCase):
def test_calcluate_association_energy_matrix(self):
"""Test to see if expected matrix is produced.
"""
expected_array = np.array([[-7.2500649680700651]])
expected_array = np.array([[-8.5730108585451354]])
actual_array = solvationcalculator.calcluate_association_energy_matrix([1.0],[1.0], 298.0, 1.0, np.array([0.5, 1.0]))
np.testing.assert_array_almost_equal(expected_array, actual_array)
actual_array2 = solvationcalculator.calcluate_association_energy_matrix(np.array([1.0]),np.array([1.0]), 298.0, 1.0, np.array([0.5, 1.0]))
......@@ -35,14 +35,14 @@ class SolvationCalculatorTestCase(unittest.TestCase):
def test_calculate_association_energy(self):
"""Test to see if expected_value is returned.
"""
expected_value = -7.2500649680700651
expected_value = -8.5730108585451354
actual_value = solvationcalculator.calculate_association_energy(1.0, 1.0,
298.0, 1.0, np.array([0.5, 1.0]))
self.assertAlmostEqual(expected_value, actual_value)
def test_calculate_binding_energy(self):
"""Test to see if expected binding energy is returned
"""
expected_value = -4.2500649680700662
expected_value = -5.5730108585451363
actual_value = solvationcalculator.calculate_binding_energy(1.0, 1.0, 298.0, 1.0)
self.assertAlmostEqual(expected_value, actual_value)
def test_calculate_solvation_energy(self):
......@@ -57,3 +57,9 @@ class SolvationCalculatorTestCase(unittest.TestCase):
expected_value = 3.200818570079262
actual_value = solvationcalculator.calculate_association_constant(1.0, 1.0, 298.0)
self.assertAlmostEqual(expected_value, actual_value)
def test_calculate_assoc_const(self):
"""Test to see if expected association constant is returned.
"""
expected_value = 3.200818570079262
actual_value = solvationcalculator.calculate_assoc_const(1.0, 1.0, 298.0)
self.assertAlmostEqual(expected_value, actual_value)
\ No newline at end of file
solventmapcreator/test/solvationcalculationtest/solvationmapgeneratortest.py
+ 6
− 6
View file @ 1fab004c
......@@ -61,9 +61,9 @@ class SolvationMapGeneratorTestCase(unittest.TestCase):
"""
expected_dict = {'x_data':np.array([[0, 0], [1, 1], [2, 2]]),
'y_data':np.array([[0, 1], [0, 1], [0, 1]]),
'z_data':np.array([[-4.332289, -5.403718],
[-5.403718, -5.757629],
[-6.475146, -6.167992]]),
'z_data':np.array([[-5.673465, -6.570665],
[-6.570665, -7.020214],
[-7.370673, -7.31281 ]]),
"x_label":"$\\beta$", "y_label":"$\\alpha$",
"plot_axis_range":(2.0, 0.0, 0.0, 1.0),
"figure_label":"water_solv_map",
......@@ -100,9 +100,9 @@ class SolvationMapGeneratorTestCase(unittest.TestCase):
"""
expected_dict = {'x_data':np.array([[0, 0], [1, 1], [2, 2]]),
'y_data':np.array([[0, 1], [0, 1], [0, 1]]),
'z_data':np.array([[-6.004792, -7.004792],
[-7.004792, -7.250065],
[-8.004792, -7.545754]]),
'z_data':np.array([[-7.396882, -8.217107],
[-8.217107, -8.573011],
[-8.93648 , -8.759837]]),
"x_label":"x_label", "y_label":"y_label",
"plot_axis_range":(2.0, 0.0, 0.0, 1.0),
"levels":[2*x -31.0 for x in range(32)]}
......
solventmapcreator/test/solvationcalculationtest/solvationplotinformationtest.py
+ 7
− 7
View file @ 1fab004c
......@@ -33,9 +33,9 @@ class SolvationPlotInformationTestCase(unittest.TestCase):
"""
expected_dict = {'x_data':np.array([[0, 0], [1, 1], [2, 2]]),
'y_data':np.array([[0, 1], [0, 1], [0, 1]]),
'z_data':np.array([[-6.004792, -7.004792],
[-7.004792, -7.250065],
[-8.004792, -7.545754]]),
'z_data':np.array([[-7.396882, -8.217107],
[-8.217107, -8.573011],
[-8.93648 , -8.759837]]),
"x_label":"x_label", "y_label":"y_label",
"plot_axis_range":(2.0, 0.0, 0.0, 1.0),
"levels":[2*x -31 for x in range(32)]}
......@@ -70,9 +70,9 @@ class SolvationPlotInformationTestCase(unittest.TestCase):
"""
expected_dict = {'x_data':np.array([[0, 0], [1, 1], [2, 2]]),
'y_data':np.array([[0, 1], [0, 1], [0, 1]]),
'z_data':np.array([[-6.004792, -7.004792],
[-7.004792, -7.250065],
[-8.004792, -7.545754]])}
'z_data':np.array([[-7.396882, -8.217107],
[-8.217107, -8.573011],
[-8.93648 , -8.759837]])}
actual_dict = solvationplotinformation.create_input_data_values(self.epsilon_i_list, self.epsilon_j_list, 298.0, 1.0, np.array([0.5, 1.0]))
self.assertListEqual(sorted(expected_dict.keys()), sorted(actual_dict.keys()))
for key in expected_dict.keys():
......@@ -94,7 +94,7 @@ class SolvationPlotInformationTestCase(unittest.TestCase):
def test_calculate_association_energy_matrix(self):
"""Test to see if expected matrix is returned.
"""
expected_array = np.array([[-7.2500649680700651]])
expected_array = np.array([[-8.5730108585451354]])
actual_array = solvationplotinformation.calculate_association_energy_matrix(np.array([1.0]),np.array([1.0]), 298.0, 1.0, np.array([0.5, 1.0]))
np.testing.assert_array_almost_equal(expected_array, actual_array)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment