diff --git a/solventmapcreator/solvationcalculation/solvationcalculator.py b/solventmapcreator/solvationcalculation/solvationcalculator.py
new file mode 100644
index 0000000000000000000000000000000000000000..1ed157c3e096d4480be80a396505b8a408a2c603
--- /dev/null
+++ b/solventmapcreator/solvationcalculation/solvationcalculator.py
@@ -0,0 +1,73 @@
+#!/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)
diff --git a/solventmapcreator/test/solvationcalculationtest/solvationcalculatortest.py b/solventmapcreator/test/solvationcalculationtest/solvationcalculatortest.py
new file mode 100644
index 0000000000000000000000000000000000000000..a479a06940b8a3466bf39a51e5e32f4c81b9983b
--- /dev/null
+++ b/solventmapcreator/test/solvationcalculationtest/solvationcalculatortest.py
@@ -0,0 +1,51 @@
+#!/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)
diff --git a/solventmapcreator/test/solvationmapcreatortests.py b/solventmapcreator/test/solvationmapcreatortests.py
index 6a71d92b381f08c22bf5bd63411287a8795a1d77..151dbdc01f853428817d1921e08877a83c98626c 100644
--- a/solventmapcreator/test/solvationmapcreatortests.py
+++ b/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__':