Takes in a list of SlabEntries and calculates the chemical potentials
            at which all slabs in the list coexists simultaneously. Useful for
            building surface phase diagrams. Note that to solve for x equations
            (x slab_entries), there must be x free variables (chemical potentials).
            Adjust delu_dict as need be to get the correct number of free variables.
        Args:
            slab_entries (array): The coefficients of the first equation
            delu_dict (Dict): Dictionary of the chemical potentials to be set as
                constant. Note the key should be a sympy Symbol object of the
                format: Symbol("delu_el") where el is the name of the element.

        Returns:
            (array): Array containing a solution to x equations with x
                variables (x-1 chemical potential and 1 surface energy)

def get_surface_equilibrium(self, slab_entries, delu_dict=None):

        """
        Takes in a list of SlabEntries and calculates the chemical potentials
            at which all slabs in the list coexists simultaneously. Useful for
            building surface phase diagrams. Note that to solve for x equations
            (x slab_entries), there must be x free variables (chemical potentials).
            Adjust delu_dict as need be to get the correct number of free variables.
        Args:
            slab_entries (array): The coefficients of the first equation
            delu_dict (Dict): Dictionary of the chemical potentials to be set as
                constant. Note the key should be a sympy Symbol object of the
                format: Symbol("delu_el") where el is the name of the element.

        Returns:
            (array): Array containing a solution to x equations with x
                variables (x-1 chemical potential and 1 surface energy)
        """
        # Generate all possible coefficients
        all_parameters = []
        all_eqns = []
        for slab_entry in slab_entries:
            se = self.surfe_dict[slab_entry]

            # remove the free chempots we wish to keep constant and
            # set the equation to 0 (subtract gamma from both sides)
            if type(se).__name__ == "float":
                all_eqns.append(se - Symbol("gamma"))
            else:

                se = sub_chempots(se, delu_dict) if delu_dict else se
                all_eqns.append(se - Symbol("gamma"))
                all_parameters.extend([p for p in list(se.free_symbols)
                                       if p not in all_parameters])

        all_parameters.append(Symbol("gamma"))
        # Now solve the system of linear eqns to find the chempot
        # where the slabs are at equilibrium with each other

        soln = linsolve(all_eqns, all_parameters)
        if not soln:
            warnings.warn("No solution")
            return soln
        return {p: list(soln)[0][i] for i, p in enumerate(all_parameters)}