diff --git a/solventmapcreator/solvationcalculation/solvationplotinformation.py b/solventmapcreator/solvationcalculation/solvationplotinformation.py
new file mode 100644
index 0000000000000000000000000000000000000000..a7db1e4cb7a8af0041a43379de625dab4987c67a
--- /dev/null
+++ b/solventmapcreator/solvationcalculation/solvationplotinformation.py
@@ -0,0 +1,51 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script for generating the input data required for plotting the contour map.
+
+@author: mark
+"""
+
+import logging
+import numpy as np
+import solventmapcreator.solvationcalculation.solvationcalculator as solvationcalculator
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.WARN)
+
+def create_plot_input_data(epsilon_i_list, epsilon_j_list, temperature, theta, polynomial_coefficients, x_label, y_label):
+    """This creates a dictionary containing the input data and labels for plotting.
+    """
+    input_data_values = create_input_data_values(epsilon_i_list, epsilon_j_list, temperature, theta, polynomial_coefficients)
+    input_data_labels = create_input_data_labels(epsilon_i_list, epsilon_j_list, x_label, y_label)
+    return {**input_data_values, **input_data_labels}
+
+def create_input_data_labels(epsilon_i_list, epsilon_j_list, x_label, y_label):
+    """This creates a dictionary containing the labels with annotations.
+    """
+    axis_range = get_axis_range(epsilon_i_list, epsilon_j_list)
+    return {"x_label":x_label, "y_label":y_label, "plot_axis_range":axis_range}
+
+def create_input_data_values(epsilon_i_list, epsilon_j_list, temperature, theta, polynomial_coefficients):
+    """This calculates the association energy matrix, and also generates the meshgrids ofthe x and y data.
+    """
+    association_matrix = calculate_association_energy_matrix(epsilon_i_list, epsilon_j_list, temperature, theta, polynomial_coefficients)
+    x_data, y_data = create_meshgrids(epsilon_i_list, epsilon_j_list)
+    return {"x_data":x_data, "y_data":y_data, "z_data":association_matrix}
+
+def get_axis_range(epsilon_i_list, epsilon_j_list):
+    """This gets the range of the values in the x and y data.
+    """
+    return ((epsilon_i_list.min(), epsilon_i_list.max()),
+            (epsilon_j_list.min(), epsilon_j_list.max()))
+
+def create_meshgrids(epsilon_i_list, epsilon_j_list):
+    """This creates the mesh grids of the x and y coordinates, to get the right dimensionality for plotting.
+    """
+    return np.meshgrid(epsilon_i_list, epsilon_j_list)
+
+def calculate_association_energy_matrix(epsilon_i_list, epsilon_j_list, temperature, theta, polynomial_coefficients):
+    """This creates the matrix of association constants.
+    """
+    return solvationcalculator.calcluate_association_energy_matrix(epsilon_i_list, epsilon_j_list, temperature, theta, polynomial_coefficients)
diff --git a/solventmapcreator/test/solvationcalculationtest/solvationplotinformationtest.py b/solventmapcreator/test/solvationcalculationtest/solvationplotinformationtest.py
new file mode 100644
index 0000000000000000000000000000000000000000..569bd61367fc80aa1a62ea6959d38f0400fc9a1b
--- /dev/null
+++ b/solventmapcreator/test/solvationcalculationtest/solvationplotinformationtest.py
@@ -0,0 +1,94 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Script containing tests for the solvationplotinformation script.
+
+@author: mark
+"""
+
+import logging
+import unittest
+import numpy as np
+import solventmapcreator.solvationcalculation.solvationplotinformation as solvationplotinformation
+
+logging.basicConfig()
+LOGGER = logging.getLogger(__name__)
+LOGGER.setLevel(logging.WARN)
+
+class SolvationPlotInformationTestCase(unittest.TestCase):
+    """Test Case for the solvationplotinformation script.
+    """
+    def setUp(self):
+        """Set up for tests.
+        """
+        self.epsilon_i_list = np.array([float(x) for x in range(3)])
+        self.epsilon_j_list = np.array([float(x) for x in range(2)])
+    def tearDown(self):
+        """clean up after tests.
+        """
+        del self.epsilon_i_list
+        del self.epsilon_j_list
+    def test_create_plot_input_data(self):
+        """Test to see if expected dictionary is produced.
+        """
+        expected_dict = {'x_data':np.array([[0, 1, 2],[0, 1, 2]]),
+                         'y_data':np.array([[0, 0, 0], [1, 1, 1]]),
+                         'z_data':np.array([[-4.004792, -3.004792],
+                                            [-3.004792, -1.250065],
+                                            [-2.004792,  0.454246]]),
+                         "x_label":"x_label", "y_label":"y_label",
+                         "plot_axis_range":((0.0, 2.0), (0.0, 1.0))}
+        actual_dict = solvationplotinformation.create_plot_input_data(self.epsilon_i_list, self.epsilon_j_list, 298.0, 1.0, np.array([0.5, 1.0]), "x_label", "y_label")
+        self.assertListEqual(sorted(expected_dict.keys()), sorted(actual_dict.keys()))
+        for key in expected_dict.keys():
+            if key == "plot_axis_range":
+                np.testing.assert_array_almost_equal(np.array(expected_dict[key]), np.array(actual_dict[key]))
+            elif "_data" in key:
+                np.testing.assert_array_almost_equal(expected_dict[key], actual_dict[key])
+            else:
+                self.assertEqual(expected_dict[key], actual_dict[key])
+    def test_create_input_data_labels(self):
+        """Test to see if expected dictionary is produced.
+        """
+        expected_dict = {"x_label":"x_label", "y_label":"y_label",
+                         "plot_axis_range":((0.0, 2.0), (0.0, 1.0))}
+        actual_dict = solvationplotinformation.create_input_data_labels(self.epsilon_i_list, self.epsilon_j_list, "x_label", "y_label")
+        self.assertListEqual(sorted(expected_dict.keys()), sorted(actual_dict.keys()))
+        for key in expected_dict.keys():
+            if key == "plot_axis_range":
+                np.testing.assert_array_almost_equal(np.array(expected_dict[key]), np.array(actual_dict[key]))
+            else:
+                self.assertEqual(expected_dict[key], actual_dict[key])
+    def test_create_input_data_values(self):
+        """Test to see if expected dictionary is produced.
+        """
+        expected_dict = {'x_data':np.array([[0, 1, 2],[0, 1, 2]]),
+                         'y_data':np.array([[0, 0, 0], [1, 1, 1]]),
+                         'z_data':np.array([[-4.004792, -3.004792],
+                                            [-3.004792, -1.250065],
+                                            [-2.004792,  0.454246]])}
+        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():
+            np.testing.assert_array_almost_equal(expected_dict[key], actual_dict[key])
+    def test_get_axis_range(self):
+        """Test to see if expected axis range is returned.
+        """
+        expected_range = ((0.0, 2.0),(0.0, 1.0))
+        actual_range = solvationplotinformation.get_axis_range(self.epsilon_i_list, self.epsilon_j_list)
+        np.testing.assert_array_almost_equal(np.array(expected_range), np.array(actual_range))
+    def test_create_meshgrids(self):
+        """Test to see if expected arrays are produced.
+        """
+        expected_x = np.array([[0, 1, 2],[0, 1, 2]])
+        expected_y = np.array([[0, 0, 0], [1, 1, 1]])
+        actual_x, actual_y = solvationplotinformation.create_meshgrids(self.epsilon_i_list, self.epsilon_j_list)
+        np.testing.assert_array_almost_equal(expected_x, actual_x)
+        np.testing.assert_array_almost_equal(expected_y, actual_y)
+    def test_calculate_association_energy_matrix(self):
+        """Test to see if expected matrix is returned.
+        """
+        expected_array = np.array([[-1.2500649680700651]])
+        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)
+    
diff --git a/solventmapcreator/test/solvationmapcreatortests.py b/solventmapcreator/test/solvationmapcreatortests.py
index d3808dc307941ddce2b64a0efa6ce56316e54b82..21684a981782d1b0ca0aab742284aa5760b07ba2 100644
--- a/solventmapcreator/test/solvationmapcreatortests.py
+++ b/solventmapcreator/test/solvationmapcreatortests.py
@@ -13,6 +13,7 @@ from solventmapcreator.test.iotest.solventxmlreadertest import SolventXMLReaderT
 from solventmapcreator.test.iotest.polynomialdatawritertest import PolynomialDataWriterTestCase
 from solventmapcreator.test.iotest.polynomialdatareadertest import PolynomialDataReaderTestCase
 from solventmapcreator.test.solvationcalculationtest.solvationcalculatortest import SolvationCalculatorTestCase
+from solventmapcreator.test.solvationcalculationtest.solvationplotinformationtest import SolvationPlotInformationTestCase
 from solventmapcreator.test.solvationcalculationtest.fractionaloccupancycalculatortest import FractionalOccupancyCalculatorTestCase
 from solventmapcreator.test.polynomialanalysistest.polynomialdataanalysistest import PolynomialDataAnalysisTestCase
 from solventmapcreator.test.polynomialanalysistest.polynomialcomparisontest import PolynomialComparisonTestCase
@@ -25,6 +26,7 @@ IO_TEST_CASES = [SolvationEnergyReaderTestCase, SolventXMLReaderTestCase,
                  PolynomialDataWriterTestCase, PolynomialDataReaderTestCase]
 
 SOLVATION_CALCULATION_TEST_CASES = [SolvationCalculatorTestCase,
+                                    SolvationPlotInformationTestCase,
                                     FractionalOccupancyCalculatorTestCase]
 
 POLYNOMIAL_ANALYSIS_TEST_CASES = [PolynomialDataAnalysisTestCase,