#!/usr/bin/env python3
# -*- coding: utf-8 -*-
__copyright__ = """ This code is licensed under the 3-clause BSD license.
Copyright ETH Zurich, Laboratory of Physical Chemistry, Reiher Group.
See LICENSE.txt for details.
"""

# Standard library imports
import time
from collections import defaultdict
from itertools import combinations, product
from typing import Dict, Iterator, List, Tuple, Optional, Generator
from json import dumps

# Third party imports
from numpy import ndarray
import scine_database as db
import scine_utilities as utils

# Local application imports
from scine_chemoton.utilities.queries import model_query, calculation_exists_in_structure, get_calculation_id
from .connectivity_analyzer import ReactionType, ConnectivityAnalyzer
from . import TrialGenerator, _sanity_check_wrapper


class FragmentBased(TrialGenerator):
    """
    Class to generate reactive complex calculations via fragment-based approaches.

    Attributes
    ----------
    options : FragmentBased.Options
        The options for generating reactive complex calculations.
    reactive_site_filter : ReactiveSiteFilter
        The filter applied to determine reactive sites.
        NOTE: Reactive pairs/coordinates for the application of the
        `filter_atom_pairs` and `filter_reaction_coordinates` methods are
        generated by probing all possible across-fragment combinations for
        unimolecular reactions (i.e. each pair has to consist of atoms belonging
        to two different fragments). If any of the combinations belonging to
        a certain fragment composition is allowed, the fragment composition
        is considered.
        For bimolecular reactions with two diatomic fragments the reactive pairs
        correspond to the atoms facing each other in the (unrotated) reactive
        complex.
    """

    class Options(TrialGenerator.Options):
        """
        The options for the fragment-based reactive complex enumeration.
        """

        __slots__ = (
            "unimolecular_dissociation_options",
            "unimolecular_association_options",
            "bimolecular_association_options",
        )

        class BimolAssociationOptions:
            """
            The options for the exploration of bimolecular reactions.
            """

            __slots__ = (
                "job",
                "job_settings",
                "complex_generator",
                "consider_diatomic_fragments",
                "max_within_fragment_graph_distance",
                "minimal_spin_multiplicity",
            )

            def __init__(self):
                self.job: db.Job = db.Job("scine_react_complex_nt")
                """
                db.Job (Scine::Database::Calculation::Job)
                    The Job used to do reactive complex calculations/trial
                    elementary step calculations.
                    `scine_react_complex_nt` and `scine_react_complex_afir`
                    jobs are supported.
                    The default is: the `scine_react_complex_nt` order on a single core.
                """
                self.job_settings: utils.ValueCollection = utils.ValueCollection({})
                """
                Additional settings added to the elementary step trial calculation.
                    Empty by default.
                """
                self.consider_diatomic_fragments = False
                """
                bool
                    Whether to also generate reactive complexes with diatomic
                    reactive fragments (i.e., reactive complex calculations with
                    the lhs and/or rhs lists containing two atoms). (default: False)
                """
                self.max_within_fragment_graph_distance = 1
                """
                int
                    Maximum number of bonds between two atoms forming a reactive
                    fragment. This option is only effective if
                    `consider_diatomic_fragments` is `True`. (default: 1)
                """

                from ....utilities.reactive_complexes.inter_reactive_complexes import InterReactiveComplexes

                self.complex_generator = InterReactiveComplexes()
                """
                InterReactiveComplexes
                    The generator used for the composition of reactive complexes
                    from given interstructural components of reactive coordinates.
                """
                self.minimal_spin_multiplicity = False
                """
                bool
                    Whether to assume max spin recombination, thus assuming minimal resulting spin, or take combination
                    of input spin multiplicities. (default: False)
                    True: | multiplicity1 - multiplicity2 | - 1
                    False: sum(multiplicities) - 1
                """

        class UnimolDissociationOptions:
            """
            The options for the generation of dissociative and
            disconnective unimolecular trial elementary step calculations.
            """

            __slots__ = (
                "enabled",
                "job",
                "max_graph_distance",
                "job_settings_dissociative",
                "job_settings_disconnective",
            )

            def __init__(self):
                self.enabled: bool = True
                """
                bool
                    If `True`, enables the exploration of these types of reactions.
                """
                self.job: db.Job = db.Job("scine_react_complex_nt")
                """
                db.Job (Scine::Database::Calculation::Job)
                    The Job used to evaluate the possible reactions.
                    The default is: the `scine_react_complex_nt` order on a single core.
                """
                self.job_settings_dissociative: utils.ValueCollection = utils.ValueCollection({})
                """
                utils.ValueCollection
                    Additional settings passed to Calculation evaluating the possible
                    reactions. These settings are passed to all calculations that are
                    expected to result in dissociations of the tested atoms only,
                    retaining a single molecule.
                    Empty by default.
                """
                self.job_settings_disconnective: utils.ValueCollection = utils.ValueCollection({})
                """
                utils.ValueCollection
                    Additional settings passed to Calculation evaluating the possible
                    reactions. These settings are passed to all calculations that are
                    expected to result in disconnection of the molecule generating
                    two new molecules in the process.
                    Empty by default.
                """

        class UnimolAssociationOptions:
            """
            The options for the generation of associative unimolecular elementary step trial calculations.
            """

            __slots__ = (
                "enabled",
                "job",
                "job_settings",
                "min_inter_fragment_graph_distance",
                "max_inter_fragment_graph_distance",
                "consider_diatomic_fragments",
                "max_within_fragment_graph_distance",
            )

            def __init__(self):
                self.enabled = True
                """
                bool
                    If `True`, enables the exploration of these types of reactions.
                """
                self.job: db.Job = db.Job("scine_react_complex_nt")
                """
                db.Job (Scine::Database::Calculation::Job)
                    The Job used to evaluate the possible reactions.
                    The default is: the `scine_react_complex_afir` order on a single core.
                """
                self.job_settings: utils.ValueCollection = utils.ValueCollection({})
                """
                utils.ValueCollection
                    Additional settings passed to Calculation evaluating the possible
                    reactions.
                    Empty by default.
                """
                self.min_inter_fragment_graph_distance = 2
                """
                int
                    The minimum graph distance between fragments (bond count) to
                    be considered for association reactions.
                    If the fragments are diatomic all possible atom pairs
                    arising from forming the cartesian product of the fragments
                    have to have a distance of at least this value.
                """
                self.max_inter_fragment_graph_distance = 6
                """
                int
                    The maximum distance between atoms (bond count) to be considered
                    for association reactions.
                    If the fragments are diatomic all possible atom pairs
                    arising from forming the cartesian product of the fragments
                    have to have a distance that is not larger than this value.
                """
                self.consider_diatomic_fragments = False
                """
                bool
                    Whether to also set up unimolecular reactions where the individual reactive fragments are composed
                    of two atoms. (default: False)
                    Note: Pairs are only made subject to an elementary step trial/reactive complex calculations if none
                    of the atoms of the one fragment is directly bound to any of the atoms in the other fragment.
                """
                self.max_within_fragment_graph_distance = 1
                """
                int
                    Maximum number of bonds between two atoms forming a reactive
                    fragment. This option is only effective if
                    `consider_diatomic_fragments` is `True`. (default: 1)
                """

        def __init__(self, parent: Optional[TrialGenerator] = None):
            super().__init__(parent)
            self.unimolecular_dissociation_options = self.UnimolDissociationOptions()
            """
            UnimolDissociationOptions
                The settings for reactions involving a single molecule, that are
                set up to be dissociative in nature.
            """
            self.unimolecular_association_options = self.UnimolAssociationOptions()
            """
            UnimolAssociationOptions
                The settings for reactions involving a single molecule, that are
                set up to be associative in nature.
            """
            self.bimolecular_association_options = self.BimolAssociationOptions()
            """
            BimolAssociationOptions
                The settings for reactions involving two molecules, that are set up
                to be associative in nature.
            """

    def clear_cache(self):
        pass

    @_sanity_check_wrapper
    def bimolecular_coordinates(self, structure_list: List[db.Structure], with_exact_settings_check: bool = False) \
            -> Dict[
        Tuple[List[Tuple[int, int]], int],
        List[Tuple[ndarray, ndarray, float, float]]
    ]:
        # Check number of compounds
        if len(structure_list) != 2:
            raise RuntimeError("Exactly two structures are needed for setting up an intermolecular reaction.")
        if self._quick_bimolecular_already_exists(structure_list, do_fast_query=not with_exact_settings_check):
            return {}
        result = defaultdict(list)
        for (
                inter_pairs,
                align1,
                align2,
                rot,
                spread,
        ) in self._generate_bimolecular_coordinates(structure_list):
            # 0 because n_diss is always 0 for fragment based
            result[(inter_pairs, 0)].append((align1, align2, rot, spread))
        return result

    def _quick_bimolecular_already_exists(self, structure_list: List[db.Structure], do_fast_query: bool) -> bool:
        """
            If there is a reactive complex calculation for the same structures, return True
        """
        if not do_fast_query:
            return False
        structure_id_list = [structure.id() for structure in structure_list]
        return calculation_exists_in_structure(
            self.get_bimolecular_job_order(),
            structure_id_list,
            self.options.model,
            self._structures,
            self._calculations
        )

    @_sanity_check_wrapper
    def bimolecular_reactions(self, structure_list: List[db.Structure], with_exact_settings_check: bool = False):
        """
        Creates reactive complex calculations corresponding to the bimolecular
        reactions between the structures if there is not already a calculation
        to search for a reaction of the same structures with the same job order.

        Parameters
        ----------
        structure_list :: List[db.Structure]
            List of the two structures to be considered.
            The Structures have to be linked to a database.
        with_exact_settings_check :: bool
            If True, more expensive queries are carried out to check if the settings of the
            calculations are exactly the same as the settings of the trial generator. This allows to add more
            inclusive additional reaction trials but the queries are less efficient, therefore this option
            should be only toggled if necessary.
        """
        # Check number of compounds
        if len(structure_list) != 2:
            raise RuntimeError("Exactly two structures are needed for setting up an intermolecular reaction.")

        structure_id_list = [s.id() for s in structure_list]

        # If there is a reactive complex calculation for the same structures, return
        if self._quick_bimolecular_already_exists(structure_list, do_fast_query=not with_exact_settings_check):
            return

        new_calculation_ids = []
        # add calculations
        for (
                inter_pairs,
                align1,
                align2,
                rot,
                spread,
        ) in self._generate_bimolecular_coordinates(structure_list):
            # First element of interstructural pair belongs to first structure and second to second
            # Indexes in the rhs and lhs list refer to the individual structures so no shifting necessary
            lhs = list(set(pair[0] for pair in inter_pairs))
            rhs = list(set(pair[1] for pair in inter_pairs))
            # Set up calculation because no further intrastructural components supported
            cid = self._add_reactive_complex_calculation(
                structure_id_list,
                ReactionType.Associative,
                lhs,
                rhs,
                self.options.bimolecular_association_options.job,
                self.options.bimolecular_association_options.job_settings,
                list(align1),
                list(align2),
                rot,
                spread,
                0.0,
                check_for_existing=with_exact_settings_check
            )
            if cid is not None:
                new_calculation_ids.append(cid)
        if new_calculation_ids:
            for s in structure_list:
                s.add_calculations(self.get_bimolecular_job_order(), [new_calculation_ids[0]])

    def _generate_bimolecular_coordinates(self, structure_list: List[db.Structure]) -> Generator[
        Tuple[List[Tuple[int, int]], ndarray, ndarray, float, float], None, None
    ]:
        # Get single reactive atoms per structure
        n_atoms1 = structure_list[0].get_atoms().size()
        n_atoms2 = structure_list[1].get_atoms().size()
        reactive_atoms = self.reactive_site_filter.filter_atoms(
            structure_list, list(range(n_atoms1)) + [i + n_atoms1 for i in range(n_atoms2)]
        )
        reactive_monoatomic1 = [i for i in reactive_atoms if i < n_atoms1]
        reactive_monoatomic2 = [i - n_atoms1 for i in reactive_atoms if i >= n_atoms1]
        # Get reactive pairs within each structure in agreement with distance settings
        reactive_diatomic1 = []
        reactive_diatomic2 = []

        if self.options.bimolecular_association_options.consider_diatomic_fragments:
            connectivity_analyzer1 = ConnectivityAnalyzer(structure_list[0])
            connectivity_analyzer2 = ConnectivityAnalyzer(structure_list[1])
            reactive_diatomic1 = self._get_intrastructural_pairs(
                connectivity_analyzer1,
                reactive_monoatomic1,
                (1, self.options.bimolecular_association_options.max_within_fragment_graph_distance),
            )
            reactive_diatomic2 = self._get_intrastructural_pairs(
                connectivity_analyzer2,
                reactive_monoatomic2,
                (1, self.options.bimolecular_association_options.max_within_fragment_graph_distance),
            )

        # Generate reactive pair combinations across the two structures from reactive sites
        # Applies the reactive_site_filter's filter_atom_pair method
        inter_coords = list(
            self._generate_inter_reactive_coords(
                structure_list, reactive_monoatomic1, reactive_monoatomic2, reactive_diatomic1, reactive_diatomic2
            )
        )
        n_atoms1 = structure_list[0].get_atoms().size()
        # TODO Unify indexing logic to avoid shifting
        # Shift to complex indices
        shifted_inter_coords = self._shift_reaction_coordinates(inter_coords, n_atoms1)
        # Filter
        filtered_inter_coords = self.reactive_site_filter.filter_reaction_coordinates(
            structure_list, shifted_inter_coords
        )
        # Shift back
        inter_coords = self._shift_reaction_coordinates(filtered_inter_coords, -n_atoms1)

        # For each coord generate the corresponding reactive complexes
        for (
            inter_pairs,
            align1,
            align2,
            rot,
            spread,
        ) in self.options.bimolecular_association_options.complex_generator.generate_reactive_complexes(
            structure_list[0], structure_list[1], inter_coords
        ):
            yield inter_pairs, align1, align2, rot, spread

    def _quick_unimolecular_already_exists(self, structure: db.Structure, do_fast_query: bool) -> bool:
        """
            If there is a reactive complex calculation for the same structures, return True
        """
        # Rule out compounds too small for intramolecular reactions right away
        atoms = structure.get_atoms()
        # No intramolecular reactions for monoatomic compounds
        if atoms.size() == 1:
            return True
        # Only consider diatomic structures if dissociations are to be considered
        if atoms.size() == 2 and not self.options.unimolecular_dissociation_options.enabled:
            return True
        # check later for complete check
        if not do_fast_query:
            return False
        return calculation_exists_in_structure(
            self.get_unimolecular_job_order(),
            [structure.id()],
            self.options.model,
            self._structures,
            self._calculations
        )

    @_sanity_check_wrapper
    def unimolecular_coordinates(self, structure: db.Structure, with_exact_settings_check: bool = False) \
            -> List[Tuple[List[List[Tuple[int, int]]], int]]:
        if self._quick_unimolecular_already_exists(structure, do_fast_query=not with_exact_settings_check):
            return []

        result = []

        atoms = structure.get_atoms()
        connectivity_analyzer = ConnectivityAnalyzer(structure)
        reactive_atoms = self.reactive_site_filter.filter_atoms([structure], list(range(atoms.size())))

        # Associative reactions
        if self.options.unimolecular_association_options.enabled:
            # Atom on atom
            # Pairs allowed wrt distance bound options
            monoatomic_fragment_pairs = self._get_intrastructural_pairs(
                connectivity_analyzer,
                reactive_atoms,
                (
                    self.options.unimolecular_association_options.min_inter_fragment_graph_distance,
                    self.options.unimolecular_association_options.max_inter_fragment_graph_distance,
                ),
            )
            filtered_monoatomic_fragment_pairs = self.reactive_site_filter.filter_atom_pairs(
                [structure], monoatomic_fragment_pairs
            )
            atom_on_atom_coordinates = [[pair] for pair in filtered_monoatomic_fragment_pairs]
            filtered_coords = self.reactive_site_filter.filter_reaction_coordinates([structure],
                                                                                    atom_on_atom_coordinates)
            for coord in filtered_coords:
                # sanity check
                if len(coord) != 1:
                    # Should not be reached
                    raise RuntimeError("Monoatomic fragment coordinate has to be of length one.")
            result.append((filtered_coords, 0))
            # Diatomic fragments
            if self.options.unimolecular_association_options.consider_diatomic_fragments:
                diatomic_fragments = self._get_intrastructural_pairs(
                    connectivity_analyzer,
                    reactive_atoms,
                    (1, self.options.unimolecular_association_options.max_within_fragment_graph_distance),
                )
                # Monoatomic to diatomic (e.g. "atom on bond" if max_within_fragment_graph_distance==1)
                for atom in reactive_atoms:
                    for frag in diatomic_fragments:
                        if atom in frag:
                            continue
                        # Check whether both atoms of the diatomic fragment respect the distance bounds wrt the
                        # monoatomic fragment
                        if not self._check_all_inter_fragment_distances(
                                connectivity_analyzer,
                                [atom],
                                list(frag),
                                (
                                    self.options.unimolecular_association_options.min_inter_fragment_graph_distance,
                                    self.options.unimolecular_association_options.max_inter_fragment_graph_distance,
                                ),
                        ):
                            continue
                        coord = [(atom, frag[0]), (atom, frag[1])]
                        # Both pairs have to be reactive
                        if len(self.reactive_site_filter.filter_atom_pairs([structure], coord)) != 2:
                            continue
                        # Coordinate has to be reactive
                        if len(self.reactive_site_filter.filter_reaction_coordinates([structure], [coord])) != 1:
                            continue
                        if connectivity_analyzer.get_reaction_type(coord) == ReactionType.Associative:
                            result.append(([coord], 0))
                # Diatomic to diatomic (e.g. "bond on bond" if max_within_fragment_graph_distance==1)
                for frag1, frag2 in combinations(diatomic_fragments, 2):
                    if len(set(frag1 + frag2)) < 4:
                        # Atoms shall not appear on lhs AND rhs at the same time
                        continue
                    # All possible pair interpretations between the two fragments shall respect the distance bounds
                    if not self._check_all_inter_fragment_distances(
                            connectivity_analyzer,
                            frag1,
                            frag2,
                            (
                                self.options.unimolecular_association_options.min_inter_fragment_graph_distance,
                                self.options.unimolecular_association_options.max_inter_fragment_graph_distance,
                            ),
                    ):
                        continue

                    coord = [(frag1[0], frag2[0]), (frag1[1], frag2[1])]
                    shuffled_coord = [(frag1[0], frag2[1]), (frag1[1], frag2[0])]
                    # Check whether any of the reactive pair interpretations is reactive
                    if (
                            len(self.reactive_site_filter.filter_atom_pairs([structure], coord)) < 2
                            and len(self.reactive_site_filter.filter_atom_pairs([structure], shuffled_coord)) < 2
                    ):
                        continue
                    # Check whether any of the reaction coordinate interpretations is reactive
                    if (
                            len(self.reactive_site_filter.filter_reaction_coordinates([structure],
                                                                                      [coord, shuffled_coord]))
                            == 0
                    ):
                        continue
                    # Only consider if there is no direct bond between fragments
                    if connectivity_analyzer.get_reaction_type(coord) != ReactionType.Associative:
                        continue
                    if connectivity_analyzer.get_reaction_type(shuffled_coord) != ReactionType.Associative:
                        continue
                    result.append(([coord], 0))

        # Dissociative reactions
        if self.options.unimolecular_dissociation_options.enabled:
            # Only adjacent (i.e. with a graph distance of 1) atoms are made subject to dissociation
            diss_pairs = self._get_intrastructural_pairs(connectivity_analyzer, reactive_atoms, (1, 1))
            # Filter wrt pair filter
            filtered_diss_pairs = self.reactive_site_filter.filter_atom_pairs([structure], diss_pairs)
            # Filter wrt coordinate filter
            filtered_coords = self.reactive_site_filter.filter_reaction_coordinates(
                [structure], [[pair] for pair in filtered_diss_pairs]
            )
            for coord in filtered_coords:
                # Dissociative coordinates only consist of one pair each
                if len(coord) > 1:
                    # Should not be reached
                    raise RuntimeError(
                        "Only one simultaneous dissociation supported in fragment based trial generator."
                    )
            result.append((filtered_coords, 1))

        return result

    def get_unimolecular_job_order(self) -> str:
        if self.options.unimolecular_dissociation_options.enabled and \
                self.options.unimolecular_association_options.enabled:
            assert self.options.unimolecular_dissociation_options.job.order == \
                self.options.unimolecular_association_options.job.order
            return self.options.unimolecular_association_options.job.order
        elif self.options.unimolecular_dissociation_options.enabled:
            return self.options.unimolecular_dissociation_options.job.order
        elif self.options.unimolecular_association_options.enabled:
            return self.options.unimolecular_association_options.job.order
        else:
            raise RuntimeError('No reactions enabled, cannot report job order.')

    def get_bimolecular_job_order(self) -> str:
        return self.options.bimolecular_association_options.job.order

    @_sanity_check_wrapper
    def unimolecular_reactions(self, structure: db.Structure, with_exact_settings_check: bool = False):
        """
        Creates reactive complex calculations corresponding to the unimolecular
        reactions of the structure if there is not already a calculation to
        search for a reaction of the same structure with the same job order.

        Parameters
        ----------
        structure :: db.Structure
            The structure to be considered. The Structure has to
            be linked to a database.
        with_exact_settings_check :: bool
            If True, more expensive queries are carried out to check if the settings of the
            calculations are exactly the same as the settings of the trial generator. This allows to add more
            inclusive additional reaction trials but the queries are less efficient, therefore this option
            should be only toggled if necessary.
        """
        if self.options.unimolecular_dissociation_options.enabled and \
                self.options.unimolecular_association_options.enabled:
            assert self.options.unimolecular_dissociation_options.job.order == \
                self.options.unimolecular_association_options.job.order
        if self._quick_unimolecular_already_exists(structure, do_fast_query=not with_exact_settings_check):
            return
        structure_id = structure.id()

        # Check whether there is a reactive complex calculation for the same structure already
        if self.options.unimolecular_dissociation_options.enabled:
            selection = {
                "$and": [
                    {"job.order": {"$eq": self.options.unimolecular_dissociation_options.job.order}},
                    {"structures": [{"$oid": structure_id.string()}]},
                    *model_query(self.options.model)
                ]
            }
            if self._calculations.get_one_calculation(dumps(selection)) is not None:
                return
        if self.options.unimolecular_association_options.enabled:
            selection = {
                "$and": [
                    {"job.order": {"$eq": self.options.unimolecular_association_options.job.order}},
                    {"structures": [{"$oid": structure_id.string()}]},
                    *model_query(self.options.model)
                ]
            }
            if self._calculations.get_one_calculation(dumps(selection)) is not None:
                return

        reaction_coordinates = self.unimolecular_coordinates(structure, with_exact_settings_check)
        connectivity_analyzer = ConnectivityAnalyzer(structure)

        def n_diss_error(r_type: ReactionType, n: int) -> str:
            return f"Analyzed reaction type to be {r_type}, but reaction coordinate " \
                   f"specifies {n} dissociations, this should not be possible in " \
                   f"{self.__class__.__name__}"

        new_calculation_ids = []
        for coordinates, n_diss in reaction_coordinates:
            for coord in coordinates:
                reaction_type = connectivity_analyzer.get_reaction_type(coord)
                if reaction_type == ReactionType.Associative:
                    if n_diss != 0:
                        raise RuntimeError(n_diss_error(reaction_type, n_diss))
                    job = self.options.unimolecular_association_options.job
                    settings = self.options.unimolecular_association_options.job_settings
                elif reaction_type == ReactionType.Dissociative:
                    if n_diss == 0:
                        raise RuntimeError(n_diss_error(reaction_type, n_diss))
                    job = self.options.unimolecular_dissociation_options.job
                    settings = self.options.unimolecular_dissociation_options.job_settings_dissociative
                elif reaction_type == ReactionType.Disconnective:
                    if n_diss == 0:
                        raise RuntimeError(n_diss_error(reaction_type, n_diss))
                    job = self.options.unimolecular_dissociation_options.job
                    settings = self.options.unimolecular_dissociation_options.job_settings_disconnective
                else:
                    raise RuntimeError(f"Unknown reaction type {reaction_type}")
                lhs = list(set([c[0] for c in coord]))
                rhs = list(set([c[1] for c in coord]))
                cid = self._add_reactive_complex_calculation(
                    [structure_id],
                    reaction_type,
                    lhs,
                    rhs,
                    job,
                    settings,
                    check_for_existing=with_exact_settings_check
                )
                if cid is not None:
                    new_calculation_ids.append(cid)
        if new_calculation_ids:
            structure.add_calculations(self.get_unimolecular_job_order(), [new_calculation_ids[0]])

    def _add_reactive_complex_calculation(
        self,
        reactive_structures: List[db.ID],
        reaction_type: ReactionType,
        lhs_list: List[int],
        rhs_list: List[int],
        job: db.Job,
        settings: utils.ValueCollection,
        lhs_alignment: Optional[List[float]] = None,
        rhs_alignment: Optional[List[float]] = None,
        x_rotation: Optional[float] = None,
        spread: Optional[float] = None,
        displacement: Optional[float] = None,
        check_for_existing: bool = False,
    ) -> Optional[db.ID]:
        """
        Adds a reactive calculation to the database and puts it on hold.

        Parameters
        ----------

        reactive_structures :: List[db.ID]
            List of the IDs of the reactants and, for intermolecular reactions,
            the reactive complex.
        reaction_type :: ReactionType
            Type of the reaction to be set up.
        lhs_list, rhs_list :: List[int]
            Indices of the reactive sites within the reactive complex.
        job :: since_database.Job
        settings :: scine_utilities.ValueCollection
        lhs_alignment :: List[float], length=9
            In case of two structures building the reactive complex, this option
            describes a rotation of the first structure (index 0) that aligns
            the reaction coordinate along the x-axis (pointing towards +x).
            The rotation assumes that the geometric mean position of all
            atoms in the reactive site (``lhs_list``) is shifted into the
            origin.
        rhs_alignment :: List[float], length=9
            In case of two structures building the reactive complex, this option
            describes a rotation of the second structure (index 1) that aligns
            the reaction coordinate along the x-axis (pointing towards -x).
            The rotation assumes that the geometric mean position of all
            atoms in the reactive site (``rhs_list``) is shifted into the
            origin.
        x_rotation :: float
            In case of two structures building the reactive complex, this option
            describes a rotation angle around the x-axis of one of the two
            structures after ``lhs_alignment`` and ``rhs_alignment`` have
            been applied.
        spread :: float
            In case of two structures building the reactive complex, this option
            gives the distance by which the two structures are moved apart along
            the x-axis after ``lhs_alignment``, ``rhs_alignment``, and
            ``x_rotation`` have been applied.
        displacement :: float
            In case of two structures building the reactive complex, this option
            adds a random displacement to all atoms (random direction, random
            length). The maximum length of this displacement (per atom) is set to
            be the value of this option.
        check_for_existing :: bool
            Whether it should be checked if a calculation with these exact
            settings and model already exists or not (default: False)
        """
        model = self.options.model
        this_settings = self._get_settings(settings)
        # Sleep a bit in order not to make the DB choke
        time.sleep(0.001)
        if lhs_alignment is not None:
            this_settings["rc_x_alignment_0"] = lhs_alignment
        if rhs_alignment is not None:
            this_settings["rc_x_alignment_1"] = rhs_alignment
        if x_rotation is not None:
            this_settings["rc_x_rotation"] = x_rotation
        if spread is not None:
            this_settings["rc_x_spread"] = spread
        if displacement is not None:
            this_settings["rc_displacement"] = displacement
        if job.order == "scine_react_complex_afir":
            this_settings["afir_afir_rhs_list"] = rhs_list
            this_settings["afir_afir_lhs_list"] = lhs_list
            # Repulsive for dissociations
            this_settings["afir_afir_attractive"] = not bool(reaction_type.value)
            # Use maximum fragment distance if disconnective
            this_settings["afir_afir_use_max_fragment_distance"] = bool(reaction_type == ReactionType.Disconnective)
            if reaction_type == ReactionType.Disconnective:
                # Set fragment distance for convergence to 3*sum of maximum covalent radii within i and j
                struct = db.Structure(reactive_structures[0])
                struct.link(self._structures)
                atoms = struct.get_atoms()
                covalent_max = [
                    max(utils.ElementInfo.covalent_radius(atoms.get_element(atom_index)) for atom_index in frag)
                    for frag in (rhs_list, lhs_list)
                ]
                this_settings["afir_afir_max_fragment_distance"] = 3 * sum(covalent_max)
            else:
                del this_settings["afir_afir_max_fragment_distance"]

        elif job.order == "scine_react_complex_nt":
            this_settings["nt_nt_rhs_list"] = rhs_list
            this_settings["nt_nt_lhs_list"] = lhs_list
            # Repulsive for dissociations; attractive (False) for dissociative/disconnective
            this_settings["nt_nt_attractive"] = not bool(reaction_type.value)
        else:
            raise RuntimeError(
                "Only 'scine_react_complex_afir' and 'scine_react_complex_nt' "
                "order supported for fragment-based reactive complex calculations."
            )

        if len(reactive_structures) > 1:
            this_settings["rc_minimal_spin_multiplicity"] = bool(
                self.options.bimolecular_association_options.minimal_spin_multiplicity
            )

        if check_for_existing and get_calculation_id(job.order, reactive_structures, model,
                                                     self._calculations, settings=this_settings) is not None:
            return None

        calculation = db.Calculation()
        calculation.link(self._calculations)
        calculation.create(model, job, reactive_structures)
        calculation.set_settings(this_settings)
        calculation.set_status(db.Status.HOLD)
        return calculation.id()

    def _generate_inter_reactive_coords(
        self,
        structure_list: List[db.Structure],
        reactive_monoatomic1: List[int],
        reactive_monoatomic2: List[int],
        reactive_diatomic1: List[Tuple[int, int]],
        reactive_diatomic2: List[Tuple[int, int]],
    ) -> Iterator[Tuple[Tuple[int, int], ...]]:
        """
        Enumerates all reactive coordinates with all pairs involving both
        structures arising from the given monoatomic and diatomic reactive
        fragments of two structures.
        Yields all coordinates whose reactive pairs pass the
        `filter_atom_pair` method of the reactive site filter.
        For two diatomic fragments the first atoms of both fragments and the
        second atoms of both fragments are interpreted as one pair each.
        If one fragment is diatomic and the other monoatomic the two pairs
        arising from combining the monoatomic fragment atom with each of the
        atoms from the diatomic fragment have to pass the pair filter.

        Parameters
        ----------
        reactive_monoatomic1 : List[int]
            Indices of reactive atoms within the first structure.
        reactive_monoatomic2 : List[int]
            Indices of reactive atoms within the second structure.
        reactive_diatomic1 : List[Tuple[int]]
            Indices of diatomic reactive fragments (e.g. bonds) within the first structure.
        reactive_diatomic2 : List[Tuple[int]]
            Indices of diatomic reactive fragments (e.g. bonds) within the second structure.

        Yields
        -------
        Tuple[Tuple[int, int], ...]
            The interstructural reactive coordinate expressed as lists of atom pairs with
            the first element of each pair referring to an atom within the first structure
            and the second to the second structure.
        """
        # Atom on atom
        n_atoms1 = structure_list[0].get_atoms().size()
        for atom1 in reactive_monoatomic1:
            for atom2 in reactive_monoatomic2:
                # Yield for easy avoidance of whole list storage when one complex at a
                # time implemented in inter_reactive_complex generation
                if len(self.reactive_site_filter.filter_atom_pairs(structure_list, [(atom1, atom2 + n_atoms1)])) > 0:
                    yield ((atom1, atom2),)

        # Monoatomic on diatomic fragment (e.g. "atom on bond" if max_within_fragment_graph_distance==1)
        # Diatomic fragment from structure 1 and atom from structure 2
        for frag1 in reactive_diatomic1:
            for atom2 in reactive_monoatomic2:
                # Both atoms of frag1 shall be reactive wrt atom2
                if (
                    len(
                        self.reactive_site_filter.filter_atom_pairs(
                            structure_list, [(frag1[0], atom2 + n_atoms1), (frag1[1], atom2 + n_atoms1)]
                        )
                    )
                    == 2
                ):
                    yield (
                        (frag1[0], atom2),
                        (frag1[1], atom2),
                    )
        # Diatomic fragment from structure 2 and atom from structure 1
        for frag2 in reactive_diatomic2:
            for atom1 in reactive_monoatomic1:
                # Both atoms of frag2 shall be reactive wrt atom1
                if (
                    len(
                        self.reactive_site_filter.filter_atom_pairs(
                            structure_list, [(atom1, frag2[0] + n_atoms1), (atom1, frag2[1] + n_atoms1)]
                        )
                    )
                    == 2
                ):
                    yield (
                        (atom1, frag2[0]),
                        (atom1, frag2[1]),
                    )
        # Diatomic on diatomic (e.g. "bond on bond" if max_within_fragment_graph_distance==1)
        # Diatomic on diatomic fragment (e.g. "bond on bond" if max_within_fragment_graph_distance==1)
        for frag1 in reactive_diatomic1:
            for frag2 in reactive_diatomic2:
                # Note: Pair matching makes difference because reactive pairs are aligned in structure generation
                # Both pairs shall be reactive
                if (
                    len(
                        self.reactive_site_filter.filter_atom_pairs(
                            structure_list, [(frag1[0], frag2[0] + n_atoms1), (frag1[1], frag2[1] + n_atoms1)]
                        )
                    )
                    == 2
                ):
                    yield (
                        (frag1[0], frag2[0]),
                        (frag1[1], frag2[1]),
                    )
                if (
                    len(
                        self.reactive_site_filter.filter_atom_pairs(
                            structure_list, [(frag1[0], frag2[1] + n_atoms1), (frag1[1], frag2[0] + n_atoms1)]
                        )
                    )
                    == 2
                ):
                    yield (
                        (frag1[0], frag2[1]),
                        (frag1[1], frag2[0]),
                    )

    @staticmethod
    def _get_intrastructural_pairs(
        connectivity_analyzer: ConnectivityAnalyzer, reactive_atoms: List[int], distance_bounds: Tuple[int, int]
    ):
        """
        Generate a list of pairs of the specified reactive atoms whose graph distance is within the given bounds.

        Parameters
        ----------
        connectivity_analyzer : ConnectivityAnalyzer
            A connectivity analyzer holding the molecule etc. corresponding to the structure.
        reactive_atoms : List[int]
            The atoms within the structure that should be combined.
        distance_bounds : Tuple[int]
            The minimum and maximum graph distance within pairs that is allowed.

        Returns
        -------
        List[Tuple[int, int]]
            A list of pairs of the provided atoms
        """
        pairs = []
        for atom1, atom2 in combinations(reactive_atoms, 2):  # Loop over unique pairs
            if distance_bounds[0] <= connectivity_analyzer.get_graph_distance(atom1, atom2) <= distance_bounds[1]:
                pairs.append((atom1, atom2))
        return pairs

    @staticmethod
    def _check_all_inter_fragment_distances(
        connectivity_analyzer: ConnectivityAnalyzer, frag1: List[int], frag2: List[int],
        distance_bounds: Tuple[int, int]
    ):
        """
        Checks whether all atom pairs arising between two reactive
        fragments (i.e. all pairs out of the product of the two fragment lists)
        have graph distances within the specified boundaries.

        Parameters
        ----------
        connectivity_analyzer : ConnectivityAnalyzer
            The connectivity analyzer of the structure of interest
        frag1 : List[int]
            The first reactive fragment.
        frag2 : List[int]
            The second reactive fragment.
        distance_bounds : Tuple[int]
            The minimum and maximum graph distance in between the fragments that
            is allowed.

        Returns
        -------
        bool
            `True` if all inter fragment pairs have a graph distance within the
            specified boundaries, otherwise `False`.
        """
        for inter_fragment_pair in product(frag1, frag2):
            dist = connectivity_analyzer.get_graph_distance(inter_fragment_pair[0], inter_fragment_pair[1])
            if not (distance_bounds[0] <= dist <= distance_bounds[1]):
                return False
        return True

    @staticmethod
    def _shift_reaction_coordinates(reaction_coordinates: List[Tuple[Tuple[int, int], ...]], shift: int):
        """
        Shifts the second member of all reactive pairs of all reaction
        coordinates by the specified value.

        Parameters
        ----------
        reaction_coordinates : List[Tuple[Tuple[int, int], ...]]
            The reaction coordinates of interest.
        shift : int
            The shift to be applied

        Returns
        -------
        List[Tuple[Tuple[int, int], ...]]
            The shifted reaction coordinate.
        """
        shifted_coords = []
        for coord in reaction_coordinates:
            shifted_coords.append(tuple((pair[0], pair[1] + shift) for pair in coord))

        return shifted_coords
