diff --git a/solventmapcreator/solvationcalculation/solvationcalculator.py b/solventmapcreator/solvationcalculation/solvationcalculator.py
index fdebc4aa8dce0ad653690d9a62275e48f5490f46..539121287ac894eea06db55522a4f6143d6979a7 100644
--- a/solventmapcreator/solvationcalculation/solvationcalculator.py
+++ b/solventmapcreator/solvationcalculation/solvationcalculator.py
@@ -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)
diff --git a/solventmapcreator/test/solvationcalculationtest/solvationcalculatortest.py b/solventmapcreator/test/solvationcalculationtest/solvationcalculatortest.py
index 9311b40ed4191f793797ff3212373f5f14993f59..227570a2e651babc5cb49b2d433df2f6884b39a8 100644
--- a/solventmapcreator/test/solvationcalculationtest/solvationcalculatortest.py
+++ b/solventmapcreator/test/solvationcalculationtest/solvationcalculatortest.py
@@ -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
diff --git a/solventmapcreator/test/solvationcalculationtest/solvationmapgeneratortest.py b/solventmapcreator/test/solvationcalculationtest/solvationmapgeneratortest.py
index 022be01e3de247626d3d01f59c6971dfae3d55f7..952beec26c8ec8f3caf5a3b8bcac30d751d4fcaa 100644
--- a/solventmapcreator/test/solvationcalculationtest/solvationmapgeneratortest.py
+++ b/solventmapcreator/test/solvationcalculationtest/solvationmapgeneratortest.py
@@ -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)]}
diff --git a/solventmapcreator/test/solvationcalculationtest/solvationplotinformationtest.py b/solventmapcreator/test/solvationcalculationtest/solvationplotinformationtest.py
index 6ff6413f4d8cbdbe6da928da4c819824030a265c..90c9f39072280183fc7ddfd398739db598717f33 100644
--- a/solventmapcreator/test/solvationcalculationtest/solvationplotinformationtest.py
+++ b/solventmapcreator/test/solvationcalculationtest/solvationplotinformationtest.py
@@ -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)