Repository: duartegroup/autodE
Pull request #225 (2023-01-11 19:04:08)

V1.3.4


Patch:
diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md
new file mode 100644
index 0000000..551bc72
--- /dev/null
+++ b/.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. -->
diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md
new file mode 100644
index 0000000..2db4fc7
--- /dev/null
+++ b/.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
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 5ca80c4..5ed79b8 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -4,4 +4,3 @@ repos:
     hooks:
       - id: black
         language_version: python3.9
-
diff --git a/autode/__init__.py b/autode/__init__.py
index eff4ba3..5789c25 100644
--- a/autode/__init__.py
+++ b/autode/__init__.py
@@ -7,6 +7,7 @@ from autode import mol_graphs
 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 @@ Bumping the version number requires following the release proceedure:
   - Merge when tests pass
 """
 
-__version__ = "1.3.3"
+__version__ = "1.3.4"
 
 
 __all__ = [
@@ -55,6 +56,7 @@ __all__ = [
     "Reactant",
     "Product",
     "Molecule",
+    "TransitionState",
     "NCIComplex",
     "Config",
     "Calculation",
diff --git a/autode/atoms.py b/autode/atoms.py
index 620b447..b441ce6 100644
--- a/autode/atoms.py
+++ b/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 @@ class Atom:
     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:
         """
diff --git a/autode/calculations/calculation.py b/autode/calculations/calculation.py
index fb4f9bd..9318f0e 100644
--- a/autode/calculations/calculation.py
+++ b/autode/calculations/calculation.py
@@ -217,6 +217,12 @@ class Calculation:
         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:
diff --git a/autode/calculations/executors.py b/autode/calculations/executors.py
index 9efb42c..a9b6c8e 100644
--- a/autode/calculations/executors.py
+++ b/autode/calculations/executors.py
@@ -199,7 +199,7 @@ class CalculationExecutor:
 
         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)
diff --git a/autode/conformers/conformer.py b/autode/conformers/conformer.py
index b32c9d8..ea5f39f 100644
--- a/autode/conformers/conformer.py
+++ b/autode/conformers/conformer.py
@@ -43,7 +43,7 @@ class Conformer(Species):
         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 @@ class Conformer(Species):
 
         if value is None:  # Clear the coordinates
             self._coordinates = None
-            return None
+            return
 
         if self._parent_atoms is None:
             self._parent_atoms = value
diff --git a/autode/hessians.py b/autode/hessians.py
index e86057c..b3f6149 100644
--- a/autode/hessians.py
+++ b/autode/hessians.py
@@ -6,7 +6,7 @@ import numpy as np
 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
diff --git a/autode/reactions/reaction.py b/autode/reactions/reaction.py
index 00d0b76..a828c13 100644
--- a/autode/reactions/reaction.py
+++ b/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.complex import ReactantComplex, ProductComplex
 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,6 +199,9 @@ class Reaction:
         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:
@@ -199,7 +209,6 @@ class Reaction:
                 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 @@ class Reaction:
 
         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 @@ class Reaction:
 
             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 @@ class Reaction:
         )
         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 @@ class Reaction:
 
         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 @@ class Reaction:
 
         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 @@ class Reaction:
         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 @@ class Reaction:
         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 @@ class Reaction:
 
         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 @@ class Reaction:
     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
diff --git a/autode/reactions/reaction_types.py b/autode/reactions/reaction_types.py
index 3b0161a..bd87c3b 100644
--- a/autode/reactions/reaction_types.py
+++ b/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 @@ Rearrangement = ReactionType(name="rearrangement")
 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 "
diff --git a/autode/species/species.py b/autode/species/species.py
index b4b3ee4..d630605 100644
--- a/autode/species/species.py
+++ b/autode/species/species.py
@@ -160,6 +160,15 @@ class Species(AtomCollection):
 
     @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 @@ class Species(AtomCollection):
 
         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:
         """
diff --git a/autode/transition_states/base.py b/autode/transition_states/base.py
index 6ce9b99..b36e111 100644
--- a/autode/transition_states/base.py
+++ b/autode/transition_states/base.py
@@ -4,6 +4,7 @@ from typing import Optional
 import autode.exceptions as ex
 from autode.atoms import metals
 from autode.config import Config
+from autode.geom import calc_rmsd
 from autode.constraints import DistanceConstraints
 from autode.log import logger
 from autode.methods import get_hmethod, get_lmethod
@@ -73,9 +74,12 @@ class TSbase(Species, ABC):
         self._init_solvent()
 
     def __eq__(self, other):
-        """Equality of this TS base to another"""
-        logger.warning("TSs types are not equal to any others")
-        return False
+        """Equality of this TS to another"""
+        return (
+            isinstance(other, TSbase)
+            and calc_rmsd(self.coordinates, other.coordinates) < 1e-6,
+            super().__eq__(other),
+        )
 
     def _init_graph(self) -> None:
         """Set the molecular graph for this TS object from the reactant"""
@@ -148,7 +152,7 @@ class TSbase(Species, ABC):
         """
         if self.bond_rearrangement is None:
             raise ValueError(
-                "Do not have a bond rearrangment - cannot "
+                "Do not have a bond rearrangement - cannot "
                 "check the imaginary mode"
             )
 
@@ -165,9 +169,7 @@ class TSbase(Species, ABC):
             )
             return False
 
-        if len(imag_freqs) > 1:
-            logger.warning(f"Hessian had {len(imag_freqs)} imaginary modes")
-
+        logger.info(f"Hessian had {len(imag_freqs)} imaginary modes")
         if imag_freqs[0] > Config.min_imag_freq:
             logger.warning("Imaginary modes were too small to be significant")
             return False
diff --git a/autode/transition_states/templates.py b/autode/transition_states/templates.py
index 3bea0e8..5bb77e1 100644
--- a/autode/transition_states/templates.py
+++ b/autode/transition_states/templates.py
@@ -38,6 +38,11 @@ def get_ts_template_folder_path(folder_path):
 
     logger.info("Folder path is not set – TS templates in the default path")
 
+    if Config.ts_template_folder_path == "":
+        raise ValueError(
+            "Cannot set ts_template_folder_path to an empty string"
+        )
+
     if Config.ts_template_folder_path is not None:
         logger.info("Configuration ts_template_folder_path is set")
         return Config.ts_template_folder_path
diff --git a/autode/transition_states/transition_state.py b/autode/transition_states/transition_state.py
index 80a17fc..1490c5a 100644
--- a/autode/transition_states/transition_state.py
+++ b/autode/transition_states/transition_state.py
@@ -3,6 +3,7 @@ from multiprocessing import Pool
 from autode.values import Frequency
 from autode.transition_states.base import displaced_species_along_mode
 from autode.transition_states.base import TSbase
+from autode.transition_states.ts_guess import TSguess
 from autode.transition_states.templates import TStemplate
 from autode.input_output import atoms_to_xyz_file
 from autode.calculations import Calculation
@@ -269,6 +270,7 @@ class TransitionState(TSbase):
         # Generate a copy of this TS on which conformers are searched, for
         # easy reversion
         _ts = self.copy()
+        _ts.hessian, _ts.gradient, _ts.energy = None, None, None
 
         hmethod = get_hmethod() if Config.hmethod_conformers else None
         _ts.find_lowest_energy_conformer(hmethod=hmethod)
@@ -358,3 +360,20 @@ class TransitionState(TSbase):
 
         logger.info("Saved TS template")
         return None
+
+    @classmethod
+    def from_species(
+        cls, species: "autode.species.Species"
+    ) -> "TransitionState":
+        """
+        Generate a TS from a species. Note this does not set the bond rearrangement
+        thus mode checking will not work from this species.
+
+        -----------------------------------------------------------------------
+        Arguments:
+            species:
+
+        Returns:
+            (autode.transition_states.transition_state.TransitionState): TS
+        """
+        return cls(ts_guess=TSguess.from_species(species), bond_rearr=None)
diff --git a/autode/utils.py b/autode/utils.py
index 74ff459..38517cb 100644
--- a/autode/utils.py
+++ b/autode/utils.py
@@ -135,7 +135,7 @@ def run_external_monitored(
     return None
 
 
-def work_in(dir_ext: str):
+def work_in(dir_ext: str) -> Callable:
     """Execute a function in a different directory"""
 
     def func_decorator(func):
@@ -173,7 +173,7 @@ def work_in_tmp_dir(
     filenames_to_copy: Optional[Sequence[str]] = None,
     kept_file_exts: Optional[Sequence[str]] = None,
     use_ll_tmp: bool = False,
-):
+) -> Callable:
     """Execute a function in a temporary directory.
 
     -----------------------------------------------------------------------
@@ -242,7 +242,7 @@ def work_in_tmp_dir(
     return func_decorator
 
 
-def log_time(prefix: str = "Executed in: ", units: str = "ms"):
+def log_time(prefix: str = "Executed in: ", units: str = "ms") -> Callable:
     """A function requiring a number of atoms to run"""
 
     if units.lower() == "s" or units.lower() == "seconds":
@@ -273,7 +273,7 @@ def log_time(prefix: str = "Executed in: ", units: str = "ms"):
     return func_decorator
 
 
-def requires_atoms(func: Callable):
+def requires_atoms(func: Callable) -> Callable:
     """A function requiring a number of atoms to run"""
 
     @wraps(func)
@@ -290,7 +290,7 @@ def requires_atoms(func: Callable):
     return wrapped_function
 
 
-def requires_graph(func: Callable):
+def requires_graph(func: Callable) -> Callable:
     """A function requiring a number of atoms to run"""
 
     @wraps(func)
@@ -307,7 +307,7 @@ def requires_graph(func: Callable):
     return wrapped_function
 
 
-def requires_conformers(func: Callable):
+def requires_conformers(func: Callable) -> Callable:
     """A function requiring the species to have a list of conformers"""
 
     @wraps(func)
@@ -324,7 +324,7 @@ def requires_conformers(func: Callable):
     return wrapped_function
 
 
-def requires_hl_level_methods(func: Callable):
+def requires_hl_level_methods(func: Callable) -> Callable:
     """A function requiring both high and low-level methods to be available"""
 
     @wraps(func)
@@ -352,7 +352,7 @@ def requires_hl_level_methods(func: Callable):
     return wrapped_function
 
 
-def requires_output(func: Callable):
+def requires_output(func: Callable) -> Callable:
     """A function requiring an output file and output file lines"""
 
     @wraps(func)
@@ -368,7 +368,7 @@ def requires_output(func: Callable):
     return wrapped_function
 
 
-def requires_output_to_exist(func):
+def requires_output_to_exist(func) -> Callable:
     """Calculation method requiring the output filename to be set"""
 
     @wraps(func)
@@ -456,7 +456,7 @@ def hashable(_method_name: str, _object: Any):
     return getattr(_object, _method_name)
 
 
-def run_in_tmp_environment(**kwargs):
+def run_in_tmp_environment(**kwargs) -> Callable:
     """
     Apply a set of environment variables, execute a function and reset them
     """
@@ -494,7 +494,7 @@ def run_in_tmp_environment(**kwargs):
     return func_decorator
 
 
-def deprecated(func):
+def deprecated(func: Callable) -> Callable:
     @wraps(func)
     def wrapped_function(*args, **kwargs):
         warnings.warn(
@@ -508,6 +508,44 @@ def deprecated(func):
     return wrapped_function
 
 
+def checkpoint_rxn_profile_step(name: str) -> Callable:
+    """
+    Decorator for a function that will save a checkpoint file with the reaction state
+    at that point in time. If the checkpoint exists then the state will be reloaded
+    and the execution skipped
+    """
+
+    def func_decorator(func: Callable[["Reaction"], Any]):
+        @wraps(func)
+        def wrapped_function(reaction: "Reaction"):
+
+            filepath = os.path.join(
+                "checkpoints", f"{str(reaction)}_{name}.chk"
+            )
+            if os.path.exists(filepath):
+                reaction.load(filepath)
+                return
+
+            start_time = time()
+            result = func(reaction)
+
+            if (
+                time() - start_time < 1.0
+            ):  # If execution is < 1s don't checkpoint
+                return result
+
+            if not os.path.exists("checkpoints"):
+                os.mkdir("checkpoints")
+
+            reaction.save(filepath)
+
+            return result
+
+        return wrapped_function
+
+    return func_decorator
+
+
 class StringDict:
     r"""
     Immutable dictionary stored as a single string. For example::
diff --git a/autode/values.py b/autode/values.py
index e611b2d..dcf4c09 100644
--- a/autode/values.py
+++ b/autode/values.py
@@ -1,6 +1,6 @@
 import numpy as np
 from abc import ABC, abstractmethod
-from typing import Any, Union, Type, Optional, Sequence
+from typing import Any, Union, Type, Optional, Sequence, Tuple
 from copy import deepcopy
 from collections.abc import Iterable
 from autode.log import logger
@@ -617,6 +617,18 @@ class ValueArray(ABC, np.ndarray):
 
         return arr
 
+    def __reduce__(self):
+        numpy_state = super().__reduce__()
+        return (
+            numpy_state[0],
+            numpy_state[1],
+            tuple(numpy_state[2]) + (self.__dict__,),
+        )
+
+    def __setstate__(self, state, *args, **kwargs):
+        self.__dict__.update(state[-1])
+        super().__setstate__(state[:-1], *args, **kwargs)
+
     def to(self, units) -> Any:
         """
         Convert this array to a new unit, returning a copy
diff --git a/autode/wrappers/methods.py b/autode/wrappers/methods.py
index 60941bc..f87ab0b 100644
--- a/autode/wrappers/methods.py
+++ b/autode/wrappers/methods.py
@@ -79,6 +79,11 @@ class Method(ABC):
             s.name for s in solvents if s.is_implicit and hasattr(s, self.name)
         ]
 
+    @property
+    def is_available(self):
+        """Is this method available?"""
+        return True
+
     def version_in(self, calc: "CalculationExecutor") -> str:
         """Determine the version of the method used in this calculation"""
         return "???"
diff --git a/doc/changelog.rst b/doc/changelog.rst
index 8376c45..43ab9f2 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -1,6 +1,30 @@
 Changelog
 =========
 
+1.3.4
+--------
+----------
+
+Feature additions.
+
+Usability improvements/Changes
+******************************
+* Throw useful exception for invalid :code:`ade.Config.ts_template_folder_path`
+* Adds the reaction temperature to the unique reaction hash
+
+
+Functionality improvements
+**************************
+- Adds :code:`ade.transition_states.TransitionState.from_species` method to construct transition states from a species or molecule
+- Adds :code:`autode.Reaction.save()` and :code:`autode.Reaction.load()` to save and reload a reaction state
+- Adds saving checkpoints of a reaction during :code:`autode.Reaction.calculate_reaction_profile`
+
+
+Bug Fixes
+*********
+-
+
+
 1.3.3
 --------
 ----------
diff --git a/doc/dev/contributing.rst b/doc/dev/contributing.rst
index 5a38ee0..6f39304 100644
--- a/doc/dev/contributing.rst
+++ b/doc/dev/contributing.rst
@@ -50,7 +50,7 @@ The high-level view of the contributing workflow is:
 2. Implement changes and tests on your own fork on a given branch
    (``<gh-username>/autode:<branch-name>``).
 3. Create a new pull request on the main autodE repository from your
-   development branch onto ``autode:master``.
+   development branch onto ``autode:v1.X``, where `X` is the latest version.
 
 To learn more about GitHub forks and pull requests, read `Fork a
 repo <https://docs.github.com/en/get-started/quickstart/fork-a-repo>`__
@@ -62,13 +62,13 @@ on the GitHub docs.
 Guidelines for pull requests
 ----------------------------
 
-First, install from source in a new envrionment and setup
+First, install from source in a new environment and setup
 `pre-commit <https://pre-commit.com/>`__ with::
 
     $ git clone https://github.com/duartegroup/autodE.git && cd autodE
     $ conda create -n ade python=3.9 --file requirements.txt --channel conda-forge
     $ conda activate ade
-    $ pip install . autode[dev]
+    $ pip install '.[dev]'
     $ pre-commit install
 
 
@@ -86,8 +86,8 @@ Several, smaller pull requests instead of one big PR
 Smaller pull requests (PRs) are reviewed faster, and more accurately. We
 therefore encourage contributors to keep the set of changes within a
 single pull request as small as possible. If your pull request modifies
-more than 5 files, and/or several hundred lines of code, you should
-probably break it down into two or more pull requests.
+more than 5 files, and/or several hundred lines of code, please break it down
+into two or more pull requests.
 
 
 Pull requests are more than code
diff --git a/examples/tutorials/t_identity_reactions.py b/examples/tutorials/t_identity_reactions.py
new file mode 100644
index 0000000..634e36a
--- /dev/null
+++ b/examples/tutorials/t_identity_reactions.py
@@ -0,0 +1,22 @@
+import autode as ade
+
+ade.Config.n_cores = 2
+ade.Config.hcode = "orca"
+
+if not ade.methods.ORCA().is_available:
+    exit("This example requires an ORCA install")
+
+# Identity reactions require the reactants to be somehow distinguished from the products
+# so that the breaking/forming bonds can be identified. For example, for
+# H-H + H -> H + H-H
+rxn = ade.Reaction(
+    ade.Reactant(atoms=[ade.Atom("H", atom_class=1)], mult=2),
+    ade.Reactant(atoms=[ade.Atom("H"), ade.Atom("H", x=0.8)]),
+    ade.Product(atoms=[ade.Atom("H")], mult=2),
+    ade.Product(atoms=[ade.Atom("H"), ade.Atom("H", atom_class=1, x=0.8)]),
+)
+rxn.calculate_reaction_profile()
+
+print("∆E_r (kcal mol-1) = ", rxn.delta("E").to("kcal mol-1"))
+print("∆E‡ (kcal mol-1)  = ", rxn.delta("E‡").to("kcal mol-1"))
+print("TS imaginary freq = ", rxn.ts.imaginary_frequencies[0])
diff --git a/examples/tutorials/u_reaction_profile_reload.py b/examples/tutorials/u_reaction_profile_reload.py
new file mode 100644
index 0000000..87cb3e2
--- /dev/null
+++ b/examples/tutorials/u_reaction_profile_reload.py
@@ -0,0 +1,23 @@
+import autode as ade
+
+ade.Config.n_cores = 4
+ade.Config.num_conformers = 1  # Use only a single conformer for speed
+ade.Config.hcode = "orca"
+
+if not ade.methods.ORCA().is_available:
+    exit("This example requires an ORCA install")
+
+# Calculate a simple reaction profile
+rxn = ade.Reaction("[F-].CCl>>CF.[Cl-]", solvent_name="water", temp=300)
+rxn.calculate_reaction_profile(free_energy=True)
+rxn.save("sn2_reaction.chk")
+
+print("∆G‡(300 K) = ", rxn.delta("G‡").to("kcal mol-1"))
+
+# Reactions can be reloaded from checkpoints and e.g. the thermal cont. recalculated
+# without recalculating any energies using external methods
+reloaded_rxn = ade.Reaction.from_checkpoint("sn2_reaction.chk")
+reloaded_rxn.temp = 400  # Set the new temperature to 400 K
+reloaded_rxn.calculate_thermochemical_cont()
+
+print("∆G‡(400 K) = ", reloaded_rxn.delta("G‡").to("kcal mol-1"))
diff --git a/setup.py b/setup.py
index b8206c6..ddfe8d1 100644
--- a/setup.py
+++ b/setup.py
@@ -24,7 +24,7 @@ extensions = [
 
 setup(
     name="autode",
-    version="1.3.3",
+    version="1.3.4",
     python_requires=">3.7",
     packages=[
         "autode",


Test patch:
diff --git a/tests/test_atoms.py b/tests/test_atoms.py
index df871b4..9c369be 100644
--- a/tests/test_atoms.py
+++ b/tests/test_atoms.py
@@ -316,11 +316,12 @@ def test_atom_coord_setting():
 
 def test_periodic_table():
 
-    with pytest.raises(ValueError):
-        _ = atoms.PeriodicTable.period(n=0)  # Periods are indexed from 1
-        _ = atoms.PeriodicTable.period(n=8)  # and don't exceed 8
+    for invalid_period in (0, 8):
+        with pytest.raises(ValueError):
+            _ = atoms.PeriodicTable.period(n=invalid_period)
 
-        _ = atoms.PeriodicTable.period(n=19)  # Groups don't exceed 18
+    with pytest.raises(ValueError):
+        _ = atoms.PeriodicTable.group(n=19)  # Groups don't exceed 18
 
     period2 = atoms.PeriodicTable.period(n=2)
     assert len(period2) == 8
@@ -452,3 +453,19 @@ def test_eqm_bond_distance():
     # Cannot determine the distance between atoms not present in the set
     with pytest.raises(ValueError):
         _ = h2_atoms.eqm_bond_distance(0, 2)
+
+
+def test_atom_equality():
+
+    assert Atom("H") == Atom("H")
+    assert Atom("H") != Atom("H", partial_charge=0.1)
+    assert Atom("H", partial_charge=0.1) != Atom("H")
+    assert Atom("H", atom_class=1) != Atom("H", atom_class=0)
+    assert Atom("C") != Atom("H")
+
+
+@pytest.mark.parametrize(
+    "element,atomic_number", [("H", 1), ("C", 6), ("F", 9), ("Cl", 17)]
+)
+def test_atomic_numbers(element: str, atomic_number: int):
+    assert Atom(element).atomic_number == atomic_number
diff --git a/tests/test_calculation.py b/tests/test_calculation.py
index 6540da9..06e1399 100644
--- a/tests/test_calculation.py
+++ b/tests/test_calculation.py
@@ -302,10 +302,10 @@ def test_exec_not_avail_method():
 
 
 @work_in_tmp_dir()
-def test_exec_too_much_memory_requested_py39():
+def test_exec_too_much_memory_requested_above_py39():
 
-    if sys.version_info.minor != 9:
-        return  # Only supported on Python 3.9
+    if sys.version_info.minor < 9:
+        return  # Only supported on Python 3.9 and above
 
     # Normal external run should be fine
     run_external(["ls"], output_filename="tmp.txt")
@@ -320,16 +320,6 @@ def test_exec_too_much_memory_requested_py39():
     Config.max_core = curr_max_core
 
 
-@work_in_tmp_dir()
-def test_exec_too_much_memory_requested_py38():
-
-    if sys.version_info.minor != 8:
-        return  # Only supported on Python 3.8
-
-    # Python 3.8 can't use os.sysconf, so check that external can still be run
-    run_external(["ls"], output_filename="tmp.txt")
-
-
 @requires_with_working_xtb_install
 @work_in_tmp_dir()
 def test_calculations_have_unique_names():
@@ -591,3 +581,14 @@ def test_non_external_io_method_can_force_cleanup():
     calc = _test_calc()
     # No exceptions should be raised
     calc.clean_up(force=True, everything=True)
+
+
+def test_init_a_calculation_without_a_valid_spin_state_throws():
+    xtb = XTB()
+    test_m = h_atom()
+    test_m.charge, test_m.mult = 0, 1
+
+    with pytest.raises(ex.CalculationException):
+        _ = Calculation(
+            name="tmp", molecule=test_m, method=xtb, keywords=xtb.keywords.sp
+        )
diff --git a/tests/test_config.py b/tests/test_config.py
index a0791d2..338a22a 100644
--- a/tests/test_config.py
+++ b/tests/test_config.py
@@ -5,6 +5,7 @@ from autode.config import Config
 from autode.values import Allocation, Distance
 from autode.wrappers.keywords import KeywordsSet
 from autode.wrappers.keywords import Keywords
+from autode.transition_states.templates import get_ts_template_folder_path
 
 
 def test_config():
@@ -76,3 +77,13 @@ def test_step_size_setter():
     # Setting in Bohr should convert to angstroms
     _config.max_step_size = Distance(0.2, units="a0")
     assert np.isclose(_config.max_step_size.to("ang"), 0.1, atol=0.02)
+
+
+def test_invalid_get_ts_template_folder_path():
+
+    Config.ts_template_folder_path = ""
+
+    with pytest.raises(ValueError):
+        _ = get_ts_template_folder_path(None)
+
+    Config.ts_template_folder_path = None
diff --git a/tests/test_hessian.py b/tests/test_hessian.py
index 1e9aa08..c611c22 100644
--- a/tests/test_hessian.py
+++ b/tests/test_hessian.py
@@ -1,5 +1,6 @@
 import os
 import pytest
+import pickle
 import numpy as np
 import autode as ade
 from autode.utils import work_in_tmp_dir
@@ -919,3 +920,19 @@ def test_serial_calculation_matches_parallel():
     serial_result = nhc.hessian
 
     assert np.allclose(parallel_result, serial_result)
+
+
+@work_in_tmp_dir()
+def test_hessian_pickle_and_unpickle():
+
+    mol = Molecule(smiles="O")
+    mol.hessian = Hessian(np.eye(3 * mol.n_atoms), atoms=mol.atoms)
+
+    with open("tmp.obj", "wb") as file:
+        pickle.dump(mol.hessian, file=file)
+
+    with open("tmp.obj", "rb") as file:
+        reloaded_hessian = pickle.load(file=file)
+
+    assert reloaded_hessian.shape == (3 * mol.n_atoms, 3 * mol.n_atoms)
+    assert reloaded_hessian.atoms == mol.atoms
diff --git a/tests/test_input_output.py b/tests/test_input_output.py
index ab4daa7..13bcaa2 100644
--- a/tests/test_input_output.py
+++ b/tests/test_input_output.py
@@ -89,6 +89,7 @@ def test_making_xyz_file():
         # Requires some atoms
         atoms_to_xyz_file(atoms=None, filename="test.xyz")
 
+    with pytest.raises(AssertionError):
         # Needs .xyz extension
         atoms_to_xyz_file(atoms, filename="test")
 
diff --git a/tests/test_log.py b/tests/test_log.py
index 5d44d41..e33bedc 100644
--- a/tests/test_log.py
+++ b/tests/test_log.py
@@ -4,9 +4,10 @@ import os
 
 def test_log_level():
 
-    # Default is debug
     if "AUTODE_LOG_LEVEL" in os.environ:
-        os.environ.pop("AUTODE_LOG_LEVEL")
+        set_log_level = os.environ.pop("AUTODE_LOG_LEVEL")
+    else:
+        set_log_level = "INFO"
 
     assert log.get_log_level() == log.logging.CRITICAL
 
@@ -22,3 +23,5 @@ def test_log_level():
     # Setting AUTODE_LOG_FILE to anything should log to a log file
     os.environ["AUTODE_LOG_FILE"] = "true"
     assert log.log_to_log_file() is True
+
+    os.environ["AUTODE_LOG_LEVEL"] = set_log_level
diff --git a/tests/test_neb.py b/tests/test_neb.py
index a7036bc..4019664 100644
--- a/tests/test_neb.py
+++ b/tests/test_neb.py
@@ -94,14 +94,14 @@ def test_full_calc_with_xtb():
         initial_species=Species(
             name="inital",
             charge=-1,
-            mult=0,
+            mult=1,
             atoms=xyz_file_to_atoms("sn2_init.xyz"),
             solvent_name="water",
         ),
         final_species=Species(
             name="final",
             charge=-1,
-            mult=0,
+            mult=1,
             atoms=xyz_file_to_atoms("sn2_final.xyz"),
             solvent_name="water",
         ),
@@ -134,7 +134,7 @@ def test_get_ts_guess_neb():
     reactant = Reactant(
         name="inital",
         charge=-1,
-        mult=0,
+        mult=1,
         solvent_name="water",
         atoms=xyz_file_to_atoms("sn2_init.xyz"),
     )
@@ -142,7 +142,7 @@ def test_get_ts_guess_neb():
     product = Reactant(
         name="final",
         charge=-1,
-        mult=0,
+        mult=1,
         solvent_name="water",
         atoms=xyz_file_to_atoms("sn2_final.xyz"),
     )
diff --git a/tests/test_opt/test_bfgs.py b/tests/test_opt/test_bfgs.py
index 1879d70..443ee9f 100644
--- a/tests/test_opt/test_bfgs.py
+++ b/tests/test_opt/test_bfgs.py
@@ -86,6 +86,9 @@ class TestBFGSOptimiser(TestBFGSOptimiser2D):
 def test_simple_quadratic_opt():
 
     optimiser = TestBFGSOptimiser()
+    method = Method()
+    assert method.is_available
+
     optimiser.run(Molecule(name="blank"), method=Method())
     assert optimiser.converged
 
diff --git a/tests/test_opt/test_dimer.py b/tests/test_opt/test_dimer.py
index c246925..9ec823f 100644
--- a/tests/test_opt/test_dimer.py
+++ b/tests/test_opt/test_dimer.py
@@ -33,7 +33,7 @@ def test_dimer_coord_init():
 def test_dimer_coord_mol_init():
 
     mol1 = Molecule()
-    mol2 = Molecule(atoms=[Atom("H")])
+    mol2 = Molecule(atoms=[Atom("H")], mult=2)
 
     # Dimer coordinates must be created from two species with the same
     # atomic composition
@@ -112,7 +112,7 @@ def test_mass_weighting_no_masses():
 
 def test_dimer_init_zero_distance():
 
-    a = Molecule(atoms=[Atom("H")])
+    a = Molecule(atoms=[Atom("H")], mult=2)
 
     dimer = Dimer(maxiter=10, coords=DimerCoordinates.from_species(a, a))
 
diff --git a/tests/test_plotting.py b/tests/test_plotting.py
index 1e0f5e7..a566764 100644
--- a/tests/test_plotting.py
+++ b/tests/test_plotting.py
@@ -145,6 +145,7 @@ def test_reaction_warnings():
     # Should be no warnings  with a TS that exists and has an energy and one
     # imaginary freq
     ts.atoms = xyz_file_to_atoms("TS.xyz")
+    ts.charge = -1
     orca = ORCA()
     ts_calc = Calculation(
         name="TS", molecule=ts, method=orca, keywords=orca.keywords.opt_ts
diff --git a/tests/test_reaction_class.py b/tests/test_reaction_class.py
index 560e9a6..8149df1 100644
--- a/tests/test_reaction_class.py
+++ b/tests/test_reaction_class.py
@@ -1,5 +1,6 @@
 import os
 import numpy as np
+from time import time
 from autode.reactions import reaction
 from autode.reactions import reaction_types
 from autode.transition_states.transition_state import TransitionState
@@ -345,7 +346,10 @@ def test_free_energy_profile():
         solvent_name="water",
     )
 
+    start_time = time()
     rxn.calculate_reaction_profile(free_energy=True)
+    full_calc_reaction_profile_time = time() - start_time
+    rxn.save("tmp.chk")
 
     # Allow ~0.5 kcal mol-1 either side of the 'true' value
 
@@ -360,6 +364,15 @@ def test_free_energy_profile():
     assert os.path.exists("enthalpy_reaction_profile.pdf")
     os.remove("enthalpy_reaction_profile.pdf")
 
+    # Rerunning the reaction should be fast
+    start_time = time()
+    rxn.calculate_reaction_profile(free_energy=True)
+    assert time() - start_time < full_calc_reaction_profile_time / 2
+
+    # Should be able to reload the entire reaction state
+    reloaded_rxn = reaction.Reaction.from_checkpoint("tmp.chk")
+    assert reloaded_rxn.ts is not None
+
     # Reset the configuration to the default values
     Config.hcode = None
     Config.G09.path = None
@@ -520,3 +533,15 @@ def test_identity_reaction_is_supported_with_labels():
 
     rxn = reaction.Reaction("[Br-:1].C[Br:2]>>C[Br:1].[Br-:2]")
     assert not reaction_is_isomorphic(rxn)
+
+
+def test_cannot_run_locate_ts_with_no_reactants_or_products():
+
+    Config.lcode = Config.hcode = "ORCA"
+    Config.ORCA.path = here
+
+    rxn = reaction.Reaction()
+    with pytest.raises(RuntimeError):
+        rxn.locate_transition_state()
+
+    Config.lcode = None
diff --git a/tests/test_reactions.py b/tests/test_reactions.py
index b493f5b..2219cae 100644
--- a/tests/test_reactions.py
+++ b/tests/test_reactions.py
@@ -24,10 +24,9 @@ def test_classify():
     rearrangement = classify([0], [0])
     assert rearrangement.name == "rearrangement"
 
-    # Needs to have at least some reactants and products
-    with pytest.raises(ReactionFormationFailed):
-        _ = classify([], [])
+    assert classify([], []) is None
 
+    # Needs to have at least some reactants or products
     with pytest.raises(ReactionFormationFailed):
         _ = classify([0], [])
 
diff --git a/tests/test_species.py b/tests/test_species.py
index 23eba98..87bbce6 100644
--- a/tests/test_species.py
+++ b/tests/test_species.py
@@ -671,3 +671,33 @@ def test_calc_thermo_not_run_calculation():
     # run() has not been called
     with pytest.raises(Exception):
         m.calc_thermo(calc=calc)
+
+
+@pytest.mark.parametrize("mult", [1, 3, 5])
+def test_argon_has_valid_spin_state(mult: int, charge: int = 0):
+    assert Molecule(
+        atoms=[Atom("Ar")], mult=mult, charge=charge
+    ).has_valid_spin_state
+
+
+@pytest.mark.parametrize("mult", [1, 3, 4])
+def test_hydrogen_has_invalid_spin_state(mult: int, charge: int = 0):
+    assert not Molecule(
+        atoms=[Atom("H")], mult=mult, charge=charge
+    ).has_valid_spin_state
+
+
+def test_has_valid_spin_state_docstring():
+
+    assert not Molecule(
+        atoms=[Atom("H")], charge=0, mult=1
+    ).has_valid_spin_state
+    assert Molecule(atoms=[Atom("H")], charge=-1, mult=1).has_valid_spin_state
+
+
+@pytest.mark.parametrize("invalid_mult", [0, -1, "a", (0, 2)])
+def test_cannot_set_multiplicity_to_invalid_value(invalid_mult):
+
+    m = Species(name="H2", atoms=[h1, h2], charge=0, mult=1)
+    with pytest.raises(Exception):
+        m.mult = invalid_mult
diff --git a/tests/test_ts/test_mode_checking.py b/tests/test_ts/test_mode_checking.py
index f4cd9e7..3089964 100644
--- a/tests/test_ts/test_mode_checking.py
+++ b/tests/test_ts/test_mode_checking.py
@@ -4,6 +4,7 @@ from autode.transition_states.base import imag_mode_generates_other_bonds
 from autode.transition_states.base import displaced_species_along_mode
 from autode.species.molecule import Reactant
 from autode.input_output import xyz_file_to_atoms
+from autode.atoms import Atom
 from autode.bond_rearrangement import BondRearrangement
 from autode.methods import ORCA
 from .. import testutils
@@ -45,7 +46,9 @@ def test_imag_modes():
 def test_graph_no_other_bonds():
 
     reac = Reactant(
-        name="r", atoms=xyz_file_to_atoms("h_shift_correct_ts_mode.xyz")
+        name="r",
+        atoms=xyz_file_to_atoms("h_shift_correct_ts_mode.xyz"),
+        mult=2,
     )
 
     calc = Calculation(
@@ -75,15 +78,16 @@ def test_graph_no_other_bonds():
 
 def has_correct_mode(name, fbonds, bbonds):
 
-    reac = Reactant(name="r", atoms=xyz_file_to_atoms(f"{name}.xyz"))
-
     calc = Calculation(
         name=name,
-        molecule=reac,
+        molecule=Reactant(atoms=[Atom("H")], mult=2),
         method=orca,
         keywords=orca.keywords.opt_ts,
         n_cores=1,
     )
+    calc.molecule = reac = Reactant(
+        name="r", atoms=xyz_file_to_atoms(f"{name}.xyz")
+    )
 
     calc.set_output_filename(f"{name}.out")
 
diff --git a/tests/test_ts/test_ts_base.py b/tests/test_ts/test_ts_base.py
new file mode 100644
index 0000000..5a48cdb
--- /dev/null
+++ b/tests/test_ts/test_ts_base.py
@@ -0,0 +1,58 @@
+import pytest
+
+from autode.config import Config
+from autode.atoms import Atom
+from autode.values import Frequency
+from autode.transition_states.ts_guess import TSguess
+from autode.species.molecule import Reactant, Product
+from autode.bond_rearrangement import BondRearrangement
+
+
+def h3_ts_guess():
+    return TSguess(atoms=[Atom("H"), Atom("H", x=0.7), Atom("H", 1.4)], mult=2)
+
+
+def test_invalid_init_different_reactants_products():
+
+    with pytest.raises(ValueError):
+        _ = TSguess(
+            atoms=[Atom("H")],
+            reactant=Reactant(atoms=[Atom("H")], solvent_name="water"),
+            product=Product(atoms=[Atom("H")], solvent_name="dcm"),
+        )
+
+
+def test_invalid_init_different_reactant():
+
+    with pytest.raises(ValueError):
+        _ = TSguess(
+            atoms=[Atom("H")],
+            reactant=Reactant(atoms=[Atom("H")], solvent_name="water"),
+            solvent_name="dcm",
+        )
+
+
+def test_mode_checking_without_bond_rearr():
+
+    ts_guess = h3_ts_guess()
+
+    with pytest.raises(Exception):
+        _ = ts_guess.could_have_correct_imag_mode
+
+
+def test_could_have_correct_imag_mode():
+    class TmpHess:
+        def __init__(self, freqs: list):
+            self.frequencies_proj = freqs
+
+    ts_guess = h3_ts_guess()
+    ts_guess.bond_rearrangement = BondRearrangement(
+        forming_bonds=[(0, 1)], breaking_bonds=[(1, 2)]
+    )
+
+    ts_guess._hess = TmpHess(freqs=[])  # No frequencies
+    assert ts_guess.imaginary_frequencies is None
+    assert not ts_guess.could_have_correct_imag_mode
+
+    ts_guess._hess = TmpHess(freqs=[Config.min_imag_freq / 2, Frequency(1000)])
+    assert not ts_guess.could_have_correct_imag_mode
diff --git a/tests/test_ts/test_ts_template.py b/tests/test_ts/test_ts_template.py
index ede99ad..7c9b5f0 100644
--- a/tests/test_ts/test_ts_template.py
+++ b/tests/test_ts/test_ts_template.py
@@ -5,10 +5,12 @@ from .. import testutils
 import pytest
 from autode.exceptions import TemplateLoadingFailed
 from autode.config import Config
+from autode.species.molecule import Molecule
 from autode.bond_rearrangement import BondRearrangement
 from autode.species.complex import ReactantComplex, ProductComplex
 from autode.species.molecule import Reactant, Product
 from autode.atoms import Atom
+from autode.utils import work_in_tmp_dir
 from autode.transition_states.templates import get_ts_templates
 from autode.transition_states.templates import get_value_from_file
 from autode.transition_states.templates import get_values_dict_from_file
@@ -214,6 +216,7 @@ def test_ts_template_parse():
         _ = get_values_dict_from_file("nodes", file_lines=["0 C", "1 F"])
 
 
+@work_in_tmp_dir()
 def test_ts_templates_find():
 
     templates = get_ts_templates(folder_path="/a/path/that/doesnt/exist")
@@ -238,3 +241,13 @@ def test_inactive_graph():
 
     template.graph = ch3f.graph.copy()
     assert not template.graph_has_correct_structure()
+
+
+@testutils.work_in_zipped_dir(os.path.join(here, "data", "ts_template.zip"))
+def test_ts_from_species_is_same_as_from_ts_guess():
+
+    ts = TransitionState(
+        TSguess(atoms=xyz_file_to_atoms("vaskas_TS.xyz"), charge=0, mult=1)
+    )
+
+    assert TransitionState.from_species(Molecule("vaskas_TS.xyz")) == ts
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 0ce056e..7df31fe 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -1,3 +1,4 @@
+import shutil
 import time
 import pytest
 import os
@@ -12,6 +13,7 @@ import multiprocessing as mp
 from autode import exceptions as ex
 from autode.mol_graphs import is_isomorphic
 from autode.utils import work_in_tmp_dir, log_time, requires_graph
+from .testutils import requires_with_working_xtb_install
 
 
 here = os.path.dirname(os.path.abspath(__file__))
@@ -244,10 +246,7 @@ def test_spawn_multiprocessing_graph():
 def test_time_units():
 
     with pytest.raises(ValueError):
-
-        @log_time(units="X")  # invalid time format
-        def f():
-            return None
+        log_time(units="X")  # invalid time format
 
 
 def test_requires_graph():
@@ -286,3 +285,19 @@ def test_string_dict():
     assert "a" in str(d)
     assert "solvent" in str(d)
     assert d["a"] == "b"
+
+
+def test_requires_xtb_install():
+
+    path_env_var = os.environ.pop("PATH")
+    assert shutil.which("xtb") is None
+
+    test_list = []
+
+    @requires_with_working_xtb_install
+    def tmp_function():
+        test_list.append("executed")
+
+    # If XTB is not in $PATH then the function should not execute
+    assert len(test_list) == 0
+    os.environ["PATH"] = path_env_var
diff --git a/tests/test_wrappers/test_gaussian.py b/tests/test_wrappers/test_gaussian.py
index c69f83c..d31da8c 100644
--- a/tests/test_wrappers/test_gaussian.py
+++ b/tests/test_wrappers/test_gaussian.py
@@ -100,7 +100,6 @@ def test_get_gradients():
 
     ester = Molecule(
         name="ester",
-        charge=-1,
         atoms=[
             Atom("C", -1.82707, 0.08502, 0.12799),
             Atom("C", -0.42971, 0.07495, -0.39721),
diff --git a/tests/test_wrappers/test_mopac.py b/tests/test_wrappers/test_mopac.py
index 44fe3fc..8f12652 100644
--- a/tests/test_wrappers/test_mopac.py
+++ b/tests/test_wrappers/test_mopac.py
@@ -1,6 +1,10 @@
 from autode.wrappers.MOPAC import MOPAC
 from autode.wrappers.MOPAC import get_keywords, _get_atoms_linear_interp
-from autode.exceptions import CouldNotGetProperty, UnsupportedCalculationInput
+from autode.exceptions import (
+    CouldNotGetProperty,
+    UnsupportedCalculationInput,
+    CalculationException,
+)
 from autode.calculations import Calculation, CalculationInput
 from autode.species.molecule import Molecule
 from autode.solvent import ImplicitSolvent
@@ -120,7 +124,7 @@ def test_other_spin_states():
     h_quin = Molecule(atoms=[Atom("H")], mult=5)
     h_quin.name = "molecule"
 
-    with pytest.raises(UnsupportedCalculationInput):
+    with pytest.raises(CalculationException):
         calc = Calculation(
             name="h",
             molecule=h_quin,
diff --git a/tests/test_wrappers/test_nwchem.py b/tests/test_wrappers/test_nwchem.py
index abbba0b..0e1bfb2 100644
--- a/tests/test_wrappers/test_nwchem.py
+++ b/tests/test_wrappers/test_nwchem.py
@@ -30,7 +30,6 @@ opt_keywords = OptKeywords(
 
 
 @testutils.work_in_zipped_dir(os.path.join(here, "data", "nwchem.zip"))
-# @utils.work_in("tmp")
 def test_opt_calc():
 
     calc = Calculation(
@@ -123,6 +122,7 @@ def test_hessian_extract_ts():
             Atom("H", -0.74001, 0.78538, 0.60979),
             Atom("H", -0.31016, -1.03356, 0.60979),
         ],
+        charge=-1,
     )
 
     calc = Calculation(
diff --git a/tests/test_wrappers/test_orca.py b/tests/test_wrappers/test_orca.py
index 773caaf..3c141ad 100644
--- a/tests/test_wrappers/test_orca.py
+++ b/tests/test_wrappers/test_orca.py
@@ -86,10 +86,9 @@ def test_orca_opt_calculation():
 @testutils.work_in_zipped_dir(os.path.join(here, "data", "orca.zip"))
 def test_orca_optts_calculation():
 
-    ts_guess = TSguess(
-        Molecule("test_ts_reopt_optts_orca.xyz", charge=-1).atoms
+    ts = TransitionState.from_species(
+        Molecule("test_ts_reopt_optts_orca.xyz", charge=-1)
     )
-    ts = TransitionState(ts_guess)
     ts.graph.add_active_edge(0, 1)
 
     optts_str = (
@@ -258,7 +257,7 @@ def test_mp2_numerical_gradients():
 
     calc = Calculation(
         name="tmp",
-        molecule=Molecule(atoms=xyz_file_to_atoms("tmp_orca.xyz")),
+        molecule=Molecule(atoms=xyz_file_to_atoms("tmp_orca.xyz"), charge=-1),
         method=method,
         keywords=method.keywords.grad,
     )
diff --git a/tests/test_wrappers/test_qchem.py b/tests/test_wrappers/test_qchem.py
index f931ba2..2356a46 100644
--- a/tests/test_wrappers/test_qchem.py
+++ b/tests/test_wrappers/test_qchem.py
@@ -24,7 +24,7 @@ def _blank_calc(name="test"):
 
     calc = Calculation(
         name=name,
-        molecule=Molecule(atoms=[Atom("H")]),
+        molecule=Molecule(atoms=[Atom("H")], mult=2),
         method=method,
         keywords=SinglePointKeywords(),
     )
diff --git a/tests/test_wrappers/test_xtb.py b/tests/test_wrappers/test_xtb.py
index 7d39a97..5430a7f 100644
--- a/tests/test_wrappers/test_xtb.py
+++ b/tests/test_wrappers/test_xtb.py
@@ -294,7 +294,7 @@ def test_xtb_with_autode_opt_method_for_a_single_atom():
 @work_in_tmp_dir()
 def test_xtb_opt_non_contiguous_range_cart_constraints():
 
-    mol = Molecule(smiles="CC", mult=2)
+    mol = Molecule(smiles="CC")
     mol.constraints.cartesian = [0, 1, 2, 5]
 
     calc = Calculation(
@@ -326,7 +326,7 @@ def test_xtb_errors_with_infinite_nuclear_repulsion():
 @work_in_tmp_dir()
 def test_xtb_did_not_terminate_normally_with_blank_output():
 
-    mol = Molecule(atoms=[Atom("H")])
+    mol = Molecule(atoms=[Atom("H")], mult=2)
     calc = Calculation(
         name="h_atom",
         molecule=mol,