docs: def reaction_scheme(self, df, temperature, pressure, pH, potential):
          """Returns a dataframe with Gibbs free reaction energies.
  
          Parameters
          ----------
          df : Pandas DataFrame generated by db_to_df
          temperature : numeric
              temperature in K
          pressure : numeric
             pressure in mbar
          pH : PH in bulk solution
          potential : Electric potential vs. SHE in eV
  
          Returns
          -------
          df : DataFrame suitable for plotting.
          """
          # set reaction scheme
          reactions = self.intermediates
          df_param = self.df_reaction_parameters.fillna(0)
  
          # set reaction labels for plotting
          reaction_labels = [w.replace('gas', '$_g$') for w in reactions]
          reaction_labels = [w.replace('star', '*') for w in reaction_labels]
          reaction_scheme = reactions[1:]
          reaction_number = np.arange(len(reactions))
  
          # get unique substrates
          labels = list(df.labels)
  
          # Gather and summarize data from data base.
          unique_system_label = []
          surface = []
          facet = []
          intermediate_label = []
          reaction_coordinate = []
          energies = []
          transition_states = []
  
          # compute proton-electron energetics
          G_PE = 0.0
          if potential is not None and pH is not None:
              ne = 1.0
              G_PE = -ne * potential + \
                     proton_hydroxide_free_energy(temperature=temperature, pressure=pressure, pH=pH)[0]
  
          for lab in list(set(labels)):
              energies_system = [0.0]
              reaction_coordinate_system = [[x, x + 0.5] for x in reaction_number]
              reaction_coordinate.append(reaction_coordinate_system)
              for step, products in enumerate(reaction_scheme):
                  df_tmp = df[df['labels'].apply(lambda x: lab in x)]
                  df_tmp = df_tmp[df_tmp['products'].apply(lambda x: products in x)]
  
                  # Set initial reactants as reference energy = 0 eV.
                  if step == 0:
                      unique_system_label.append(df_tmp.iloc[0]['labels'])
                      surface.append(df_tmp.iloc[0]['system'])
                      facet.append(df_tmp.iloc[0]['facet'])
                      intermediate_label.append(reaction_labels)
                      previous_reaction_energy = 0
                      ts = [False]
  
                  # Get free energy corrections (FEC).
                  all_reactants = list(df_tmp['reactants'])[0]
                  all_products = list(df_tmp['products'])[0]
                  reactant_FEC = get_FEC(all_reactants, temperature, pressure)
                  product_FEC = get_FEC(all_products, temperature, pressure)
                  FEC = product_FEC - reactant_FEC
  
                  # Get betas, and corrections.
                  intermediate_name = reactions[step + 1]
                  beta = df_param[df_param['intermediate'] == intermediate_name].betas.values[0]
                  net_correction = df_param[df_param['intermediate'] == intermediate_name].net_corrections.values[0]
                  ts.append(df_param[df_param['intermediate'] == intermediate_name].transition_state.values[0])
  
                  # Calculate the energy with respect to the initial reactants.
                  reaction_energy = previous_reaction_energy + df_tmp.iloc[0]['reaction_energy'] \
                                    + FEC - beta * G_PE + net_correction
  
                  # print('FEC : '+str(FEC))
                  previous_reaction_energy = reaction_energy
                  energies_system.append(reaction_energy)
              energies.append(energies_system)
              transition_states.append(ts)
  
  
          # Dump data into data frame.
          df_new = pd.DataFrame()
          df_new['facet'] = facet
          df_new['intermediate_labels'] = intermediate_label
          df_new['reaction_coordinate'] = reaction_coordinate
          df_new['reaction_energy'] = energies
          df_new['system_label'] = unique_system_label
          df_new['system'] = surface
          df_new['transition_states'] = transition_states
          df_new = df_new.sort_values(by=['facet', 'system'])
          df_new = df_new.reset_index(drop=True)
          self.df_react = df_new
          return(True)