diff --git a/examples/ideal_gas_law/ideal_gas_law.ipynb b/examples/ideal_gas_law/ideal_gas_law.ipynb
index f55e986..4304eb6 100644
--- a/examples/ideal_gas_law/ideal_gas_law.ipynb
+++ b/examples/ideal_gas_law/ideal_gas_law.ipynb
@@ -450,7 +450,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.2"
+   "version": "3.10.12"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,
diff --git a/examples/simple_examples/molecular_dynamics.ipynb b/examples/simple_examples/molecular_dynamics.ipynb
index e111b50..c748151 100644
--- a/examples/simple_examples/molecular_dynamics.ipynb
+++ b/examples/simple_examples/molecular_dynamics.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -23,7 +23,8 @@
    "source": [
     "def md_simulation(number_of_particles, temperature, box_length, number_of_steps, sample_frequency):\n",
     "    # Initialise the system\n",
-    "    system = md.initialise(number_of_particles, temperature, box_length, 'square')\n",
+    "    constants=[[1.363e-134, 9.273e-78],[1.365e-130, 9.278e-77],[1.368e-130, 9.278e-77]]\n",
+    "    system = md.initialise(number_of_particles, temperature, box_length, 'square', constants=constants)\n",
     "    # This sets the sampling class\n",
     "    sample_system = sample.MaxBolt(system)# Start at time 0\n",
     "    system.time = 0\n",
@@ -63,7 +64,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "md_simulation(20, 1000, 100, 5000, 100)"
+    "system = md_simulation(20, 1000, 100, 5000, 25)"
    ]
   },
   {
@@ -90,7 +91,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.2"
+   "version": "3.12.4"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,
diff --git a/examples/simple_examples/monte_carlo.ipynb b/examples/simple_examples/monte_carlo.ipynb
index f728d8b..d92a6b1 100644
--- a/examples/simple_examples/monte_carlo.ipynb
+++ b/examples/simple_examples/monte_carlo.ipynb
@@ -22,7 +22,8 @@
    "source": [
     "def mc_simulation(number_of_particles, temperature, box_length, number_of_steps, sample_frequency):\n",
     "    # Initialise the system placing the particles on a square lattice\n",
-    "    system = mc.initialise(number_of_particles, temperature, box_length, 'square')\n",
+    "    constants=[[1.363e-134, 9.273e-78],[1.363e-134, 9.273e-78]]\n",
+    "    system = mc.initialise(number_of_particles, temperature, box_length, 'square', constants=constants, mixing = False)\n",
     "    # This sets the sampling class as Energy, which shows the energy of the system\n",
     "    sample_system = sample.Energy(system)\n",
     "    # Compute the energy of the system\n",
@@ -71,7 +72,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "system = mc_simulation(100, 273.15, 45, 5000, 25)"
+    "system = mc_simulation(50, 273.15, 45, 1000, 25)"
    ]
   },
   {
@@ -105,7 +106,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.2"
+   "version": "3.12.4"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,
diff --git a/pylj/forcefields.py b/pylj/forcefields.py
index ddc59e2..e54fa41 100644
--- a/pylj/forcefields.py
+++ b/pylj/forcefields.py
@@ -1,61 +1,6 @@
 import numpy as np
 
 
-class lennard_jones(object):
-    r"""Calculate the energy or force for a pair of particles using the
-    Lennard-Jones (A/B variant) forcefield.
-
-    Parameters
-    ----------
-    constants: float, array_like
-        An array of length two consisting of the A and B parameters for the
-        12-6 Lennard-Jones function
-
-    """
-    def __init__(self, constants):
-        self.a = constants[0]
-        self.b = constants[1]
-
-    def energy(self, dr ):
-        r"""Calculate the energy for a pair of particles using the
-        Lennard-Jones (A/B variant) forcefield.
-
-        .. math::
-            E = \frac{A}{dr^{12}} - \frac{B}{dr^6}
-
-        Attributes:
-        ----------
-        dr (float): The distance between particles.
-
-        Returns
-        -------
-        float: array_like
-        The potential energy between the particles.
-        """
-        self.energy = self.a * np.power(dr, -12) - (self.b * np.power(dr, -6))
-        return self.energy
-    
-    def force(self, dr):
-        r"""Calculate the force for a pair of particles using the
-        Lennard-Jones (A/B variant) forcefield.
-
-        .. math::
-            f = \frac{12A}{dr^{13}} - \frac{6B}{dr^7}
-
-        Attributes:
-        ----------
-        dr (float): The distance between particles.
-        
-        Returns
-        -------
-        float: array_like
-        The force between the particles.
-        """
-        self.force = 12 * self.a * np.power(dr, -13) - (6 * self.b * np.power(dr, -7))
-        return self.force
-
-
-
 class lennard_jones_sigma_epsilon(object):
     r"""Calculate the energy or force for a pair of particles using the
     Lennard-Jones (sigma/epsilon variant) forcefield.
@@ -68,8 +13,12 @@ class lennard_jones_sigma_epsilon(object):
 
     """
     def __init__(self, constants):
+        if len(constants) != 2:
+            raise IndexError(f'There should be two constants per set, not {len(constants)}')
+        
         self.sigma = constants[0]
         self.epsilon = constants[1]
+        self.point_size = 1.3e10 * (self.sigma*(2**(1/6)))
     
     def energy(self, dr):
         r"""Calculate the energy for a pair of particles using the
@@ -110,8 +59,70 @@ def force(self, dr):
         self.force = 48 * self.epsilon * np.power(self.sigma, 12) * np.power(
             dr, -13) - (24 * self.epsilon * np.power(self.sigma, 6) * np.power(dr, -7))
         return self.force
+    
+    def mixing(self, constants_2):
+        r""" Calculates mixing for two sets of constants
+        
+        ..math::
+            \sigma_{12} = \frac{\sigma_1 + \sigma_2}{2}
+            \epsilon{12} = \sqrt{\epsilon_1 * \epsilon_2}
+        
+        Parameters:
+        ----------
+        constants_2: float, array_like
+            The second set of constantss
+        """
+        sigma2 = constants_2[0]
+        epsilon2 = constants_2[1]
+        self.sigma = (self.sigma+sigma2)/2
+        self.epsilon = np.sqrt(self.epsilon * epsilon2)
 
 
+class lennard_jones(lennard_jones_sigma_epsilon):
+    r"""Converts a/b variant values to sigma/epsilon variant
+    then maps to lennard_jones_sigma_epsilon class
+
+    ..math::
+        \sigma = \frac{a}{b}^(\frac{1}{6})
+        \sigma = \frace{b^2}{4*a}
+
+    Parameters
+    ----------
+    constants: float, array_like
+        An array of length two consisting of the A and B
+        parameters for the 12-6 Lennard-Jones function
+    """   
+    def __init__(self, constants):
+        if len(constants) != 2:
+            raise IndexError(f'There should be two constants per set, not {len(constants)}')
+        self.a = constants[0]
+        self.b = constants[1]
+        sigma = (self.a / self.b)**(1/6)
+        epsilon = (self.b**2)/(4*self.a)
+        super().__init__([sigma, epsilon])
+
+    def mixing(self, constants_2):
+        r"""Converts second set of a/b constants into sigma/epsilon
+        for use in mixing method. Then converts changed self sigma/epsilon
+        values back to a/b
+
+        ..math::
+            a = 4*\epsilon*(\sigma^12)
+            b = 4*\epsilon*(\sigma^6)
+
+        Parameters:
+        ----------
+        constants_2: float, array_like
+            The second set of constantss
+        """
+        a2 = constants_2[0]
+        b2 = constants_2[1]
+        sigma2 = (a2 / b2)**(1/6)
+        epsilon2 = (b2**2)/(4*a2)
+        super().mixing([sigma2,epsilon2])
+        self.a = 4 * self.epsilon * (self.sigma**12)
+        self.b = 4 * self.epsilon * (self.sigma**6)
+
 
 class buckingham(object):
     r""" Calculate the energy or force for a pair of particles using the
@@ -125,9 +136,12 @@ class buckingham(object):
 
     """
     def __init__(self, constants):
+        if len(constants) != 3:
+            raise IndexError(f'There should be three constants per set, not {len(constants)}')
         self.a = constants[0]
         self.b = constants[1]
         self.c = constants[2]
+        self.point_size = 8 # Needs better solution relevant to constants
     
     def energy(self, dr):
         r"""Calculate the energy for a pair of particles using the
@@ -145,7 +159,13 @@ def energy(self, dr):
         float: array_like
         The potential energy between the particles.
         """
-        self.energy = self.a * np.exp(- np.multiply(self.b, dr)) - self.c / np.power(dr, 6)
+        energy = self.a * np.exp(- np.multiply(self.b, dr)) - self.c / np.power(dr, 6)
+        # Cut out infinite values where r = 0
+        if type(dr) != float:
+            energy = np.array(energy)
+            energy[np.where(energy > 10e300)] = 0
+            energy[np.where(energy < -10e300)] = 0
+        self.energy = energy
         return self.energy
     
     def force(self, dr):
@@ -164,9 +184,31 @@ def force(self, dr):
         float: array_like
         The force between the particles.
         """
-        self.force = self.a * self.b * np.exp(- np.multiply(self.b, dr)) - 6 * self.c / np.power(dr, 7)
+        force = self.a * self.b * np.exp(- np.multiply(self.b, dr)) - 6 * self.c / np.power(dr, 7)
+        # Cut out infinite values where r = 0
+        if type(dr) != float:
+            force = np.array(force)
+            force[np.where(force > 10e300)] = 0
+            force[np.where(force < -10e300)] = 0
+        self.force = force
         return self.force
 
+    def mixing(self, constants2):
+        r""" Calculates mixing for two sets of constants
+        
+        ..math::
+            a_{12} = \sqrt{a_1 * a_2}
+            b_{12} = \sqrt{b_1 * b_2}
+            c_{12} = \sqrt{c_1 * c_2}
+        
+        Attributes:
+        ----------
+        constants_2: float, array_like
+            The second set of constantss
+        """
+        self.a = np.sqrt(self.a*constants2[0])
+        self.b = np.sqrt(self.b*constants2[1])
+        self.c = np.sqrt(self.c*constants2[2])
 
 
 class square_well(object):
@@ -183,11 +225,13 @@ class square_well(object):
         
     '''
     def __init__(self, constants, max_val=np.inf):
-
+        if len(constants) != 3:
+            raise IndexError(f'There should be three constants per set, not {len(constants)}')
         self.epsilon = constants[0]
         self.sigma = constants[1]
         self.lamda = constants[2] #Spelling as lamda not lambda to avoid calling python lambda function
         self.max_val = max_val
+        self.point_size = 10
 
     def energy(self, dr):
         r'''Calculate the energy for a pair of particles using a
diff --git a/pylj/mc.py b/pylj/mc.py
index 8dc41cc..2532148 100644
--- a/pylj/mc.py
+++ b/pylj/mc.py
@@ -8,8 +8,8 @@ def initialise(
     box_length,
     init_conf,
     mass=39.948,
-    constants=[1.363e-134, 9.273e-78],
-    forcefield=ff.lennard_jones,
+    constants=[[1.363e-134, 9.273e-78]],
+    forcefield=ff.lennard_jones
 ):
     """Initialise the particle positions (this can be either as a square or
     random arrangement), velocities (based on the temperature defined, and #
@@ -48,7 +48,7 @@ def initialise(
         constants,
         forcefield,
         mass,
-        init_conf=init_conf,
+        init_conf=init_conf
     )
     system.particles["xvelocity"] = 0
     system.particles["yvelocity"] = 0
diff --git a/pylj/md.py b/pylj/md.py
index 86e9268..f640f2d 100644
--- a/pylj/md.py
+++ b/pylj/md.py
@@ -10,8 +10,8 @@ def initialise(
     init_conf,
     timestep_length=1e-14,
     mass=39.948,
-    constants=[1.363e-134, 9.273e-78],
-    forcefield=ff.lennard_jones,
+    constants=[[1.363e-134, 9.273e-78]],
+    forcefield=ff.lennard_jones
 ):
     """Initialise the particle positions (this can be either as a square or
     random arrangement), velocities (based on the temperature defined, and
@@ -53,7 +53,7 @@ def initialise(
         forcefield,
         mass,
         init_conf=init_conf,
-        timestep_length=timestep_length,
+        timestep_length=timestep_length
     )
     v = np.random.rand(system.particles.size, 2, 12)
     v = np.sum(v, axis=2) - 6.0
@@ -153,6 +153,7 @@ def sample(particles, box_length, initial_particles, system):
         system.cut_off,
         system.constants,
         system.forcefield,
+        system.mass
     )
     msd_new = calculate_msd(particles, initial_particles, box_length)
     system.pressure_sample = np.append(system.pressure_sample, pressure_new)
@@ -329,37 +330,6 @@ def compute_force(particles, box_length, cut_off, constants, forcefield, mass):
     return part, dist, forces, energies
 
 
-def compute_energy(particles, box_length, cut_off, constants, forcefield):
-    r"""Calculates the total energy of the simulation.
-    Parameters
-    ----------
-    particles: util.particle_dt, array_like
-        Information about the particles.
-    box_length: float
-        Length of a single dimension of the simulation square, in Angstrom.
-    cut_off: float
-        The distance greater than which the energies between particles is
-        taken as zero.
-    constants: float, array_like (optional)
-        The constants associated with the particular forcefield used, e.g. for
-        the function forcefields.lennard_jones, theses are [A, B]
-    forcefield: function (optional)
-        The particular forcefield to be used to find the energy and forces.
-    Returns
-    -------
-    util.particle_dt, array_like
-        Information about particles, with updated accelerations and forces.
-    float, array_like
-        Current distances between pairs of particles in the simulation.
-    float, array_like
-        Current energies between pairs of particles in the simulation.
-    """
-    dist, energies = heavy.compute_energy(
-        particles, box_length, cut_off, constants, forcefield
-    )
-    return dist, energies
-
-
 def heat_bath(particles, temperature_sample, bath_temperature):
     r"""Rescales the velocities of the particles in the system to control the
     temperature of the simulation. Thereby allowing for an NVT ensemble. The
diff --git a/pylj/pairwise.py b/pylj/pairwise.py
index 72f2f60..06b885c 100644
--- a/pylj/pairwise.py
+++ b/pylj/pairwise.py
@@ -39,17 +39,32 @@ def compute_force(particles, box_length, cut_off, constants, forcefield, mass):
         (particles["xacceleration"].size - 1) * particles["xacceleration"].size / 2
     )
     forces = np.zeros(pairs)
-    distances = np.zeros(pairs)
     energies = np.zeros(pairs)
     atomic_mass_unit = 1.660539e-27  # kilograms
     mass_amu = mass  # amu
     mass_kg = mass_amu * atomic_mass_unit  # kilograms
-    distances, dx, dy = heavy.dist(
-        particles["xposition"], particles["yposition"], box_length
+    distances, dx, dy, pair_types = heavy.dist(
+        particles["xposition"], particles["yposition"], box_length, particles['types']
     )
-    ff = forcefield(constants)
-    forces = ff.force(distances)
-    energies = ff.energy(distances)
+    unique_pairs = list(set(pair_types))
+    for pair in unique_pairs:
+        type_distances = distances.copy()
+        for i in range(len(distances)):
+            if pair != pair_types[i]:
+                type_distances[i] = 0
+        type_1 = pair.split(',')[0]
+        type_2 = pair.split(',')[1]
+        constants_1 = np.array(constants[int(pair.split(',')[0])])
+        ff = forcefield(constants_1)
+        if type_1 != type_2:
+            constants_2 = np.array(constants[int(pair.split(',')[1])])
+            ff.mixing(constants_2)
+        type_forces = ff.force(type_distances)
+        type_energies = ff.energy(distances)
+        type_forces = np.nan_to_num(type_forces)
+        type_energies = np.nan_to_num(type_energies)
+        forces+=type_forces
+        energies+=type_energies
     forces[np.where(distances > cut_off)] = 0.0
     energies[np.where(distances > cut_off)] = 0.0
     particles = update_accelerations(particles, forces, mass_kg, dx, dy, distances)
@@ -95,11 +110,14 @@ def update_accelerations(particles, f, m, dx, dy, dr):
     k = 0
     for i in range(0, particles.size - 1):
         for j in range(i + 1, particles.size):
-            particles["xacceleration"][i] += second_law(f[k], m, dx[k], dr[k])
-            particles["yacceleration"][i] += second_law(f[k], m, dy[k], dr[k])
-            particles["xacceleration"][j] -= second_law(f[k], m, dx[k], dr[k])
-            particles["yacceleration"][j] -= second_law(f[k], m, dy[k], dr[k])
+
+            particles["xacceleration"][i] += second_law(f[k], m, dx[k], dr[k]) if f[k]!=0 else 0
+            particles["yacceleration"][i] += second_law(f[k], m, dy[k], dr[k]) if f[k]!=0 else 0
+            
+            particles["xacceleration"][j] -= second_law(f[k], m, dx[k], dr[k]) if f[k]!=0 else 0
+            particles["yacceleration"][j] -= second_law(f[k], m, dy[k], dr[k]) if f[k]!=0 else 0
             k += 1
+        
     return particles
 
 
@@ -178,47 +196,8 @@ def lennard_jones_force(A, B, dr):
     return 12 * A * np.power(dr, -13) - 6 * B * np.power(dr, -7)
 
 
-def compute_energy(particles, box_length, cut_off, constants, forcefield):
-    r"""Calculates the total energy of the simulation.
-    Parameters
-    ----------
-    particles: util.particle_dt, array_like
-        Information about the particles.
-    box_length: float
-        Length of a single dimension of the simulation square, in Angstrom.
-    cut_off: float
-        The distance greater than which the energies between particles is
-        taken as zero.
-    constants: float, array_like (optional)
-        The constants associated with the particular forcefield used, e.g. for
-        the function forcefields.lennard_jones, theses are [A, B]
-    forcefield: function (optional)
-        The particular forcefield to be used to find the energy and forces.
-    Returns
-    -------
-    util.particle_dt, array_like
-        Information about particles, with updated accelerations and forces.
-    float, array_like
-        Current distances between pairs of particles in the simulation.
-    float, array_like
-        Current energies between pairs of particles in the simulation.
-    """
-    pairs = int(
-        (particles["xacceleration"].size - 1) * particles["xacceleration"].size / 2
-    )
-    distances = np.zeros(pairs)
-    energies = np.zeros(pairs)
-    distances, dx, dy = heavy.dist(
-        particles["xposition"], particles["yposition"], box_length
-    )
-    ff = forcefield(constants)
-    energies = ff.energy(distances)
-    energies[np.where(distances > cut_off)] = 0.0
-    return distances, energies
-
-
 def calculate_pressure(
-    particles, box_length, temperature, cut_off, constants, forcefield
+    particles, box_length, temperature, cut_off, constants, forcefield, mass
 ):
     r"""Calculates the instantaneous pressure of the simulation cell, found
     with the following relationship:
@@ -241,17 +220,16 @@ def calculate_pressure(
         the function forcefields.lennard_jones, theses are [A, B]
     forcefield: function (optional)
         The particular forcefield to be used to find the energy and forces.
+    mass: float (optional)
+        The mass of the particle being simulated (units of atomic mass units).
     Returns
     -------
     float:
         Instantaneous pressure of the simulation.
     """
-    distances, dx, dy = heavy.dist(
-        particles["xposition"], particles["yposition"], box_length
-    )
-    ff = forcefield(constants)
-    forces = ff.force(distances)
-    forces[np.where(distances > cut_off)] = 0.0
+    particles, distances, forces, energies = heavy.compute_force(
+        particles, box_length, cut_off, constants, forcefield, mass
+        )
     pres = np.sum(forces * distances)
     boltzmann_constant = 1.3806e-23  # joules / kelvin
     pres = 1.0 / (2 * box_length * box_length) * pres + (
@@ -290,7 +268,7 @@ def heat_bath(particles, temperature_sample, bath_temp):
 
 
 #Jit tag here had to be removed
-def dist(xposition, yposition, box_length):
+def dist(xposition, yposition, box_length, types):
     """Returns the distance array for the set of particles.
     Parameters
     ----------
@@ -302,18 +280,25 @@ def dist(xposition, yposition, box_length):
         y-dimension positions of the particles.
     box_length: float
         The box length of the simulation cell.
+    types: str, array_like (N)
+        Array of length N, where N is the number of particles, providing the
+        type of each particle.
     Returns
     -------
-    drr float, array_like ((N - 1) * N / 2))
+    drr: float, array_like ((N - 1) * N / 2))
         The pairs of distances between the particles.
-    dxr float, array_like ((N - 1) * N / 2))
+    dxr: float, array_like ((N - 1) * N / 2))
         The pairs of distances between the particles, in only the x-dimension.
-    dyr float, array_like ((N - 1) * N / 2))
+    dyr: float, array_like ((N - 1) * N / 2))
         The pairs of distances between the particles, in only the y-dimension.
+    pair_types: str, array_like ((N - 1) * N / 2))
+        The types of the two particles in each interaction, saved as 
+        'type1,type2' for each
     """
     drr = np.zeros(int((xposition.size - 1) * xposition.size / 2))
     dxr = np.zeros(int((xposition.size - 1) * xposition.size / 2))
     dyr = np.zeros(int((xposition.size - 1) * xposition.size / 2))
+    pair_types = []
     k = 0
     for i in range(0, xposition.size - 1):
         for j in range(i + 1, xposition.size):
@@ -325,8 +310,9 @@ def dist(xposition, yposition, box_length):
             drr[k] = dr
             dxr[k] = dx
             dyr[k] = dy
+            pair_types.append(types[i] + ',' + types[j])
             k += 1
-    return drr, dxr, dyr
+    return drr, dxr, dyr, pair_types
 
 
 #Jit tag here had to be removed
@@ -345,3 +331,28 @@ def pbc_correction(position, cell):
     if np.abs(position) > 0.5 * cell:
         position *= 1 - cell / np.abs(position)
     return position
+
+def create_dist_identifiers(type_identifier):
+    '''
+    Creates correct distance identifier matrix for particular type
+    of particle
+    Parameters
+    -------
+    type identifiers:
+        the identifier array listing 1 for particles of that type
+        or 0 for particles of a different type
+    Returns
+    -------
+    distances: float, array_like 
+        the distance identifier for interactions between each particle 
+        of that type, or 0 for interactions involving particles of a
+        different type
+
+    '''
+    distance_type_identifier = np.array([])
+    for index in range(len(type_identifier)):
+        if type_identifier[index]:
+            distance_type_identifier = np.append(distance_type_identifier,type_identifier[index+1:])
+        else:
+            distance_type_identifier = np.append(distance_type_identifier,np.zeros(len(type_identifier[index+1:])))
+    return distance_type_identifier
\ No newline at end of file
diff --git a/pylj/sample.py b/pylj/sample.py
index 3ad023e..b1de4b2 100644
--- a/pylj/sample.py
+++ b/pylj/sample.py
@@ -494,12 +494,13 @@ def setup_cellview(ax, system, scale=1):  # pragma: no cover
     scale: float (optional)
         The amount by which the particle size should be scaled down.
     """
-    xpos = system.particles["xposition"]
-    ypos = system.particles["yposition"]
     # These numbers are chosen to scale the size of the particles nicely to
     # allow the particles to appear to interact appropriately
     mk = 6.00555e-8 / (system.box_length - 2.2727e-10) - 1e-10 / scale
-    ax.plot(xpos, ypos, "o", markersize=mk, markeredgecolor="black", color="#34a5daff")
+    for index, type in enumerate(system.type_identifiers):
+        xpos = system.particles["xposition"] * type
+        ypos = system.particles["yposition"] * type   
+        ax.plot(xpos, ypos, "o", markersize=system.point_sizes[index], markeredgecolor="black")
     ax.set_xlim([0, system.box_length])
     ax.set_ylim([0, system.box_length])
     ax.set_xticks([])
@@ -612,11 +613,12 @@ def update_cellview(ax, system):  # pragma: no cover
     system: System
         The whole system information.
     """
-    x3 = system.particles["xposition"]
-    y3 = system.particles["yposition"]
-    line = ax.lines[0]
-    line.set_ydata(y3)
-    line.set_xdata(x3)
+    for index, type in enumerate(system.type_identifiers):
+        x3 = system.particles["xposition"] * type
+        y3 = system.particles["yposition"] * type
+        line = ax.lines[index]
+        line.set_ydata(y3)
+        line.set_xdata(x3)
 
 
 def update_rdfview(ax, system, average_rdf, r):  # pragma: no cover
diff --git a/pylj/tests/test_forcefields.py b/pylj/tests/test_forcefields.py
index a1376d8..5d57326 100644
--- a/pylj/tests/test_forcefields.py
+++ b/pylj/tests/test_forcefields.py
@@ -4,67 +4,69 @@
 
 
 class TestForcefields(unittest.TestCase):
-    def test_lennard_jones_energy(self):
+    def test_lennard_jones(self):
         a = forcefields.lennard_jones([1.0, 1.0])
         assert_almost_equal(a.energy(2.0), -0.015380859)
+        assert_almost_equal(a.force(2.0), -0.045410156)
+        a.mixing([2.0, 3.0])
+        assert_almost_equal([a.a, a.b], [1.4239324, 1.7379922])
         b = forcefields.lennard_jones([1.0, 1.0])
         assert_almost_equal(b.energy([2.0, 4.5],), [-0.015380859, -0.00012041])
+        assert_almost_equal(b.force([2.0, 4.0]), [-0.045410156, -3.66032124e-04])
         c = forcefields.lennard_jones([0.5, 3.5])
         assert_almost_equal(c.energy([2.0, 4.5],), [-0.05456543, -0.00042149])
+        assert_almost_equal(c.force([2.0, 4.5],), [-0.1633301, -0.000562])
         d = forcefields.lennard_jones( [5.0, 3.5])
         assert_almost_equal(d.energy([1.0, 1.5, 20.0]), [1.50000000, -0.268733500, -5.46874988e-08])
+        assert_almost_equal(d.force([1.0, 1.5, 20.0]), [ 3.9000000e+01, -9.2078707e-01, -1.6406249e-08])
         e = forcefields.lennard_jones([100.0, 300.0])
         assert_almost_equal(e.energy([100.0, 200.0, 500.0]), [0, 0, 0])
+        assert_almost_equal(e.force([100.0, 200.0, 500.0]), [0, 0, 0])
+        with self.assertRaises(IndexError):
+            f = forcefields.lennard_jones([1.0, 1.0, 1.0])
 
-    def test_lennard_jones_force(self):
-        a = forcefields.lennard_jones([1.0, 1.0])
-        assert_almost_equal(a.force(2.0), -0.045410156)
-        b = forcefields.lennard_jones([1.0, 1.0])
-        assert_almost_equal(b.force([2.0, 4.0, 6.0]), [-0.045410156, -3.66032124e-04, -2.14325517e-05])
-        c = forcefields.lennard_jones([1.5, 4.0])
-        assert_almost_equal(c.force([2.0, 4.0, 6.0]), [-0.185302734, -1.46457553e-03, -8.57325038e-05])
-        d = forcefields.lennard_jones([200.0, 500.0])
-        assert_almost_equal(d.force([150.0, 300.0, 500.0]), [-1.7558299e-12, -1.3717421e-14, -3.8400000e-16])
-
-    def test_lennard_jones_sigma_epsilon_energy(self):
+    def test_lennard_jones_sigma_epsilon(self):
         a = forcefields.lennard_jones_sigma_epsilon([1.0, 0.25])
         assert_almost_equal(a.energy(2.0), -0.015380859)
+        assert_almost_equal(a.force(2.0), -0.045410156)
         b = forcefields.lennard_jones_sigma_epsilon([1.0, 0.25])
         assert_almost_equal(b.energy([2.0, 1.0]), [-0.015380859, 0])
+        assert_almost_equal(b.force([2.0, 4.0]), [-0.0454102, -0.000366])
         c = forcefields.lennard_jones_sigma_epsilon([0.5, 0.75])
         assert_almost_equal(c.energy([2.0, 1.0, 1.5]), [-0.0007322, -0.0461425, -0.0041095])
+        assert_almost_equal(c.force([2.0, 1.0, 1.5]), [-0.0021962, -0.2724609, -0.0164157])
         d = forcefields.lennard_jones_sigma_epsilon([5e-10, 9e-9])
         assert_almost_equal(d.energy([400.0, 500.0, 600.0]), [0, 0, 0])
+        assert_almost_equal(d.force([400.0, 500.0, 600.0]), [0, 0, 0])
+        e = forcefields.lennard_jones_sigma_epsilon([1.0, 1.0])
+        e.mixing([4.0, 4.0])
+        assert_almost_equal([e.sigma, e.epsilon], [2.5, 2.0])
+        with self.assertRaises(IndexError):
+            f = forcefields.lennard_jones_sigma_epsilon([1.0, 1.0, 1.0])
 
-    def test_lennard_jones_sigma_epsilon_force(self):
-        a = forcefields.lennard_jones_sigma_epsilon([1.0, 0.25])
-        assert_almost_equal(a.force(2.0), -0.045410156)
-        b = forcefields.lennard_jones_sigma_epsilon([1.0, 0.25])
-        assert_almost_equal(b.force([2.0, 4.0]), [-0.0454102, -0.000366])
-        c = forcefields.lennard_jones_sigma_epsilon([3.0, 1.0])
-        assert_almost_equal(c.force([3.0, 4.0]), [8.0, -0.6877549])
-
-    def test_buckingham_energy(self):
+    def test_buckingham(self):
         a = forcefields.buckingham([1.0, 1.0, 1.0])
         assert_almost_equal(a.energy(2.0), 0.1197103832)
+        assert_almost_equal(a.force(2.0), 0.08846028324)
+        a.mixing([4.0, 4.0, 4.0])
+        assert_almost_equal([a.a, a.b, a.c], [2.0, 2.0, 2.0])
         b = forcefields.buckingham([1.0, 1.0, 1.0])
         assert_almost_equal(b.energy([2.0]), 0.1197103832)
+        assert_almost_equal(b.force([2.0]), 0.08846028324)
         c = forcefields.buckingham([1.0, 1.0, 1.0])
         assert_almost_equal(c.energy([2.0, 4.0]), [0.1197103832, 0.0180715])
+        assert_almost_equal(c.force([2.0, 4.0]), [0.0884603, 0.0179494])
         d = forcefields.buckingham([0.01, 0.01, 0.01])
         assert_almost_equal(d.energy([2.0, 4.0, 5.0]), [0.0096457, 0.0096055, 0.0095117])
+        assert_almost_equal(d.force([2.0, 4.0, 5.0]), [-3.7073013e-04,  9.2416835e-05,  9.4354942e-05])
+        with self.assertRaises(IndexError):
+            e = forcefields.buckingham([1.0, 1.0])
 
-    def test_buckingham_force(self):
-        a = forcefields.buckingham([1.0, 1.0, 1.0])
-        assert_almost_equal(a.force(2.0), 0.08846028324)
-        b = forcefields.buckingham([1.0, 1.0, 1.0])
-        assert_almost_equal(b.force([2.0]), 0.08846028324)
-        c = forcefields.buckingham([1.5, 0.1, 2.0])
-        assert_almost_equal(c.force([2.0, 1.0, 4.0]), [0.0290596, -11.8642744, 0.0998156])
-
-    def test_square_well_energy(self):
+    def test_square_well(self):
         a = forcefields.square_well([1.0, 1.5, 2.0])
         assert_equal(a.energy(2.0), -1.0)
+        with self.assertRaises(ValueError):
+            a.force()
         b = forcefields.square_well([1.0, 2.0, 1.25])
         assert_equal(b.energy(0.5), float('inf'))
         c = forcefields.square_well([0.5, 1.5, 1.25])
@@ -75,15 +77,8 @@ def test_square_well_energy(self):
         assert_equal(e.energy([3.0, 3.0, 0.25]), [0, 0, float('inf')])
         f = forcefields.square_well([1.0, 1.5, 1.25], max_val=5000)
         assert_equal(f.energy([3.0, 3.0, 0.25]), [0, 0, 5000])
-
-    def test_square_well_force(self):
-        a = forcefields.square_well([1.0, 1.5, 2.0])
-        with self.assertRaises(ValueError):
-            a.force()
-        b = forcefields.square_well([1.0, 1.5, 2.0])
-        with self.assertRaises(ValueError):
-            b.force()
-
+        with self.assertRaises(IndexError):
+            g = forcefields.square_well([1.0, 1.0])
 
 if __name__ == '__main__':
     unittest.main(exit=False)
diff --git a/pylj/tests/test_pairwise.py b/pylj/tests/test_pairwise.py
index 8efbe07..d1dc3e2 100644
--- a/pylj/tests/test_pairwise.py
+++ b/pylj/tests/test_pairwise.py
@@ -30,13 +30,14 @@ def test_compute_forces(self):
         particles = np.zeros(2, dtype=part_dt)
         particles["xposition"][0] = 1e-10
         particles["xposition"][1] = 5e-10
+        particles['types'] = ['0','0']
         particles, distances, forces, energies = pairwise.compute_force(
             particles,
             30,
             15,
-            constants=[1.363e-134, 9.273e-78],
+            constants=[[1.363e-134, 9.273e-78]],
             forcefield=ff.lennard_jones,
-            mass=39.948,
+            mass=39.948
         )
         assert_almost_equal(distances, [4e-10])
         assert_almost_equal(energies, [-1.4515047e-21])
@@ -45,38 +46,47 @@ def test_compute_forces(self):
         assert_almost_equal(particles["xacceleration"][0] / 1e14, 1.4451452)
         assert_almost_equal(particles["xacceleration"][1] / 1e14, -1.4451452)
 
-    def test_compute_energy(self):
-        part_dt = util.particle_dt()
-        particles = np.zeros(2, dtype=part_dt)
-        particles["xposition"][0] = 1e-10
-        particles["xposition"][1] = 5e-10
-        d, e = pairwise.compute_energy(
-            particles,
-            30,
-            15,
-            constants=[1.363e-134, 9.273e-78],
-            forcefield=ff.lennard_jones,
-        )
-        assert_almost_equal(d, [4e-10])
-        assert_almost_equal(e, [-1.4515047e-21])
-
     def test_calculate_pressure(self):
         part_dt = util.particle_dt()
         particles = np.zeros(2, dtype=part_dt)
         particles["xposition"][0] = 1e-10
         particles["xposition"][1] = 5e-10
+        particles['types'] = ['0','0']
         p = pairwise.calculate_pressure(
             particles,
             30,
             300,
             15,
-            constants=[1.363e-134, 9.273e-78],
+            constants=[[1.363e-134, 9.273e-78]],
             forcefield=ff.lennard_jones,
+            mass = 39.948
         )
-        assert_almost_equal(p * 1e24, 7.07368869)
+        assert_almost_equal(p * 1e24, 7.07368867)
 
     def test_pbc_correction(self):
         a = pairwise.pbc_correction(1, 10)
         assert_almost_equal(a, 1)
         b = pairwise.pbc_correction(11, 10)
         assert_almost_equal(b, 1)
+    
+    def test_multiple_particles(self):
+        part_dt = util.particle_dt()
+        particles = np.zeros(3, dtype=part_dt)
+        particles["xposition"][0] = 1e-10
+        particles["xposition"][1] = 5e-10
+        particles["yposition"][2] = 5e-10
+        particles['types'] = ['0','1','0']
+        particles, distances, forces, energies = pairwise.compute_force(
+            particles,
+            30,
+            15,
+            constants=[[1.363e-134, 9.273e-78],[1.363e-133, 9.273e-77]],
+            forcefield=ff.lennard_jones,
+            mass=39.948
+        )
+        assert_almost_equal(distances, [4.0000000e-10, 5.0990195e-10, 7.0710678e-10])
+        assert_almost_equal(energies, [-3.0626388e-20, -1.0201147e-20, -1.5468582e-21])
+        assert_almost_equal(forces, [-9.6342138e-11, -5.1698213e-12, -6.1773405e-12])
+        assert_almost_equal(particles["yacceleration"], [7.6421357e+13,  2.07196175e+13,-9.71409740e+13], decimal=-7)
+        assert_almost_equal(particles["xacceleration"][0] / 1e14, 4.4171075)
+        assert_almost_equal(particles["xacceleration"][1] / 1e14, -4.7771464)
\ No newline at end of file
diff --git a/pylj/tests/test_util.py b/pylj/tests/test_util.py
index df8c474..4c7c849 100644
--- a/pylj/tests/test_util.py
+++ b/pylj/tests/test_util.py
@@ -11,7 +11,7 @@ def test_system_square(self):
             300,
             8,
             mass=39.948,
-            constants=[1.363e-134, 9.273e-78],
+            constants=[[1.363e-134, 9.273e-78]],
             forcefield=ff.lennard_jones,
         )
         assert_equal(a.number_of_particles, 2)
@@ -34,7 +34,7 @@ def test_system_random(self):
             8,
             init_conf="random",
             mass=39.948,
-            constants=[1.363e-134, 9.273e-78],
+            constants=[[1.363e-134, 9.273e-78]],
             forcefield=ff.lennard_jones,
         )
         assert_equal(a.number_of_particles, 2)
@@ -57,7 +57,7 @@ def test_system_too_big(self):
                 300,
                 1000,
                 mass=39.948,
-                constants=[1.363e-134, 9.273e-78],
+                constants=[[1.363e-134, 9.273e-78]],
                 forcefield=ff.lennard_jones,
             )
         self.assertTrue(
@@ -73,7 +73,7 @@ def test_system_too_small(self):
                 300,
                 2,
                 mass=39.948,
-                constants=[1.363e-134, 9.273e-78],
+                constants=[[1.363e-134, 9.273e-78]],
                 forcefield=ff.lennard_jones,
             )
         self.assertTrue(
@@ -89,7 +89,7 @@ def test_system_init_conf(self):
                 100,
                 init_conf="horseradish",
                 mass=39.948,
-                constants=[1.363e-134, 9.273e-78],
+                constants=[[1.363e-134, 9.273e-78]],
                 forcefield=ff.lennard_jones,
             )
         self.assertTrue(
diff --git a/pylj/util.py b/pylj/util.py
index eda0122..bcf2cbc 100644
--- a/pylj/util.py
+++ b/pylj/util.py
@@ -45,13 +45,21 @@ def __init__(
         mass,
         init_conf="square",
         timestep_length=1e-14,
-        cut_off=15,
+        cut_off=15
     ):
         self.number_of_particles = number_of_particles
         self.init_temp = temperature
         self.constants = constants
-        self.forcefield = forcefield
         self.mass = mass
+        self.forcefield = forcefield
+        self.type_identifiers = None
+        self.particle_list = None
+        self.long_const = None
+        self.types = None
+        self.point_sizes = None
+        self.setup_point_sizes()
+        self.setup_type_identifiers()
+        self.setup_types()
         if box_length <= 600:
             self.box_length = box_length * 1e-10
         else:
@@ -115,6 +123,7 @@ def square(self):
         """
         part_dt = particle_dt()
         self.particles = np.zeros(self.number_of_particles, dtype=part_dt)
+        self.particles['types'] = self.types
         m = int(np.ceil(np.sqrt(self.number_of_particles)))
         d = self.box_length / m
         n = 0
@@ -130,10 +139,53 @@ def random(self):
         """
         part_dt = particle_dt()
         self.particles = np.zeros(self.number_of_particles, dtype=part_dt)
+        self.particles['types'] = self.types
         num_part = self.number_of_particles
         self.particles["xposition"] = np.random.uniform(0, self.box_length, num_part)
         self.particles["yposition"] = np.random.uniform(0, self.box_length, num_part)
 
+    def setup_types(self):
+        """Sets the long constants and types arrays of the particles
+        """
+        long_const = []
+        types = []
+        particle_list = []
+        for i in range(self.number_of_particles):
+            # Get set of constants and index
+            constants_type = self.constants[i%len(self.constants)]
+            particle_type = f'{i%len(self.constants)}'
+            # Append to lists
+            long_const.append(constants_type)
+            types.append(particle_type)
+            particle = Particle(constants_type, i, self.mass, particle_type)
+            particle.add(particle_list)
+        self.particle_list = particle_list
+        self.long_const = long_const
+        self.types = types
+
+    def setup_type_identifiers(self):
+        """Sets type-identifers arrays - legacy method now only used for plotting
+        """
+        # Creates arrays to identify which particle is in which type
+        number_of_types = len(self.constants)
+        type_identifiers = np.zeros((number_of_types,self.number_of_particles))
+        i = 0
+        while i < self.number_of_particles:
+            for k in range(number_of_types):
+                if i < self.number_of_particles:
+                    type_identifiers[k][i] = 1
+                    i+=1
+        self.type_identifiers = type_identifiers
+
+    def setup_point_sizes(self):
+        """Sets point sizes for use in plotting
+        """
+        point_sizes = []
+        for pair in self.constants:
+            size = self.forcefield(pair).point_size
+            point_sizes.append(size)
+        self.point_sizes = point_sizes
+
     def compute_force(self):
         """Maps to the compute_force function in either the comp (if Cython is
         installed) or the pairwise module and allows for a cleaner interface.
@@ -152,17 +204,9 @@ def compute_force(self):
         self.energies = energies
 
     def compute_energy(self):
-        """Maps to the compute_energy function in either the comp (if Cython
-        is installed) or the pairwise module and allows for a cleaner
-        interface.
+        """Maps to the compute_force function, as this also calculates energy
         """
-        self.distances, self.energies = md.compute_energy(
-            self.particles,
-            self.box_length,
-            self.cut_off,
-            self.constants,
-            self.forcefield,
-        )
+        self.compute_force()
 
     #Jit tag here had to be removed
     def integrate(self, method):
@@ -179,7 +223,7 @@ def integrate(self, method):
             self.cut_off,
             self.constants,
             self.forcefield,
-            self.mass,
+            self.mass
         )
 
     def md_sample(self):
@@ -234,6 +278,21 @@ def reject(self):
             self.position_store, self.particles, self.random_particle
         )
 
+class Particle:
+    def __init__(self,
+               constants,
+               index,
+               mass,
+               particle_type
+               ):
+        self.constants = constants
+        self.index = index
+        self.mass = mass
+        self.type = particle_type
+
+    def add(self, particles):
+        particles.append(self)
+        return particles
 
 def __cite__():  # pragma: no cover
     """This function will launch the website for the JOSE publication on
@@ -257,6 +316,7 @@ def particle_dt():
     - xprevious_position and yprevious_position
     - xforce and yforce
     - energy
+    - types
     """
     return np.dtype(
         [
@@ -271,5 +331,6 @@ def particle_dt():
             ("energy", np.float64),
             ("xpbccount", int),
             ("ypbccount", int),
+            ("types", list)
         ]
     )