Match adsorbate structures with reaction entries

def read_energies(self, root):
        self.key_value_pairs_reaction = None

        slab_structures = collect_structures(root)

        if len(slab_structures) == 0:
            self.raise_warning('No structure files in {root}: Skipping this folder'
                               .format(root=root))
            return

        # Remove old species from ase_ids
        all_reaction_species = list(self.reaction.values())[0] \
            + list(self.reaction.values())[1]
        all_reaction_species = [clear_prefactor(rs) for rs in
                                all_reaction_species]

        for ase_id in list(self.ase_ids.keys()):
            if ase_id == 'star' or 'bulk' in ase_id:
                continue
            if not ase_id in all_reaction_species:
                del self.ase_ids[ase_id]

        if 'TS' in self.structures:  # Delete old TS
            del self.structures['TS']
            del self.structures['TSempty']
            del self.prefactors['TS']

        for k in list(self.structures.keys()):
            if 'neb' in k:
                del self.structures[k]


        neb_indices = [i for i, slab in enumerate(slab_structures) if 'neb' in
                       slab[-1].info['filename']]
        neb_names = {}

        if len(neb_indices) == 1:
            index = neb_indices[0]
            slab = slab_structures[index]
            f = slab[-1].info['filename']
            del slab_structures[index]
            neb_indices = []
            for i, s in enumerate(slab):
                s.info['filename'] = f
                slab_structures.append(s)
                index = len(slab_structures) - 1
                neb_indices += [index]
                neb_names.update({str(i): 'neb' + str(i)})

        elif len(neb_indices) > 1:
            for i in neb_indices:
                f = slab_structures[i][-1].info['filename']
                neb_names.update({str(i): os.path.basename(f).split('.')[0]})

        for i, slab in enumerate(slab_structures):
            if isinstance(slab, list):
                slab_structures[i] = slab[-1]

        empty = self.empty

        reactant_entries = self.reaction['reactants'] + \
            self.reaction['products']
        if not empty:
            if 'star' in reactant_entries and len(neb_indices) == 0:
                message = 'Empty slab needed for reaction!'
                self.raise_error(message)
                return
            else:
                empty0 = slab_structures[0]
                chemical_composition = ase_tools.get_chemical_formula(empty0)
                self.raise_warning("Using '{}' as a reference instead of empty slab"
                                   .format(empty0.info['filename']))
                empty_atn = list(empty0.get_atomic_numbers())
        else:
            chemical_composition = ase_tools.get_chemical_formula(empty)
            empty_atn = list(empty.get_atomic_numbers())

        prefactor_scale = copy.deepcopy(self.prefactors)
        for key1, values in prefactor_scale.items():
            prefactor_scale[key1] = [1 for v in values]

        key_value_pairs = {}

        key_value_pairs.update({'name':
                                chemical_composition,
                                'facet': self.ase_facet,
                                'state': 'star'})

        """ Match adsorbate structures with reaction entries"""
        for i, slab in enumerate(slab_structures):
            f = slab.info['filename']
            atns = list(slab.get_atomic_numbers())
            if not (np.array(atns) > 8).any() and \
               (np.array(empty_atn) > 8).any():
                self.raise_warning("Only molecular species for structure: {}"
                                   .format(f))
                continue

            """Get supercell size relative to empty slab"""
            supercell_factor = 1
            if len(atns) > len(empty_atn) * 2:  # different supercells
                supercell_factor = len(res_slab_atn) // len(empty_atn)

            """Atomic numbers of adsorbate"""
            ads_atn = []
            if 'star' in reactant_entries and len(neb_indices) == 0:
                ads_atn = copy.copy(atns)
                for atn in empty_atn * supercell_factor:
                    try:
                        ads_atn.remove(atn)
                    except ValueError as e:
                        self.raise_error(
                            'Empty slab: {} contains species not in: {}' \
                                .format(empty.info['filename'], f))
                ads_atn = sorted(ads_atn)
                if ads_atn == []:
                    self.raise_warning("No adsorbates for structure: {}"
                                   .format(f))
                    continue

            ase_id = None
            id, ase_id = ase_tools.check_in_ase(slab, self.cathub_db)
            key_value_pairs.update({'epot': ase_tools.get_energies([slab])})

            if 'empty' in f and 'TS' in f:  # empty slab for transition state
                self.structures.update({'TSempty': [slab]})
                self.prefactors.update({'TSempty': [1]})
                prefactor_scale.update({'TSempty': [1]})
                key_value_pairs.update({'species': ''})
                if ase_id is None:
                    ase_id = ase_tools.write_ase(slab, self.cathub_db,
                                                 self.stdout,
                                                 self.user,
                                                 **key_value_pairs)
                elif self.update:
                    ase_tools.update_ase(self.cathub_db, id, self.stdout,
                                         **key_value_pairs)
                self.ase_ids.update({'TSemptystar': ase_id})
                continue

            elif 'TS' in f:  # transition state
                self.structures.update({'TS': [slab]})
                self.prefactors.update({'TS': [1]})
                prefactor_scale.update({'TS': [1]})
                key_value_pairs.update({'species': 'TS'})
                if ase_id is None:
                    ase_id = ase_tools.write_ase(slab, self.cathub_db,
                                                 self.stdout,
                                                 self.user,
                                                 **key_value_pairs)
                elif self.update:
                    ase_tools.update_ase(self.cathub_db, id, self.stdout,
                                         **key_value_pairs)
                self.ase_ids.update({'TSstar': ase_id})
                continue

            if i in neb_indices:
                self.structures.update({neb_names[str(i)]: [slab]})
                key_value_pairs.update({'species': 'neb'})
                if ase_id is None:
                    ase_id = ase_tools.write_ase(slab,
                                                 self.cathub_db,
                                                 self.stdout,
                                                 self.user,
                                                 **key_value_pairs)
                elif self.update:
                    ase_tools.update_ase(self.cathub_db, id, self.stdout,
                                         **key_value_pairs)
                self.ase_ids.update({neb_names[str(i)]: ase_id})
                continue
            

            found = False
            for key, mollist in self.reaction.items():
                if found:
                    break
                for n, molecule in enumerate(mollist):
                    if found:
                        break
                    if not self.states[key][n] == 'star':  # only slab stuff
                        continue
                    if not self.structures[key][n] == '':  # allready found
                        continue
                    molecule = clear_state(clear_prefactor(molecule))
                    if molecule == '':
                        continue
                    molecule_atn = sorted(
                        ase_tools.get_numbers_from_formula(molecule))
                    if slab.info['filename'] == molecule:
                        if ads_atn == molecule_atn:
                            found = True
                            match_key = key
                            match_n = n
                            break
                        else:
                            raise_warning('Name of file does not match chemimcal formula: {}'
                                          .format(slab.info['filename']))

                    for n_ads in range(1, 5):
                        mol_atn = sorted(molecule_atn * n_ads)
                        if (ads_atn == mol_atn or len(ads_atn) == 0):
                            match_n_ads = n_ads
                            found = True
                            match_key = key
                            match_n = n
                            break

            if found:
                key = match_key
                n = match_n
                n_ads = match_n_ads
                self.structures[key][n] = slab
                species = clear_prefactor(
                    self.reaction[key][n])
                id, ase_id = ase_tools.check_in_ase(
                    slab, self.cathub_db)
                key_value_pairs.update(
                    {'species':
                     clear_state(
                         species),
                     'n': n_ads,
                     'site': str(self.sites.get(species, ''))})
                if ase_id is None:
                    ase_id = ase_tools.write_ase(
                        slab, self.cathub_db, self.stdout,
                        self.user,
                        **key_value_pairs)
                elif self.update:
                    ase_tools.update_ase(
                        self.cathub_db, id, self.stdout,
                        **key_value_pairs)
                self.ase_ids.update({species: ase_id})

            if n_ads > 1:
                for key1, values in prefactor_scale.items():
                    for mol_i in range(len(values)):
                        if self.states[key1][mol_i] == 'gas':
                            prefactor_scale[key1][mol_i] = n_ads

            if supercell_factor > 1:
                for key2, values in prefactor_scale.items():
                    for mol_i in range(len(values)):
                        if self.reaction[key2][mol_i] == 'star':
                            prefactor_scale[key2][mol_i] *= supercell_factor


        # Check that all structures have been found
        structurenames = [s for s in list(self.structures.keys())
                          if s not in ['reactants', 'products']]
        for k in ['reactants', 'products']:
            structurenames += [s for s in self.structures[k] if s != ''
                               and s is not None]
        only_neb = np.all(['neb' in s for s in structurenames])
        surface_composition = self.metal

        if only_neb:
            if not self.empty:
                for ads in self.reaction_atoms['reactants']:
                    ads_atn = ase_tools.get_numbers_from_formula(ads)
                    for atn in ads_atn:
                        empty_atn.remove(atn)
                chemical_composition = \
                    ase_tools.get_formula_from_numbers(empty_atn, mode='metal')
            neb_numbers = []
            neb_energies = []
            for key, structure in self.structures.items():
                if key in ['reactants', 'products']:
                    continue

                neb_no = int(key.split('.')[0].replace('neb', ''))
                neb_numbers += [neb_no]
                neb_energies += [structure[0].get_potential_energy()]

            initial = neb_energies[np.argmin(neb_numbers)]
            final = neb_energies[np.argmax(neb_numbers)]
            TS = np.max(neb_energies)

            reaction_energy = final - initial
            activation_energy = TS - initial

            if activation_energy == 0 or activation_energy == reaction_energy:
                activation_energy = None
        else:
            for key, structurelist in self.structures.items():
                if '' in structurelist:
                    index = structurelist.index('')
                    molecule = clear_state(
                        clear_prefactor(self.reaction[key][index]))
                    if self.states[key][index] == 'star':
                        message = "Adsorbate '{}' not found for any structure files in '{}'."\
                            .format(molecule, root) + \
                            " Please check your adsorbate structures and the empty slab."
                    if self.states[key][index] == 'gas':
                        message = "Gas phase molecule '{}' not found for any structure files in '{}'."\
                            .format(molecule, self.gas_folder) + \
                            " Please check your gas phase references."
                    self.raise_error(message)
                    return

            prefactors_final = copy.deepcopy(self.prefactors)
            for key in self.prefactors:
                for i, v in enumerate(self.prefactors[key]):
                    prefactors_final[key][i] = self.prefactors[key][i] * \
                        prefactor_scale[key][i]

            reaction_energy = None
            activation_energy = None
            try:
                reaction_energy, activation_energy = \
                    ase_tools.get_reaction_energy(
                        self.structures, self.reaction,
                        self.reaction_atoms,
                        self.states, prefactors_final,
                        self.prefactors_TS,
                        self.energy_corrections)
            except BaseException as e:
                message = "reaction energy failed for files in '{}'"\
                    .format(root)
                self.raise_error(message + '\n' + str(e))

        if not -self.energy_limit < reaction_energy < self.energy_limit:
            self.raise_error('reaction energy is very large: {} eV \n  '
                             .format(reaction_energy) +
                             'Folder: {}. \n  '.format(root) +
                             'If the value is correct, you can reset the limit with cathub folder2db --energy-limit <value>. Default is --energy-limit=5 (eV)'
                             )
        if activation_energy is not None:
            if activation_energy < reaction_energy:
                self.raise_warning('activation energy is smaller than reaction energy: {} vs {} eV \n  Folder: {}'.format(
                    activation_energy, reaction_energy, root))
            if not activation_energy < self.energy_limit:
                self.raise_error(' Very large activation energy: {} eV \n  Folder: {}'
                                 .format(activation_energy, root))

        reaction_info = {'reactants': {},
                         'products': {}}

        for key in ['reactants', 'products']:
            for i, r in enumerate(self.reaction[key]):
                r = clear_prefactor(r)
                reaction_info[key].update({r: self.prefactors[key][i]})

        self.key_value_pairs_reaction = {
            'chemical_composition': chemical_composition,
            'surface_composition': surface_composition,
            'facet': self.facet,
            'sites': self.sites,
            'coverages': self.coverages,
            'reactants': reaction_info['reactants'],
            'products': reaction_info['products'],
            'reaction_energy': float(reaction_energy),
            'activation_energy': activation_energy,
            'dft_code': self.DFT_code,
            'dft_functional': self.DFT_functional,
            'pub_id': self.pub_id,
            'doi': self.doi,
            'year': int(self.year),
            'ase_ids': self.ase_ids,
            'energy_corrections': self.energy_corrections,
            'username': self.user}