def calculate_excess(adsorbant, slab_cations, area, bulk, nspecies=1, check
    =False):
    """Calculates the excess of a given species at the surface.
    Depending on the nature of the species, there are two ways to do this.
    If the species is a constituent part of the surface, e.g.
    Oxygen in :math:`TiO_2` then the calculation must account for
    the stoichiometry of that material. Using the :math:`TiO_2` example

    .. math::
        \\Gamma_O = \\frac{1}{2A} \\Bigg( nO_{Slab} - \\frac{nO_{Bulk}}
        {nTi_{Bulk}}nTi_{Slab}  \\Bigg)

    where :math:`nO_{Slab}` is the number of oxygen in the slab,
    :math:`nO_{Bulk}` is the number of oxygen in the bulk,
    A is the surface area, :math:`nTi_{Bulk}` is the number of Ti in
    the bulk and :math:`nTi_{Slab}` is the number of Ti in the slab. 
    If the species is just an external adsorbant, e.g. water or carbon dioxide
    then one does not need to consider the state of the surface,
    as there was none there to begin with.

    .. math::
        \\Gamma_{H_2O} = \\frac{nH_2O}{2A}

    where :math:`nH_2O` is the number of water molecules and A is the
    surface area.

    Parameters
    ----------
    adsorbant : int
        Number of species
    slab_cations : int
        Number of cations
    area : float
        Area of surface
    bulk : dic
        Dictonary of bulk properties
    nspecies : int (optional)
        number of external species
    check : bool (optional)
        Check if this is an external or constituent species.

    Returns
    -------
    float:
        Surface excess of given species.
    """
    if check is True and nspecies == 1:
        return (adsorbant - bulk['Anion'] / bulk['Cation'] * slab_cations) / (
            2 * area)
    else:
        return adsorbant / (area * 2)
