Thermochemistry calculation exception (duartegroup#166)

* Throw useful exception in calc_thermo with non run calc

* Tweak sim axis tolerance for correct methane sn

autode/species/species.py

@@ -1149,6 +1149,7 @@ def calc_thermo(self,
             :meth:`autode.thermochemistry.igm.calculate_thermo_cont` for
             additional kwargs
         """
+        logger.info(f"Calculating thermochemical contributions for {self.name}")
 
         if 'lfm_method' in kwargs:
             try:
@@ -1158,8 +1159,14 @@ def calc_thermo(self,
                                  f'be one of: {[m for m in LFMethod]}')
 
         if calc is not None and calc.output.exists:
-            self.energy = calc.molecule.energy
+            logger.info("Setting the atoms, energy and Hessian from an "
+                        "existing calculation")
+            if calc.molecule.hessian is None:
+                raise ValueError(f"Failed to set the Hessian from {calc.name}."
+                                 f" Maybe run() hasn't been called?")
+
             self.atoms = calc.molecule.atoms.copy()
+            self.energy = calc.molecule.energy
             self.hessian = calc.molecule.hessian
 
         elif self.hessian is None or (calc is not None and not calc.output.exists):

autode/thermochemistry/symmetry.py

@@ -32,7 +32,7 @@ def strip_identical_and_inv_axes(axes, sim_axis_tol):
     return unique_possible_axes
 
 
-def get_possible_axes(coords, max_triple_dist=2.0, sim_axis_tol=0.05):
+def get_possible_axes(coords, max_triple_dist=2.0, sim_axis_tol=0.1):
     r"""
     Possible rotation axes in a molecule. Currently limited to average vectors
     and cross products i.e.::
@@ -215,7 +215,7 @@ def create_pcoords(species):
 def symmetry_number(species, max_n_fold_rot_searched=6,  dist_tol=0.25):
     """
     Calculate the symmetry number of a molecule. See:
-    Theor Chem Account (2007) 118:813–826
+    Theor Chem Account (2007) 118:813–826. 10.1007/s00214-007-0328-0
 
     ---------------------------------------------------------------------------
     Arguments:

doc/changelog.rst

@@ -47,6 +47,8 @@ Bug Fixes
 - Fixes variable harmonic frequencies (<2 cm-1 differences) due to projection vectors becoming close to rotational axes
 - Fixes the extraction of atomic partial charges from ORCA output files
 - Fixes gradients and Hessians not being reset on a molecule where the coordinates change
+- Fixes unhelpful exception when calculating thermochemistry with EST methods without implemented "get_hessian" methods
+
 
 See the table below for a quick benchmark of constrained optimisations in autodE
 compared to ORCA. In all cases the structures were generated from SMILES strings (RDKit)

tests/test_species.py

@@ -599,3 +599,14 @@ def test_species_does_not_have_reasonable_coordinates():
 
     x = ch4_flat.coordinates
     assert np.min(distance_matrix(x, x)+np.eye(5)) > 0.7
+
+
+@testutils.requires_with_working_xtb_install
+def test_calc_thermo_not_run_calculation():
+
+    m = Molecule(smiles="O")
+    calc = Calculation(name="water", molecule=m, method=xtb,
+                       keywords=xtb.keywords.hess)
+    # run() has not been called
+    with pytest.raises(Exception):
+        m.calc_thermo(calc=calc)

tests/test_thermochem.py

@@ -28,6 +28,7 @@ def test_symmetry_number():
     assert Molecule('CO2.xyz').symmetry_number == 2
     assert Molecule('H2O.xyz').symmetry_number == 2
     assert Molecule('H3N.xyz').symmetry_number == 3
+    assert Molecule(smiles="C").symmetry_number == 12
 
     # Symmetry numbers aren't calculated for large molecules
     h_100 = Species('tmp', atoms=100*[Atom('H')], charge=1, mult=1)