added class for the calculation of the solvation energy.

solventmapcreator/solvationcalculation/solvationcalculator.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+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
+
+where R is gas constant,
+T is temperature,
+K is association constant,
+theta is fractional occupancy of the phase,
+S1, S2 are the solvation energy of SSIP 1 and 2 respectively.
+
+@author: mark
+"""
+
+import logging
+import numpy as np
+import numpy.polynomial.polynomial as poly
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.WARN)
+
+CONVERSION_FROM_J_TO_KJ = 0.001
+"""
+value for the gas constant, R in JK^-1mol^-1.
+https://en.wikipedia.org/wiki/Gas_constant
+Multiply by conversion factor to express this in kJK^-1mol^-1, so that
+output values are in kJmol^-1.
+"""
+GAS_CONSTANT = 8.3144598
+
+"""
+value of the van der Waals interaction of 2 SSIPs interacting in a
+pairwise manner. As detailed in DOI: http://dx.doi.org/10.1039/c3sc22124e.
+"""
+EVDW = -5.6
+
+def calculate_association_energy(epsilon_i, epsilon_j, temperature, theta, polynomial_coefficients):
+    """This calculates the energy described above
+    """
+    solvation_energy_i = calculate_solvation_energy(epsilon_i, polynomial_coefficients)
+    solvation_energy_j = calculate_solvation_energy(epsilon_j, polynomial_coefficients)
+    
+    binding_energy = calculate_binding_energy(epsilon_i, epsilon_j, temperature, theta)
+    return binding_energy + solvation_energy_i + solvation_energy_j
+
+def calculate_binding_energy(epsilon_i, epsilon_j, temperature, theta):
+    """This evaluates the first term in the association energy calculation.
+    This corresponds to the energy of binding between the 2 SSIPs. 
+    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
+    return 2 * (GAS_CONSTANT * CONVERSION_FROM_J_TO_KJ) * temperature * np.log(fraction_numerator/fraction_denominator)
+
+def calculate_solvation_energy(epsilon, polynomial_coefficients):
+    """This calcultes the solvation energy for the given SSIP value, using the
+    polynomial coefficients given.
+    """
+    return poly.polyval(epsilon, polynomial_coefficients)
+
+
+def calculate_association_constant(epsilon_i, epsilon_j, temperature):
+    """This returns the association constant
+    """
+    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

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script for testing solvationcalculator script.
+
+@author: mark
+"""
+
+import logging
+import unittest
+import numpy as np
+import solventmapcreator.solvationcalculation.solvationcalculator as solvationcalculator
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.WARN)
+
+class SolvationCalculatorTestCase(unittest.TestCase):
+    """Test case for the solvationcalculator script
+    """
+    def setUp(self):
+        """Set up for tests.
+        """
+    def tearDown(self):
+        """tear down after tests.
+        """
+    def test_calculate_association_energy(self):
+        """Test to see if expected_value is returned.
+        """
+        expected_value = -1.2500649680700651
+        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
+        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):
+        """Test to see if expected solvation energy is returned.
+        """
+        expected_value = 1.0
+        actual_value = solvationcalculator.calculate_solvation_energy(0.5, np.array([0.5, 1.0]))
+        self.assertAlmostEqual(expected_value, actual_value)
+    def test_calculate_association_constant(self):
+        """Test to see if expected value is returned.
+        """
+        expected_value = 3.200818570079262
+        actual_value = solvationcalculator.calculate_association_constant(1.0, 1.0, 298.0)
+        self.assertAlmostEqual(expected_value, actual_value)

solventmapcreator/test/solvationmapcreatortests.py

@@ -9,6 +9,7 @@ Script for running the tests.
 import logging
 import unittest
 from solventmapcreator.test.iotest.solvationenergyreadertest import SolvationEnergyReaderTestCase
+from solventmapcreator.test.solvationcalculationtest.solvationcalculatortest import SolvationCalculatorTestCase
 
 logging.basicConfig()
 LOGGER = logging.getLogger(__name__)
@@ -16,6 +17,7 @@ LOGGER.setLevel(logging.WARN)
 
 IO_TEST_CASES = [SolvationEnergyReaderTestCase]
 
+SOLVATION_CALCULATION_TEST_CASES = [SolvationCalculatorTestCase]
 def test_suite():
     """Function creates a test suite and then loads all the tests from the
     different test cases.
@@ -26,6 +28,9 @@ def test_suite():
     for test_case in IO_TEST_CASES:
         LOGGER.debug("Adding %s", test_case)
         suite.addTests(loader.loadTestsFromTestCase(test_case))
+    for test_case in SOLVATION_CALCULATION_TEST_CASES:
+        LOGGER.debug("Adding %s", test_case)
+        suite.addTests(loader.loadTestsFromTestCase(test_case))
     return suite
 
 if __name__ == '__main__':