diff --git a/pylj/forcefields.py b/pylj/forcefields.py
index a644ea0..b4ea541 100644
--- a/pylj/forcefields.py
+++ b/pylj/forcefields.py
@@ -6,12 +6,6 @@ class lennard_jones(object):
     r"""Calculate the energy or force for a pair of particles using the
     Lennard-Jones (A/B variant) forcefield.
 
-    .. math::
-        E = \frac{A}{dr^{12}} - \frac{B}{dr^6}
-
-    .. math::
-        f = \frac{12A}{dr^{13}} - \frac{6B}{dr^7}
-
     Parameters
     ----------
     dr: float, array_like
@@ -19,13 +13,6 @@ class lennard_jones(object):
     constants: float, array_like
         An array of length two consisting of the A and B parameters for the
         12-6 Lennard-Jones function
-    force: bool (optional)
-        If true, the negative first derivative will be found.
-
-    Returns
-    -------
-    float: array_like
-        The potential energy or force between the particles.
     """
     def __init__(self, dr, constants):
         self.dr = dr
@@ -35,23 +22,41 @@ def __init__(self, dr, constants):
         self.force = self.update_force()
 
     def update_energy(self):
+        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}
+
+        Returns
+        -------
+        float: array_like
+        The potential energy between the particles.
+        """
         self.energy = self.a * np.power(self.dr, -12) - (self.b * np.power(self.dr, -6))
         return self.energy
     
     def update_force(self):
+        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}
+
+        Returns
+        -------
+        float: array_like
+        The force between the particles.
+        """
         self.force = 12 * self.a * np.power(self.dr, -13) - (6 * self.b * np.power(self.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.
 
-    .. math::
-        E = \frac{4e*a^{12}}{dr^{12}} - \frac{4e*a^{6}}{dr^6}
-
-    .. math::
-        f = \frac{48e*a^{12}}{dr^{13}} - \frac{24e*a^{6}}{dr^7}
-
     Parameters
     ----------
     dr: float, array_like
@@ -59,13 +64,6 @@ class lennard_jones_sigma_epsilon(object):
     constants: float, array_like
         An array of length two consisting of the sigma (a) and epsilon (e)
         parameters for the 12-6 Lennard-Jones function
-    force: bool (optional)
-        If true, the negative first derivative will be found.
-
-    Returns
-    -------
-    float: array_like
-        The potential energy or force between the particles.
     """
     def __init__(self, dr, constants):
         self.dr = dr
@@ -75,25 +73,43 @@ def __init__(self, dr, constants):
         self.force = self.update_force()
     
     def update_energy(self):
+        r"""Calculate the energy for a pair of particles using the
+        Lennard-Jones (sigma/epsilon variant) forcefield.
+
+        .. math::
+            E = \frac{4e*a^{12}}{dr^{12}} - \frac{4e*a^{6}}{dr^6}
+
+        Returns
+        -------
+        float: array_like
+        The potential energy between the particles.
+        """
         self.energy = 4 * self.epsilon * np.power(self.sigma, 12) * np.power(self.dr, -12) - (
                         4 * self.epsilon * np.power(self.sigma, 6) * np.power(self.dr, -6)) 
         return self.energy 
     
     def update_force(self):
+        r"""Calculate the force for a pair of particles using the
+        Lennard-Jones (sigma/epsilon variant) forcefield.
+
+        .. math::
+            f = \frac{48e*a^{12}}{dr^{13}} - \frac{24e*a^{6}}{dr^7}
+
+        Returns
+        -------
+        float: array_like
+        The force between the particles.
+        """
         self.force = 48 * self.epsilon * np.power(self.sigma, 12) * np.power(
             self.dr, -13) - (24 * self.epsilon * np.power(self.sigma, 6) * np.power(self.dr, -7))
         return self.force
 
+
+
 class buckingham(object):
     r""" Calculate the energy or force for a pair of particles using the
     Buckingham forcefield.
 
-    .. math::
-        E = Ae^{(-Bdr)} - \frac{C}{dr^6}
-
-    .. math::
-        f = ABe^{(-Bdr)} - \frac{6C}{dr^7}
-
     Parameters
     ----------
     dr: float, array_like
@@ -101,13 +117,6 @@ class buckingham(object):
     constants: float, array_like
         An array of length three consisting of the A, B and C parameters for
         the Buckingham function.
-    force: bool (optional)
-        If true, the negative first derivative will be found.
-
-    Returns
-    -------
-    float: array_like
-        The potential energy or force between the particles.
     """
     def __init__(self, dr, constants):
         self.dr = dr
@@ -118,33 +127,41 @@ def __init__(self, dr, constants):
         self.force = self.update_force()
     
     def update_energy(self):
+        r"""Calculate the energy for a pair of particles using the
+        Buckingham forcefield.
+
+        .. math::
+            E = Ae^{(-Bdr)} - \frac{C}{dr^6}
+
+        Returns
+        -------
+        float: array_like
+        The potential energy between the particles.
+        """
         self.energy = self.a * np.exp(- np.multiply(self.b, self.dr)) - self.c / np.power(self.dr, 6)
         return self.energy
     
     def update_force(self):
+        r"""Calculate the force for a pair of particles using the
+        Buckingham forcefield.
+
+        .. math::
+            f = ABe^{(-Bdr)} - \frac{6C}{dr^7}
+
+        Returns
+        -------
+        float: array_like
+        The force between the particles.
+        """
         self.force = self.a * self.b * np.exp(- np.multiply(self.b, self.dr)) - 6 * self.c / np.power(self.dr, 7)
         return self.force
 
+
+
 class square_well(object):
     r'''Calculate the energy or force for a pair of particles using a
     square well model.
 
-    .. math::
-        E = {
-        if dr < sigma:
-            E = max_val
-        elif sigma <= dr < lambda * sigma:
-            E = -epsilon
-        elif r >= lambda * sigma:
-            E = 0
-        }
-    .. math::
-        f = {
-        if sigma <= dr < lambda * sigma:
-            f = inf
-        else:
-            f = 0
-        }
     Parameters
     ----------
     dr: float, array_like
@@ -154,13 +171,6 @@ class square_well(object):
         parameters for the square well model.
     max_val: int (optional)
         Upper bound for values in square well - replaces usual infinite values
-    force: bool (optional)
-        If true, the negative first derivative will be found.
-
-    Returns
-    -------
-    float: array_like
-        The potential energy between the particles.
     '''
     def __init__(self, dr, constants, max_val=np.inf):
 
@@ -179,22 +189,52 @@ def __init__(self, dr, constants, max_val=np.inf):
         self.force = 0
 
     def update_energy(self):
-            E = np.zeros_like(self.dr)
-            E = np.zeros_like(self.dr)
-            E[np.where(self.dr < self.epsilon)] = self.max_val
-            E[np.where(self.dr >= self.lamda * self.sigma)] = 0
-
-            # apply mask for sigma <= self.dr < lambda * sigma
-            a = self.sigma <= self.dr
-            b = self.dr < self.lamda * self.sigma
-            E[np.where(a & b)] = -self.epsilon
-
-            if len(E) == 1:
-                self.energy = float(E[0])
-            else:
-                self.energy = np.array(E, dtype='float')
+        r'''Calculate the energy for a pair of particles using a
+        square well model.
+
+        .. math::
+            E = {
+            if dr < sigma:
+                E = max_val
+            elif sigma <= dr < lambda * sigma:
+                E = -epsilon
+            elif r >= lambda * sigma:
+                E = 0
+            }
+
+        Returns
+        -------
+        float: array_like
+            The potential energy between the particles.
+        '''
+        E = np.zeros_like(self.dr)
+        E = np.zeros_like(self.dr)
+        E[np.where(self.dr < self.epsilon)] = self.max_val
+        E[np.where(self.dr >= self.lamda * self.sigma)] = 0
+
+        # apply mask for sigma <= self.dr < lambda * sigma
+        a = self.sigma <= self.dr
+        b = self.dr < self.lamda * self.sigma
+        E[np.where(a & b)] = -self.epsilon
+
+        if len(E) == 1:
+            self.energy = float(E[0])
+        else:
+            self.energy = np.array(E, dtype='float')
             
-            return self.energy
+        return self.energy
         
     def update_force(self):
-            raise ValueError("Force is infinite at sigma <= dr < lambda * sigma")
\ No newline at end of file
+        r'''The force of a pair of particles using a square well model is given by:
+
+        .. math::
+        f = {
+        if sigma <= dr < lambda * sigma:
+            f = inf
+        else:
+            f = 0
+        }
+
+        Therefore the force here will always be infinite, and therefore not possible to simulate
+        '''
+        raise ValueError("Force is infinite at sigma <= dr < lambda * sigma")
\ No newline at end of file