docs: 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)