Merge pull request #225 from duartegroup/v1.3.4

V1.3.4

.github/ISSUE_TEMPLATE/feature_request.md

@@ -0,0 +1,11 @@
+---
+name: Feature request
+about: Request a feature to be added
+title: ''
+labels: ''
+assignees: ''
+
+---
+
+**Describe the feature**
+<!-- A clear and concise description of what the feature is. -->

.github/pull_request_template.md

@@ -0,0 +1,13 @@
+<!--
+Please replace _#ISSUE_ with just e.g. #12 if this PR resolves issue 12.
+-->
+## Resolves _#ISSUE_
+
+
+***
+
+## Checklist
+
+* [ ] The changes include an associated explanation of how/why
+* [ ] Test pass
+* [ ] Documentation has been updated

.pre-commit-config.yaml

@@ -4,4 +4,3 @@ repos:
     hooks:
       - id: black
         language_version: python3.9
-

autode/__init__.py

@@ -7,6 +7,7 @@
 from autode import hessians
 from autode.reactions.reaction import Reaction
 from autode.reactions.multistep import MultiStepReaction
+from autode.transition_states.transition_state import TransitionState
 from autode.atoms import Atom
 from autode.species.molecule import Reactant, Product, Molecule, Species
 from autode.species.complex import NCIComplex
@@ -38,7 +39,7 @@
   - Merge when tests pass
 """
 
-__version__ = "1.3.3"
+__version__ = "1.3.4"
 
 
 __all__ = [
@@ -55,6 +56,7 @@
     "Reactant",
     "Product",
     "Molecule",
+    "TransitionState",
     "NCIComplex",
     "Config",
     "Calculation",

autode/atoms.py

@@ -1,6 +1,6 @@
 import numpy as np
 from copy import deepcopy
-from typing import Union, Optional, List, Sequence
+from typing import Union, Optional, List, Sequence, Any
 from autode.log import logger
 from autode.geom import get_rot_mat_euler
 from autode.values import (
@@ -79,6 +79,21 @@ def __repr__(self):
     def __str__(self):
         return self.__repr__()
 
+    def __eq__(self, other: Any):
+        """Equality of another atom to this one"""
+        are_equal = (
+            isinstance(other, Atom)
+            and other.label == self.label
+            and other.atom_class == self.atom_class
+            and isinstance(other.partial_charge, type(self.partial_charge))
+            and (
+                (other.partial_charge is None and self.partial_charge is None)
+                or np.isclose(other.partial_charge, self.partial_charge)
+            )
+            and np.allclose(other._coord, self._coord)
+        )
+        return are_equal
+
     @property
     def atomic_number(self) -> int:
         """

autode/calculations/calculation.py

@@ -217,6 +217,12 @@ def _check(self) -> None:
         if self.molecule.atoms is None or self.molecule.n_atoms == 0:
             raise ex.NoInputError("Have no atoms. Can't form a calculation")
 
+        if not self.molecule.has_valid_spin_state:
+            raise ex.CalculationException(
+                f"Cannot execute a calculation without a valid spin state: "
+                f"Spin multiplicity (2S+1) = {self.molecule.mult}"
+            )
+
         return None
 
     def _add_to_comp_methods(self) -> None:

autode/calculations/executors.py

@@ -199,7 +199,7 @@ def clean_up(self, force: bool = False, everything: bool = False) -> None:
 
         logger.info(f"Deleting: {set(filenames)}")
 
-        for filename in set(filenames):
+        for filename in [fn for fn in set(filenames) if fn is not None]:
 
             try:
                 os.remove(filename)

autode/conformers/conformer.py

@@ -43,7 +43,7 @@ def __init__(
         self.constraints.update(distance=dist_consts)
 
     def __repr__(self):
-        """Representation of"""
+        """Representation of a conformer"""
         return self._repr(prefix="Conformer")
 
     def __eq__(self, other):
@@ -158,7 +158,7 @@ def atoms(self, value: Optional[Atoms]):
 
         if value is None:  # Clear the coordinates
             self._coordinates = None
-            return None
+            return
 
         if self._parent_atoms is None:
             self._parent_atoms = value

autode/hessians.py

@@ -6,7 +6,7 @@
 import multiprocessing as mp
 
 from functools import cached_property
-from typing import List, Tuple, Iterator, Optional, Sequence, Union
+from typing import List, Tuple, Iterator, Optional, Sequence, Union, Any
 from autode.wrappers.keywords import Functional, GradientKeywords
 from autode.log import logger
 from autode.config import Config

autode/reactions/reaction.py

@@ -1,6 +1,9 @@
+import os
 import base64
 import hashlib
-from typing import Union, Optional, List
+import pickle
+
+from typing import Union, Optional, List, Generator
 from datetime import date
 from autode.config import Config
 from autode.solvent.solvents import get_solvent
@@ -13,7 +16,11 @@
 from autode.species.molecule import Reactant, Product
 from autode.plotting import plot_reaction_profile
 from autode.values import Energy, PotentialEnergy, Enthalpy, FreeEnergy
-from autode.utils import work_in, requires_hl_level_methods
+from autode.utils import (
+    work_in,
+    requires_hl_level_methods,
+    checkpoint_rxn_profile_step,
+)
 from autode.reactions import reaction_types
 
 
@@ -192,14 +199,16 @@ def _check_solvent(self) -> None:
         If self.solvent is set then override the reactants and products
         """
         molecules = self.reacs + self.prods
+        if len(molecules) == 0:
+            return  # No molecules thus no solvent needs to be checked
+
         first_solvent = self.reacs[0].solvent
 
         if self.solvent is None:
             if all([mol.solvent is None for mol in molecules]):
                 logger.info("Reaction is in the gas phase")
                 return
 
-            # Are solvents defined for all molecules?
             elif all([mol.solvent is not None for mol in molecules]):
 
                 if not all(
@@ -302,16 +311,21 @@ def _init_from_molecules(self, molecules) -> None:
 
         return None
 
-    def _reasonable_components_with_energy(self) -> None:
-        """Generator for components of a reaction that have sensible geometries
-        and also energies"""
+    def _components(self) -> Generator:
+        """Components of this reaction"""
 
         for mol in (
             self.reacs
             + self.prods
             + [self.ts, self._reactant_complex, self._product_complex]
         ):
+            yield mol
 
+    def _reasonable_components_with_energy(self) -> Generator:
+        """Generator for components of a reaction that have sensible geometries
+        and also energies"""
+
+        for mol in self._components():
             if mol is None:
                 continue
 
@@ -324,8 +338,6 @@ def _reasonable_components_with_energy(self) -> None:
 
             yield mol
 
-        return None
-
     def _estimated_barrierless_delta(self, e_type: str) -> Optional[Energy]:
         """
         Assume an effective free energy barrier = 4.35 kcal mol-1 calcd.
@@ -571,6 +583,7 @@ def switch_reactants_products(self) -> None:
         )
         return None
 
+    @checkpoint_rxn_profile_step("reactant_product_conformers")
     def find_lowest_energy_conformers(self) -> None:
         """Try and locate the lowest energy conformation using simulated
         annealing, then optimise them with xtb, then optimise the unique
@@ -584,6 +597,7 @@ def find_lowest_energy_conformers(self) -> None:
 
         return None
 
+    @checkpoint_rxn_profile_step("reactants_and_products")
     @work_in("reactants_and_products")
     def optimise_reacs_prods(self) -> None:
         """Perform a geometry optimisation on all the reactants and products
@@ -596,6 +610,7 @@ def optimise_reacs_prods(self) -> None:
 
         return None
 
+    @checkpoint_rxn_profile_step("complexes")
     @work_in("complexes")
     def calculate_complexes(self) -> None:
         """Find the lowest energy conformers of reactant and product complexes
@@ -618,24 +633,29 @@ def calculate_complexes(self) -> None:
         return None
 
     @requires_hl_level_methods
+    @checkpoint_rxn_profile_step("transition_states")
     @work_in("transition_states")
     def locate_transition_state(self) -> None:
 
+        if self.type is None:
+            raise RuntimeError(
+                "Cannot invoke locate_transition_state without a reaction type"
+            )
+
         # If there are more bonds in the product e.g. an addition reaction then
         # switch as the TS is then easier to find
         if sum(p.graph.number_of_edges() for p in self.prods) > sum(
             r.graph.number_of_edges() for r in self.reacs
         ):
-
             self.switch_reactants_products()
             self.tss = find_tss(self)
             self.switch_reactants_products()
-
         else:
             self.tss = find_tss(self)
 
         return None
 
+    @checkpoint_rxn_profile_step("transition_state_conformers")
     @work_in("transition_states")
     def find_lowest_energy_ts_conformer(self) -> None:
         """Find the lowest energy conformer of the transition state"""
@@ -646,6 +666,7 @@ def find_lowest_energy_ts_conformer(self) -> None:
         else:
             return self.ts.find_lowest_energy_ts_conformer()
 
+    @checkpoint_rxn_profile_step("single_points")
     @work_in("single_points")
     def calculate_single_points(self) -> None:
         """Perform a single point energy evaluations on all the reactants and
@@ -719,6 +740,7 @@ def print_energies_to_csv(_mol):
 
         return None
 
+    @checkpoint_rxn_profile_step("thermal")
     @work_in("thermal")
     def calculate_thermochemical_cont(
         self, free_energy: bool = True, enthalpy: bool = True
@@ -823,3 +845,24 @@ def atomic_symbols(self) -> List[str]:
     def has_identical_composition_as(self, reaction: "Reaction") -> bool:
         """Does this reaction have the same chemical identity as another?"""
         return self.atomic_symbols == reaction.atomic_symbols
+
+    def save(self, filepath: str) -> None:
+        """Save the state of this reaction to a binary file that can be reloaded"""
+
+        with open(filepath, "wb") as file:
+            pickle.dump(self.__dict__, file)
+
+    def load(self, filepath: str) -> None:
+        """Load a reaction state from a binary file"""
+
+        with open(filepath, "rb") as file:
+            for attr, value in dict(pickle.load(file)).items():
+                setattr(self, attr, value)
+
+    @classmethod
+    def from_checkpoint(cls, filepath: str) -> "Reaction":
+        """Create a reaction from a checkpoint file"""
+        logger.info(f"Loading a reaction object from {filepath}")
+        rxn = cls()
+        rxn.load(filepath)
+        return rxn

autode/reactions/reaction_types.py

@@ -1,4 +1,4 @@
-from typing import Sequence
+from typing import Sequence, Optional
 from autode.exceptions import ReactionFormationFailed
 from autode.log import logger
 
@@ -31,7 +31,7 @@ def __eq__(self, other) -> bool:
 def classify(
     reactants: Sequence["autode.species.molecule.Reactant"],
     products: Sequence["autode.species.molecule.Product"],
-) -> ReactionType:
+) -> Optional[ReactionType]:
     """
     Classify a reaction into a type given some reactants and products
 
@@ -45,6 +45,9 @@ def classify(
     """
     n_reactants, n_products = len(reactants), len(products)
 
+    if n_reactants == n_products == 0:
+        return None
+
     if n_reactants == 0:
         raise ReactionFormationFailed(
             f"Reaction had 0 reactants and "

autode/species/species.py

@@ -160,6 +160,15 @@ def mult(self) -> int:
 
     @mult.setter
     def mult(self, value: Any) -> None:
+
+        try:
+            assert int(value) > 0
+        except (ValueError, AssertionError, TypeError):
+            raise ValueError(
+                f"Failed to set the spin multiplicity to {value}. "
+                f"Must be a non-zero positive integer"
+            )
+
         self._mult = int(value)
 
     @property
@@ -774,6 +783,31 @@ def has_reasonable_coordinates(self) -> bool:
 
         return True
 
+    @property
+    def has_valid_spin_state(self) -> bool:
+        """
+        Does this species have a valid spin state given the atomic composition and
+        charge state?
+
+        .. code-block:: Python
+
+            >>> import autode as ade
+            >>> h = ade.Molecule(atoms=[ade.Atom('H')], charge=0, mult=1)
+            >>> h.has_valid_spin_state
+            False
+            >>> hydride = ade.Molecule(atoms=[ade.Atom('H')], charge=-1, mult=1)
+            >>> hydride.has_valid_spin_state
+            True
+        """
+        num_electrons = (
+            sum(atom.atomic_number for atom in self.atoms) - self.charge
+        )
+        num_unpaired_electrons = self.mult - 1
+        return (
+            num_unpaired_electrons <= num_electrons
+            and num_electrons % 2 == num_unpaired_electrons % 2
+        )
+
     @property
     def n_conformers(self) -> int:
         """