deprotonation

@dataclass
class CRESTDeprotonation[InputType: Molecule, OutputType: Molecule](
    DeprotonationMaker, CRESTCalculator, PymatGenMaker[InputType, OutputType]
):
    """Generate deprotonated structures using CREST.

    Uses CREST's automated deprotonation workflow to identify acidic sites
    and generate low-energy deprotonated structures. The method systematically
    explores different deprotonation sites and optimizes the resulting structures
    using GFN2-xTB.

    Attributes:
        name: Name of the job (default: "CREST Deprotonation").
        ewin: Energy window [kcal/mol] for selecting deprotonated structures
            (default: None, uses CREST default). Structures within ewin of the
            lowest energy structure are retained.

    References:
        - CREST Documentation: https://crest-lab.github.io/crest-docs/

    Examples:
        >>> from jfchemistry.modification import CRESTDeprotonation # doctest: +SKIP
        >>> from pymatgen.core import Molecule # doctest: +SKIP
        >>> from ase.build import molecule # doctest: +SKIP
        >>> ethane = Molecule.from_ase_atoms(molecule("C2H6")) # doctest: +SKIP
        >>> # Deprotonate a ethane
        >>> deprot = CRESTDeprotonation(ewin=6.0) # doctest: +SKIP
        >>> job = deprot.make(ethane) # doctest: +SKIP
        >>> deprotonated_structures = job.output["structure"] # doctest: +SKIP
        >>>
        >>> deprot_default = CRESTDeprotonation() # doctest: +SKIP
        >>> job = deprot_default.make(ethane) # doctest: +SKIP
    """

    name: str = "CREST Deprotonation"

    # INTERNAL
    _runtype: Literal["deprotonate"] = "deprotonate"
    _output_filename: str = "deprotonated.xyz"

    def _make_commands(self):
        """Make the CLI for the CREST input."""
        super()._make_commands()
        self._commands.append(f"--{self._runtype}")
        self._commands.append("--newversion")

    def _operation(
        self, input: InputType, **kwargs
    ) -> tuple[OutputType | list[OutputType], Properties | list[Properties] | None]:
        """Generate deprotonated structures using CREST.

        Runs CREST's deprotonation workflow to identify acidic sites and
        generate optimized deprotonated structures. The calculation uses
        GFN2-xTB with Wiberg bond order analysis.

        Args:
            input: Input molecular structure with 3D coordinates. The
                molecule's charge is used for the CREST calculation.
            **kwargs: Additional kwargs to pass to the operation.

        Returns:
            Tuple containing:
                - List of deprotonated structures sorted by energy
                - None (no additional properties)

        Examples:
            >>> from pymatgen.core import Molecule # doctest: +SKIP
            >>> from ase.build import molecule # doctest: +SKIP
            >>> ethane = Molecule.from_ase_atoms(molecule("C2H6")) # doctest: +SKIP
            >>> deprot = CRESTDeprotonation(ewin=8.0) # doctest: +SKIP
            >>> structures, props = deprot.operation(ethane) # doctest: +SKIP
            >>> print(f"Generated {len(structures)} deprotonated structures") # doctest: +SKIP
        """
        input.to("input.xyz", fmt="xyz")
        if self.charge is None and input.charge is not None:
            self.charge = input.charge
        super()._make_dict()
        super()._write_toml()
        self._make_commands()
        super()._run()
        if not os.path.exists(self._output_filename):
            raise FileNotFoundError(
                "No deprotonated structures found. Please check your CREST settings and log file."
            ) from None
        molecules = XYZ.from_file(self._output_filename).all_molecules
        for i, deprotonated_structure in enumerate(molecules):
            molecules[i] = deprotonated_structure.set_charge_and_spin(
                charge=deprotonated_structure.charge - 1,
                spin_multiplicity=int(deprotonated_structure.charge - 1) // 2 + 1,
            )
        return cast("list[OutputType]", molecules), cast("list[Properties]", Properties())