Returns a matplotlib object with the plotted reaction path.

    Parameters
    ----------
    df : Pandas DataFrame generated by reaction_network
    temperature : numeric
        temperature in K
    pressure : numeric
       pressure in mbar
    pH : PH in bulk solution
    potential : Electric potential vs. SHE in eV
    e_lim: Limits for the energy axis.

    Returns
    -------
    fig: matplotlib object.

def plot_reaction_scheme(df, temperature, pressure, potential, pH, e_lim=None):
    """Returns a matplotlib object with the plotted reaction path.

    Parameters
    ----------
    df : Pandas DataFrame generated by reaction_network
    temperature : numeric
        temperature in K
    pressure : numeric
       pressure in mbar
    pH : PH in bulk solution
    potential : Electric potential vs. SHE in eV
    e_lim: Limits for the energy axis.

    Returns
    -------
    fig: matplotlib object.
    """
    ncols = int((df.shape[0]/20)) +1
    fig_width = ncols + 1.5*len(df['intermediate_labels'][0])
    figsize = (fig_width, 6)
    fig, ax = plt.subplots(figsize=figsize)

    if pressure == None:
        pressure_label = '0'
    else:
        pressure_label = str(pressure)

    lines = []
    for j, energy_list in enumerate(df['reaction_energy']):
        ts = df['transition_states'][j]
        R = df['reaction_coordinate'][j]
        E = [[x, x] for x in energy_list]
        labels = df['system_label']



        for i, n in enumerate(R):
            if i == 0:
                line = Line2D([0], [0], color=colors[j], lw=4)
                lines.append(line)
                ax.plot(n, E[i], ls='-', color=colors[j], linewidth=3.25, solid_capstyle='round',
                        path_effects=[pe.Stroke(linewidth=6, foreground=edge_colors[j]), pe.Normal()],
                        label=labels[j])
                ax.plot([n[1], n[1] + 0.5], [E[i], E[i + 1]], ls='--',
                        dashes=(3, 2), color=colors[j], linewidth=1.)
            else:
                if ts[i]:
                    xts = [R[i-1][1], R[i][0], R[i+1][0]]
                    yts = [energy_list[i-1], energy_list[i], energy_list[i+1]]
                    z1 = np.polyfit(xts, yts, 2)
                    xp1 = np.linspace(xts[0], xts[2], 100)
                    p1 = np.poly1d(z1)
                    ax.plot(xp1, p1(xp1), ls='--', color=colors[j], linewidth=2.)
                    ax.plot(xts[1], yts[1], marker = 'o', c=colors[j], mec = edge_colors[j], lw=1.5,
                            markersize=7)

                else:
                    ax.plot(n, E[i], ls='-', color=colors[j], linewidth=3.25, solid_capstyle='round',
                            path_effects=[pe.Stroke(linewidth=6, foreground=edge_colors[j]), pe.Normal()])
                    if i < len(R) - 1:
                        ax.plot([n[1], n[1] + 0.5], [E[i], E[i + 1]], ls='--',
                                dashes=(3, 2), color=colors[j], linewidth=1.)

    ax.legend(handlelength=0.4, ncol=ncols, loc=2, frameon=False,
              bbox_to_anchor=(1.05, 1), borderaxespad=0., fontsize=12)
    if e_lim:
        ax.set_ylim(e_lim)
    ax.set_xlabel('Reaction coordinate')
    ax.set_ylabel('Reaction free energy (eV)')
    reaction_labels = df['intermediate_labels'][0]
    reaction_labels = [sub(w) for w in reaction_labels]
    plt.xticks(np.arange(len(reaction_labels)) + 0.25, tuple(reaction_labels), rotation=45)
    # plt.tight_layout()

    a = ax.get_xlim()[1]+0.05*ax.get_xlim()[1]
    b = ax.get_ylim()[0]+0.05*ax.get_ylim()[1]

    if potential is not None and pH is not None:
        ax.text(a,b, 'U = '+str(potential)+' eV vs. SHE \n pH = '
                     +str(pH)+' \n T = '+str(temperature)
                     +' K \n p = '+pressure_label+' mbar',fontsize=12)
    else:
        ax.text(a,b, 'T = '+str(temperature)+' \n p = '+pressure_label+' mbar',fontsize=12)
    plt.tight_layout()
    return(fig)