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)