Throw exception on invalid molecule spin state (duartegroup#220)

* Add exception on calculation init

* Add test

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/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:
         """

tests/test_atoms.py

@@ -462,3 +462,10 @@ def test_atom_equality():
     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

tests/test_calculation.py

@@ -581,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
+        )

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,15 +134,15 @@ 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"),
     )
 
     product = Reactant(
         name="final",
         charge=-1,
-        mult=0,
+        mult=1,
         solvent_name="water",
         atoms=xyz_file_to_atoms("sn2_final.xyz"),
     )

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))
 

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

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

tests/test_ts/test_mode_checking.py

@@ -4,6 +4,7 @@
 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")
 

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),

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,

tests/test_wrappers/test_nwchem.py

@@ -122,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(

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,
     )

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(),
     )

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,