Repository: duartegroup/autodE
Pull request #351 (2024-08-02 15:12:41)

More comprehensive optimiser convergence checks
Optimiser convergence checks now use ORCA/Gaussian style parameters (energy change, RMS and max gradient and step sizes).

Partial fix of #300 

***

## Checklist

* [x] The changes include an associated explanation of how/why
* [x] Test pass
* [x] Documentation has been updated
* [x] Changelog has been updated


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/duartegroup/autodE/pull/351?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) Report
Attention: Patch coverage is `82.21154%` with `37 lines` in your changes missing coverage. Please review.
> Project coverage is 88.77%. Comparing base [(`8be1e8a`)](https://app.codecov.io/gh/duartegroup/autodE/commit/8be1e8a66b106dc30f8fb31703da473e81c7c952?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) to head [(`009d86e`)](https://app.codecov.io/gh/duartegroup/autodE/commit/009d86e0929683608c18f61cb552a2cb487736eb?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup).

| [Files](https://app.codecov.io/gh/duartegroup/autodE/pull/351?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) | Patch % | Lines |
|---|---|---|
| [autode/bracket/dhs.py](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&filepath=autode%2Fbracket%2Fdhs.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL2JyYWNrZXQvZGhzLnB5) | 20.00% | [12 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [autode/opt/optimisers/base.py](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&filepath=autode%2Fopt%2Foptimisers%2Fbase.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL29wdC9vcHRpbWlzZXJzL2Jhc2UucHk=) | 92.10% | [9 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [tests/test\_opt/test\_crfo.py](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&filepath=tests%2Ftest_opt%2Ftest_crfo.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-dGVzdHMvdGVzdF9vcHQvdGVzdF9jcmZvLnB5) | 46.15% | [7 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [tests/test\_opt/test\_opt.py](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&filepath=tests%2Ftest_opt%2Ftest_opt.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-dGVzdHMvdGVzdF9vcHQvdGVzdF9vcHQucHk=) | 64.28% | [5 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [autode/opt/coordinates/dic.py](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&filepath=autode%2Fopt%2Fcoordinates%2Fdic.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL29wdC9jb29yZGluYXRlcy9kaWMucHk=) | 84.61% | [2 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [tests/test\_opt/test\_prfo.py](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&filepath=tests%2Ftest_opt%2Ftest_prfo.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-dGVzdHMvdGVzdF9vcHQvdGVzdF9wcmZvLnB5) | 0.00% | [1 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [tests/test\_opt/test\_qa.py](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&filepath=tests%2Ftest_opt%2Ftest_qa.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-dGVzdHMvdGVzdF9vcHQvdGVzdF9xYS5weQ==) | 50.00% | [1 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/351?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |

> :exclamation:  There is a different number of reports uploaded between BASE (8be1e8a) and HEAD (009d86e). Click for more details.
> 
> <details><summary>HEAD has 1 upload less than BASE</summary>
>
>| Flag | BASE (8be1e8a) | HEAD (009d86e) |
>|------|------|------|
>|unittests|2|1|
></details>

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           v1.4.4     #351      +/-   ##
==========================================
- Coverage   97.46%   88.77%   -8.70%     
==========================================
  Files         204      204              
  Lines       23632    23738     +106     
==========================================
- Hits        23034    21074    -1960     
- Misses        598     2664    +2066     
```

| [Flag](https://app.codecov.io/gh/duartegroup/autodE/pull/351/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) | Coverage Δ | |
|---|---|---|
| [unittests](https://app.codecov.io/gh/duartegroup/autodE/pull/351/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) | `88.77% <82.21%> (-8.70%)` | :arrow_down: |

Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#carryforward-flags-in-the-pull-request-comment) to find out more.

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/duartegroup/autodE/pull/351?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup).

shoubhikraj: @t-young31 Could you have one final check before merge please?

Patch:
diff --git a/autode/bracket/dhs.py b/autode/bracket/dhs.py
index 0a9c39b..bb3df01 100644
--- a/autode/bracket/dhs.py
+++ b/autode/bracket/dhs.py
@@ -15,13 +15,14 @@ from autode.opt.coordinates import CartesianCoordinates
 from autode.opt.optimisers.utils import TruncatedTaylor
 from autode.opt.optimisers.hessian_update import BFGSSR1Update
 from autode.bracket.base import BaseBracketMethod
-from autode.opt.optimisers import RFOptimiser
+from autode.opt.optimisers import RFOptimiser, ConvergenceParams
 from autode.exceptions import OptimiserStepError
 from autode.log import logger
 
 if TYPE_CHECKING:
     from autode.species.species import Species
     from autode.wrappers.methods import Method
+    from autode.opt.optimisers.base import ConvergenceTolStr
 
 
 class DistanceConstrainedOptimiser(RFOptimiser):
@@ -60,10 +61,10 @@ class DistanceConstrainedOptimiser(RFOptimiser):
             pivot_point: Coordinates of the pivot point
             line_search: Whether to use linear search
             angle_thresh: An angle threshold above which linear search
-                          will be rejected (in Degrees)
+                        will be rejected (in Degrees)
             old_coords_read_hess: Old coordinate with hessian which will
-                                  be used to obtain initial hessian by
-                                  a Hessian update scheme
+                        be used to obtain the initial hessian by a
+                        Hessian update scheme
         """
         kwargs.pop("init_alpha", None)
         super().__init__(*args, init_alpha=init_trust, **kwargs)
@@ -103,15 +104,22 @@ class DistanceConstrainedOptimiser(RFOptimiser):
     @property
     def converged(self) -> bool:
         """Has the optimisation converged"""
-        # The tangential gradient should be close to zero
-        return (
-            self.rms_tangent_grad < self.gtol and self._abs_delta_e < self.etol
-        )
+        assert self._coords is not None
+
+        # Check only the tangential component of gradient
+        g_tau = self.tangent_grad
+        rms_g_tau = np.sqrt(np.mean(np.square(g_tau)))
+        max_g_tau = np.max(np.abs(g_tau))
+
+        curr_params = self._history.conv_params()
+        curr_params.rms_g = GradientRMS(rms_g_tau, "Ha/ang")
+        curr_params.max_g = GradientRMS(max_g_tau, "Ha/ang")
+        return self.conv_tol.meets_criteria(curr_params)
 
     @property
-    def rms_tangent_grad(self) -> GradientRMS:
+    def tangent_grad(self) -> np.ndarray:
         """
-        Obtain the RMS of the gradient tangent to the distance
+        Obtain the component of atomic gradients tangent to the distance
         vector between current coords and pivot point
         """
         assert self._coords is not None and self._coords.g is not None
@@ -120,8 +128,7 @@ class DistanceConstrainedOptimiser(RFOptimiser):
         # unit vector in the direction of distance vector
         d_hat = self.dist_vec / np.linalg.norm(self.dist_vec)
         tangent_grad = grad - (grad.dot(d_hat)) * d_hat
-        rms_grad = np.sqrt(np.mean(np.square(tangent_grad)))
-        return GradientRMS(rms_grad)
+        return tangent_grad
 
     @property
     def dist_vec(self) -> np.ndarray:
@@ -491,6 +498,7 @@ class DHS(BaseBracketMethod):
         large_step: Union[Distance, float] = Distance(0.2, "ang"),
         small_step: Union[Distance, float] = Distance(0.05, "ang"),
         switch_thresh: Union[Distance, float] = Distance(1.5, "ang"),
+        conv_tol: Union["ConvergenceParams", "ConvergenceTolStr"] = "loose",
         **kwargs,
     ):
         """
@@ -513,6 +521,9 @@ class DHS(BaseBracketMethod):
             switch_thresh: When distance between the two images is less than
                         this cutoff, smaller DHS extrapolation steps are taken
 
+            conv_tol: Convergence tolerance for the distance-constrained
+                      optimiser
+
         Keyword Args:
 
             maxiter: Maximum number of en/grad evaluations
@@ -521,9 +532,6 @@ class DHS(BaseBracketMethod):
                       stop, values less than 0.5 Angstrom are not
                       recommended.
 
-            gtol: Gradient tolerance for the optimiser micro-iterations
-                  in DHS (Hartree/angstrom)
-
             cineb_at_conv: Whether to run CI-NEB calculation from the end
                            points after the DHS is converged
         """
@@ -542,6 +550,7 @@ class DHS(BaseBracketMethod):
         self._small_step = Distance(abs(small_step), "ang")
         self._sw_thresh = Distance(abs(switch_thresh), "ang")
         assert self._small_step < self._large_step
+        self._conv_tol = conv_tol
 
         self._step_size: Optional[Distance] = None
         if self._large_step > self.imgpair.dist:
@@ -597,8 +606,7 @@ class DHS(BaseBracketMethod):
 
         opt = DistanceConstrainedOptimiser(
             maxiter=curr_maxiter,
-            gtol=self._gtol,
-            etol=1.0e-3,  # seems like a reasonable etol
+            conv_tol=self._conv_tol,
             init_trust=opt_trust,
             pivot_point=pivot,
             old_coords_read_hess=old_coords,
@@ -612,9 +620,10 @@ class DHS(BaseBracketMethod):
         if not opt.converged:
             return None
 
+        rms_g_tau = np.sqrt(np.mean(np.square(opt.tangent_grad)))
         logger.info(
             "Successful optimization after DHS step, final RMS of "
-            f"tangential gradient = {opt.rms_tangent_grad:.6f} "
+            f"tangential gradient = {rms_g_tau:.6f} "
             f"Ha/angstrom"
         )
 
diff --git a/autode/calculations/executors.py b/autode/calculations/executors.py
index 2e7ebe9..68ac1b3 100644
--- a/autode/calculations/executors.py
+++ b/autode/calculations/executors.py
@@ -346,10 +346,7 @@ class CalculationExecutorO(_IndirectCalculationExecutor):
     def __init__(self, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
-        self.etol = PotentialEnergy(3e-5, units="Ha")
-        self.gtol = GradientRMS(
-            1e-3, units="Ha Å^-1"
-        )  # TODO: A better number here
+        self.conv_tol = "normal"
         self._fix_unique()
 
     def run(self) -> None:
@@ -367,8 +364,7 @@ class CalculationExecutorO(_IndirectCalculationExecutor):
         self.optimiser: "NDOptimiser" = type_(
             init_alpha=self._step_size,
             maxiter=self._max_opt_cycles,
-            etol=self.etol,
-            gtol=self.gtol,
+            conv_tol=self.conv_tol,
         )
         method = self.method.copy()
         method.keywords.grad = kws.GradientKeywords(self.input.keywords)
diff --git a/autode/opt/coordinates/base.py b/autode/opt/coordinates/base.py
index d6cb73e..8b26086 100644
--- a/autode/opt/coordinates/base.py
+++ b/autode/opt/coordinates/base.py
@@ -343,6 +343,16 @@ class OptCoordinates(ValueArray, ABC):
     def iadd(self, value: np.ndarray) -> "OptCoordinates":
         """Inplace addition of some coordinates"""
 
+    @property
+    @abstractmethod
+    def n_constraints(self) -> int:
+        """Number of constraints in these coordinates"""
+
+    @property
+    @abstractmethod
+    def n_satisfied_constraints(self) -> int:
+        """Number of constraints that are satisfied in these coordinates"""
+
     @property
     @abstractmethod
     def active_indexes(self) -> List[int]:
@@ -358,6 +368,13 @@ class OptCoordinates(ValueArray, ABC):
     def inactive_indexes(self) -> List[int]:
         """A list of indexes which are non-active in this coordinate set"""
 
+    @property
+    @abstractmethod
+    def cart_proj_g(self) -> Optional[np.ndarray]:
+        """
+        The Cartesian gradient with any constraints projected out
+        """
+
     def __eq__(self, other):
         """Coordinates can never be identical..."""
         return False
diff --git a/autode/opt/coordinates/cartesian.py b/autode/opt/coordinates/cartesian.py
index 08cc192..7b9c00f 100644
--- a/autode/opt/coordinates/cartesian.py
+++ b/autode/opt/coordinates/cartesian.py
@@ -9,7 +9,6 @@ from autode.opt.coordinates.dic import DIC
 if TYPE_CHECKING:
     from autode.values import Gradient
     from autode.hessians import Hessian
-    from autode.opt.coordinates.primitives import ConstrainedPrimitive
 
 
 class CartesianCoordinates(OptCoordinates):
@@ -61,6 +60,14 @@ class CartesianCoordinates(OptCoordinates):
         assert self.h_or_h_inv_has_correct_shape(arr)
         self._h = None if arr is None else np.array(arr)
 
+    @property
+    def n_constraints(self) -> int:
+        return 0
+
+    @property
+    def n_satisfied_constraints(self) -> int:
+        return 0
+
     @property
     def active_indexes(self) -> List[int]:
         return list(range(len(self)))
@@ -106,6 +113,10 @@ class CartesianCoordinates(OptCoordinates):
                 f"Cannot convert Cartesian coordinates to {value}"
             )
 
+    @property
+    def cart_proj_g(self) -> Optional[np.ndarray]:
+        return self.g
+
     @property
     def expected_number_of_dof(self) -> int:
         """Expected number of degrees of freedom for the system"""
diff --git a/autode/opt/coordinates/dic.py b/autode/opt/coordinates/dic.py
index 44d4418..e0a39d9 100644
--- a/autode/opt/coordinates/dic.py
+++ b/autode/opt/coordinates/dic.py
@@ -139,6 +139,10 @@ class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
         logger.info(f"Transformed in      ...{time() - start_time:.4f} s")
         return dic
 
+    @property
+    def cart_proj_g(self) -> Optional[np.ndarray]:
+        return self.to("cart").g
+
     def _update_g_from_cart_g(self, arr: Optional["Gradient"]) -> None:
         """
         Updates the gradient from a calculated Cartesian gradient
@@ -422,6 +426,20 @@ class DICWithConstraints(DIC):
 
         return super().iadd(delta_s)
 
+    @property
+    def cart_proj_g(self) -> Optional[np.ndarray]:
+        """Obtain Cartesian gradient with constraints projected out"""
+        if self.g is None:
+            return None
+        # constrained gradient with inactive terms set to zero
+        g_s = self.g
+        g_s[self.inactive_indexes] = 0.0
+        g_s = g_s[: len(self)]
+        # back to Cartesian
+        g_x = np.matmul(self.B.T, g_s)
+        assert len(g_x) == len(self.to("cart"))
+        return g_x
+
     @property
     def g(self):
         """
diff --git a/autode/opt/optimisers/__init__.py b/autode/opt/optimisers/__init__.py
index 3b355dc..23d5f25 100644
--- a/autode/opt/optimisers/__init__.py
+++ b/autode/opt/optimisers/__init__.py
@@ -1,7 +1,8 @@
-from autode.opt.optimisers.base import NDOptimiser
+from autode.opt.optimisers.base import NDOptimiser, ConvergenceParams
 from autode.opt.optimisers.rfo import RFOptimiser
 from autode.opt.optimisers.prfo import PRFOptimiser
 from autode.opt.optimisers.crfo import CRFOptimiser
+from autode.opt.optimisers.qa import QAOptimiser
 from autode.opt.optimisers.steepest_descent import (
     CartesianSDOptimiser,
     DIC_SD_Optimiser,
@@ -9,9 +10,11 @@ from autode.opt.optimisers.steepest_descent import (
 
 __all__ = [
     "NDOptimiser",
+    "ConvergenceParams",
     "RFOptimiser",
     "PRFOptimiser",
     "CRFOptimiser",
+    "QAOptimiser",
     "CartesianSDOptimiser",
     "DIC_SD_Optimiser",
 ]
diff --git a/autode/opt/optimisers/base.py b/autode/opt/optimisers/base.py
index 0e3468c..d01a306 100644
--- a/autode/opt/optimisers/base.py
+++ b/autode/opt/optimisers/base.py
@@ -1,25 +1,17 @@
 import os
 import pickle
-
+from dataclasses import dataclass, fields
 import numpy as np
 
 from abc import ABC, abstractmethod
 from zipfile import ZipFile, is_zipfile
 from collections import deque
-from typing import (
-    Type,
-    List,
-    Union,
-    Optional,
-    Callable,
-    Any,
-    TYPE_CHECKING,
-    Iterator,
-)
+from typing import Type, List, Union, Optional, Callable, Any
+from typing import TYPE_CHECKING, Iterator, Literal
 
 from autode.log import logger
 from autode.config import Config
-from autode.values import GradientRMS, PotentialEnergy, method_string
+from autode.values import GradientRMS, PotentialEnergy, method_string, Distance
 from autode.opt.coordinates.base import OptCoordinates
 from autode.opt.optimisers.hessian_update import NullUpdate
 from autode.exceptions import CalculationException
@@ -170,7 +162,6 @@ class Optimiser(BaseOptimiser, ABC):
             f"with {self._n_cores} cores using {self._maxiter} max "
             f"iterations"
         )
-        logger.info("Iteration\t|∆E| / \\kcal mol-1 \t||∇E|| / Ha Å-1")
 
         while not self.converged:
             self._callback(self._coords)
@@ -466,14 +457,220 @@ class NullOptimiser(BaseOptimiser):
         raise RuntimeError("A NullOptimiser has no coordinates")
 
 
+ConvergenceTolStr = Literal["loose", "normal", "tight", "verytight"]
+
+
+@dataclass
+class ConvergenceParams:
+    """
+    Various convergence parameters for optimisers and some common
+    preset convergence tolerances
+
+    Args:
+        abs_d_e: Absolute change in energy, |E_i - E_i-1|
+        rms_g: RMS of the gradient, RMS(∇E)
+        max_g: Maximum component of gradient, max(∇E)
+        rms_s: RMS of the last step, RMS(x_i - x_i-1)
+        max_s: Maximum component of last step, max(x_i - x_i-1)
+        strict: Whether all criteria must be converged strictly.
+                If False, convergence is signalled when some criteria
+                are overachieved and others are close to convergence
+
+    """
+
+    abs_d_e: Optional[PotentialEnergy] = None
+    rms_g: Optional[GradientRMS] = None
+    max_g: Optional[GradientRMS] = None
+    rms_s: Optional[Distance] = None
+    max_s: Optional[Distance] = None
+    strict: bool = False
+
+    @property
+    def _num_attrs(self) -> List[str]:
+        """Numerical attributes of this dataclass, in order"""
+        return ["abs_d_e", "rms_g", "max_g", "rms_s", "max_s"]
+
+    def __post_init__(self):
+        """Type checking and sanity checks on parameters"""
+
+        # convert units for easier comparison
+        self._to_base_units()
+        self.strict = bool(self.strict)
+        # RMS(g) is the most basic criteria that is always needed
+        if self.rms_g is None:
+            raise ValueError(
+                "At least the RMS gradient criteria has to be defined!"
+            )
+
+        for attr in self._num_attrs:
+            if getattr(self, attr) is None:
+                continue
+            if not getattr(self, attr) > 0:
+                raise ValueError(
+                    f"Value of {attr} should be positive"
+                    f" but set to {getattr(self, attr)}!"
+                )
+
+    def _to_base_units(self) -> None:
+        """
+        Convert all set criteria to the default units in terms of
+        Hartree and Angstrom, and also ensure everything has units
+        """
+        if self.abs_d_e is not None:
+            self.abs_d_e = PotentialEnergy(self.abs_d_e).to("Ha")
+        if self.rms_g is not None:
+            self.rms_g = GradientRMS(self.rms_g).to("Ha/ang")
+        if self.max_g is not None:
+            self.max_g = GradientRMS(self.max_g).to("Ha/ang")
+        if self.rms_s is not None:
+            self.rms_s = Distance(self.rms_s).to("ang")
+        if self.max_s is not None:
+            self.max_s = Distance(self.max_s).to("ang")
+        return None
+
+    @classmethod
+    def from_preset(cls, preset_name: str) -> "ConvergenceParams":
+        """
+        Obtains preset values of convergence criteria - given as
+        "loose", "normal", "tight" and "verytight".
+
+        Args:
+            preset_name: Must be one of the strings "loose", "normal"
+                    "tight" or "verytight"
+
+        Returns:
+            (ConvergenceCriteria): Optimiser convergence criteria, with
+                    preset values
+        """
+        # NOTE: Taken from ORCA
+        preset_dicts = {
+            "loose": {
+                "abs_d_e": PotentialEnergy(3e-5, "Ha"),
+                "rms_g": GradientRMS(5e-4, "Ha/bohr").to("Ha/ang"),
+                "max_g": GradientRMS(2e-3, "Ha/bohr").to("Ha/ang"),
+                "rms_s": Distance(7e-3, "bohr").to("ang"),
+                "max_s": Distance(1e-2, "bohr").to("ang"),
+            },
+            "normal": {
+                "abs_d_e": PotentialEnergy(5e-6, "Ha"),
+                "rms_g": GradientRMS(1e-4, "Ha/bohr").to("Ha/ang"),
+                "max_g": GradientRMS(3e-4, "Ha/bohr").to("Ha/ang"),
+                "rms_s": Distance(2e-3, "bohr").to("ang"),
+                "max_s": Distance(4e-3, "bohr").to("ang"),
+            },
+            "tight": {
+                "abs_d_e": PotentialEnergy(1e-6, "Ha"),
+                "rms_g": GradientRMS(3e-5, "Ha/bohr").to("Ha/ang"),
+                "max_g": GradientRMS(1e-4, "Ha/bohr").to("Ha/ang"),
+                "rms_s": Distance(6e-4, "bohr").to("ang"),
+                "max_s": Distance(1e-3, "bohr").to("ang"),
+            },
+            "verytight": {
+                "abs_d_e": PotentialEnergy(2e-7, "Ha"),
+                "rms_g": GradientRMS(8e-6, "Ha/bohr").to("Ha/ang"),
+                "max_g": GradientRMS(3e-5, "Ha/bohr").to("Ha/ang"),
+                "rms_s": Distance(1e-4, "bohr").to("ang"),
+                "max_s": Distance(2e-4, "bohr").to("ang"),
+            },
+        }
+
+        allowed_strs = list(preset_dicts.keys())
+        preset_name = preset_name.strip().lower()
+        if preset_name not in allowed_strs:
+            raise ValueError(
+                f"Unknown preset convergence: {preset_name}, please select"
+                f" from {allowed_strs}"
+            )
+
+        return cls(**preset_dicts[preset_name])
+
+    def __mul__(self, factors: List[float]):
+        """Multiply a set of criteria with ordered list of numerical factors"""
+        assert len(factors) == len(self._num_attrs)
+        kwargs = {}
+        for idx, attr in enumerate(self._num_attrs):
+            c = getattr(self, attr)
+            if c is not None:
+                kwargs[attr] = getattr(self, attr) * factors[idx]
+            else:
+                kwargs[attr] = None
+        return ConvergenceParams(**kwargs, strict=self.strict)
+
+    def are_satisfied(self, other: "ConvergenceParams") -> List[bool]:
+        """
+        Return an elementwise comparison between the current criteria
+        and another set of parameters (comparing only numerical attributes)
+
+        Args:
+            other: Another set of parameters
+
+        Returns:
+            (list[bool]): List containing True or False
+        """
+        are_satisfied = []
+
+        # unset criteria are always satisfied
+        for attr in self._num_attrs:
+            c = getattr(self, attr)
+            v = getattr(other, attr)
+            if c is None:
+                are_satisfied.append(True)
+            else:
+                are_satisfied.append(float(v) <= float(c))
+        return are_satisfied
+
+    def meets_criteria(self, other: "ConvergenceParams") -> bool:
+        """
+        Does a set of parameters satisfy the current convergence criteria?
+        Will signal convergence if gradient or energy change are overachieved
+        or all other criteria except energy is satisfied
+
+        Args:
+            other (ConvergenceParams): Another set of parameters to be
+                            checked against the current set
+
+        Returns:
+            (bool):
+        """
+        # everything satisfied - simplest case
+        if all(self.are_satisfied(other)):
+            return True
+
+        # strict = everything must be converged
+        elif self.strict:
+            return False
+
+        if all((self * [0.5, 0.5, 0.8, 3, 3]).are_satisfied(other)):
+            logger.warning(
+                "Overachieved gradient and energy convergence, reasonable "
+                "convergence on step size."
+            )
+            return True
+
+        if all((self * [1.5, 0.1, 0.2, 2, 2]).are_satisfied(other)):
+            logger.warning(
+                "Gradient is one order of magnitude below convergence, "
+                "other parameter(s) are almost converged."
+            )
+            return True
+
+        if all((self * [3, 0.7, 0.7, 1, 1]).are_satisfied(other)):
+            logger.warning(
+                "Everything except energy has been converged. Reasonable"
+                " convergence on energy"
+            )
+            return True
+
+        return False
+
+
 class NDOptimiser(Optimiser, ABC):
     """Abstract base class for an optimiser in N-dimensions"""
 
     def __init__(
         self,
         maxiter: int,
-        gtol: GradientRMS,
-        etol: PotentialEnergy,
+        conv_tol: Union[ConvergenceParams, ConvergenceTolStr],
         coords: Optional[OptCoordinates] = None,
         **kwargs,
     ):
@@ -486,70 +683,56 @@ class NDOptimiser(Optimiser, ABC):
         Arguments:
             maxiter (int): Maximum number of iterations to perform
 
-            gtol (autode.values.GradientRMS): Tolerance on RMS(|∇E|)
-
-            etol (autode.values.PotentialEnergy): Tolerance on |E_i+1 - E_i|
+            conv_tol (ConvergenceParams|ConvergenceTolStr): Convergence tolerances,
+                    indicating thresholds for absolute energy change (|E_i+1 - E_i|),
+                    RMS and max. gradients (∇E) and RMS and max. step size (Δx)
+                    Either supplied as a dictionary or a ConvergenceParams object
 
         See Also:
 
             :py:meth:`Optimiser <Optimiser.__init__>`
+            :py:meth:`ConvergenceParams <ConvergenceParams.__init__>`
         """
         super().__init__(maxiter=maxiter, coords=coords, **kwargs)
 
-        self.etol = etol
-        self.gtol = gtol
-
+        if isinstance(conv_tol, str):
+            conv_tol = ConvergenceParams.from_preset(conv_tol)
+        self.conv_tol = conv_tol
         self._hessian_update_types: List[Type[HessianUpdater]] = [NullUpdate]
 
     @property
-    def gtol(self) -> GradientRMS:
+    def conv_tol(self) -> "ConvergenceParams":
         """
-        Gradient tolerance on |∇E| i.e. the root mean square of each component
+        All convergence parameters for this optimiser. If numerical
+        values are unset, they appear as infinity.
 
-        -----------------------------------------------------------------------
         Returns:
-            (autode.values.GradientRMS):
+            (ConvergenceParams):
         """
-        return self._gtol
+        return self._conv_tol
 
-    @gtol.setter
-    def gtol(self, value: Union[int, float, GradientRMS]):
-        """Set the gradient tolerance"""
-
-        if float(value) <= 0:
-            raise ValueError(
-                "Tolerance on the gradient (||∇E||) must be "
-                f"positive. Had: gtol={value}"
-            )
-
-        self._gtol = GradientRMS(value)
-
-    @property
-    def etol(self) -> PotentialEnergy:
+    @conv_tol.setter
+    def conv_tol(self, value: Union["ConvergenceParams", ConvergenceTolStr]):
         """
-        Energy tolerance between two consecutive steps of the optimisation
+        Set the convergence parameters for this optimiser.
 
-        -----------------------------------------------------------------------
-        Returns:
-            (autode.values.PotentialEnergy): Energy tolerance
+        Args:
+            value (ConvergenceParams|str):
         """
-        return self._etol
-
-    @etol.setter
-    def etol(self, value: Union[int, float, PotentialEnergy]):
-        """Set the energy tolerance"""
-        if float(value) <= 0:
+        if isinstance(value, str):
+            self._conv_tol = ConvergenceParams.from_preset(value)
+        elif isinstance(value, ConvergenceParams):
+            self._conv_tol = value
+        else:
             raise ValueError(
-                "Tolerance on the energy change is absolute so "
-                f"must be positive. Had etol = {value}"
+                "Convergence tolerance should be of type ConvergenceParams"
+                f" or a preset string, but assigned {type(value)}"
             )
 
-        self._etol = PotentialEnergy(value)
-
     @property
     def optimiser_params(self):
         """Optimiser params to save"""
-        return {"maxiter": self._maxiter, "gtol": self.gtol, "etol": self.etol}
+        return {"maxiter": self._maxiter, "conv_tol": self.conv_tol}
 
     @classmethod
     def optimise(
@@ -559,8 +742,7 @@ class NDOptimiser(Optimiser, ABC):
         n_cores: Optional[int] = None,
         coords: Optional[OptCoordinates] = None,
         maxiter: int = 100,
-        gtol: Any = GradientRMS(1e-3, units="Ha Å-1"),
-        etol: Any = PotentialEnergy(1e-4, units="Ha"),
+        conv_tol: Union[ConvergenceParams, ConvergenceTolStr] = "normal",
         **kwargs,
     ) -> None:
         """
@@ -574,12 +756,9 @@ class NDOptimiser(Optimiser, ABC):
 
             maxiter (int): Maximum number of iteration to perform
 
-            gtol (float | autode.values.GradientNorm): Tolerance on RMS(|∇E|)
-                 i.e. the root mean square of the gradient components. If
-                 a float then assume units of Ha Å^-1
-
-            etol (float | autode.values.PotentialEnergy): Tolerance on |∆E|
-                 between two consecutive iterations of the optimiser
+            conv_tol (ConvergenceParams|ConvergenceTolStr): Convergence parameters
+                        for the absolute energy change, RMS and max gradient,
+                        and RMS and max step sizes.
 
             coords (OptCoordinates | None): Coordinates to optimise in
 
@@ -594,7 +773,7 @@ class NDOptimiser(Optimiser, ABC):
         """
 
         optimiser = cls(
-            maxiter=maxiter, gtol=gtol, etol=etol, coords=coords, **kwargs
+            maxiter=maxiter, conv_tol=conv_tol, coords=coords, **kwargs
         )
         optimiser.run(species, method, n_cores=n_cores)
 
@@ -607,8 +786,8 @@ class NDOptimiser(Optimiser, ABC):
     @property
     def converged(self) -> bool:
         """
-        Is this optimisation converged? Must be converged based on both energy
-        and gradient tolerance.
+        Is this optimisation converged? Must be converged based on energy, gradient
+        and step size criteria.
 
         -----------------------------------------------------------------------
         Returns:
@@ -617,15 +796,13 @@ class NDOptimiser(Optimiser, ABC):
         if self._species is not None and self._species.n_atoms == 1:
             return True  # Optimisation 0 DOF is always converged
 
-        if self._abs_delta_e < self.etol / 10:
-            logger.warning(
-                f"Energy change is overachieved. "
-                f'{self.etol.to("kcal") / 10:.3E} kcal mol-1. '
-                f"Signaling convergence"
-            )
-            return True
+        assert self._coords is not None, "Must have coordinates!"
+        curr_params = self._history.conv_params()
 
-        return self._abs_delta_e < self.etol and self._g_norm < self.gtol
+        constrs_met = (
+            self._coords.n_constraints == self._coords.n_satisfied_constraints
+        )
+        return constrs_met and self.conv_tol.meets_criteria(curr_params)
 
     def clean_up(self) -> None:
         """
@@ -643,69 +820,29 @@ class NDOptimiser(Optimiser, ABC):
         optimiser._history = hist
         return optimiser
 
-    @property
-    def _abs_delta_e(self) -> PotentialEnergy:
-        """
-        Calculate the absolute energy difference
-
-        .. math::
-            |∆E| = |E_i - E_{i-1}|   for a step i
-
-        -----------------------------------------------------------------------
-        Returns:
-            (autode.values.PotentialEnergy): Energy difference. Infinity if
-                                  an energy difference cannot be calculated
-        """
-        assert (
-            self._coords is not None
-        ), "Must have coordinates to calculate ∆E"
-
-        if len(self._history) < 2:
-            logger.info("First iteration - returning |∆E| = ∞")
-            return PotentialEnergy(np.inf)
-
-        e1, e2 = self._coords.e, self._history.penultimate.e
-
-        if e1 is None or e2 is None:
-            logger.error(
-                "Cannot determine absolute energy difference. Using |∆E| = ∞"
-            )
-            return PotentialEnergy(np.inf)
-
-        return PotentialEnergy(abs(e1 - e2))  # type: ignore
-
-    @property
-    def _g_norm(self) -> GradientRMS:
-        """
-        Calculate RMS(∇E) based on the current Cartesian gradient.
-
-        -----------------------------------------------------------------------
-        Returns:
-            (autode.values.GradientRMS): Gradient norm. Infinity if the
-                                          gradient is not defined
-        """
-        if self._coords is None:
-            logger.warning("Had no coordinates - cannot determine ||∇E||")
-            return GradientRMS(np.inf)
-
-        if self._coords.g is None:
-            return GradientRMS(np.inf)
-
-        return GradientRMS(np.sqrt(np.mean(np.square(self._coords.g))))
-
     def _log_convergence(self) -> None:
-        """Log the convergence of the energy"""
-        assert (
-            self._coords is not None
-        ), "Must have coordinates to log convergence"
-        log_string = f"{self.iteration}\t"
-
-        if len(self._history) > 1:
-            assert self._coords.e and self._history.penultimate.e, "Need ∆E"
-            de: PotentialEnergy = self._coords.e - self._history.penultimate.e
-            log_string += f'{de.to("kcal mol-1"):.3f}\t{self._g_norm:.5f}'
-
-        logger.info(log_string)
+        """Log the convergence of the all convergence parameters"""
+        assert self._coords is not None, "Must have coordinates!"
+
+        curr_params = self._history.conv_params()
+        assert curr_params.abs_d_e is not None
+
+        conv_msgs = [
+            "(YES)" if param_converged else "(NO)"
+            for param_converged in self.conv_tol.are_satisfied(curr_params)
+        ]
+        log_string1 = (
+            f"iter# {self.iteration}   |dE|="
+            f"{curr_params.abs_d_e.to('kcalmol'):.5f} kcal/mol {conv_msgs[0]}"
+            f" RMS(g)={curr_params.rms_g:.5f} Ha/Å {conv_msgs[1]} "
+        )
+        log_string2 = (
+            f"max(g)={curr_params.max_g:.5f} Ha/Å {conv_msgs[2]}  "
+            f"RMS(dx)={curr_params.rms_s:.5f} Å {conv_msgs[3]}  "
+            f"max(dx)={curr_params.max_s:.5f} Å  {conv_msgs[4]}"
+        )
+        logger.info(log_string1)
+        logger.info(log_string2)
         return None
 
     def plot_optimisation(
@@ -1070,6 +1207,43 @@ class OptimiserHistory:
         for i in reversed(range(len(self))):
             yield self[i]
 
+    def conv_params(self, idx: int = -1) -> ConvergenceParams:
+        """
+        Calculate the convergence parameters for the coordinates at
+        specified index (default -1 i.e. the last set of coordinates)
+
+        Args:
+            idx (int): Index of the set of coordinates for which to
+                calculate the parameter
+
+        Returns:
+            (ConvergenceParams):
+        """
+        # NOTE: Internal coordinates have inconsistent units, so we
+        # calculate step sizes and gradients in Cartesian coordinates
+        coords_l = self[idx]
+        assert coords_l is not None
+        g_x = coords_l.cart_proj_g
+        if g_x is not None:
+            rms_g = np.sqrt(np.mean(np.square(g_x)))
+            max_g = np.max(np.abs(g_x))
+        else:
+            rms_g = max_g = np.inf
+
+        if len(self) > 1:
+            coords_k = self[idx - 1]
+            assert coords_k is not None
+            assert coords_l.e is not None and coords_k.e is not None
+            abs_d_e = PotentialEnergy(abs(coords_l.e - coords_k.e))
+            delta_x = coords_l.to("cart") - coords_k.to("cart")
+            rms_s = np.sqrt(np.mean(np.square(delta_x)))
+            max_s = np.max(np.abs(delta_x))
+        else:
+            abs_d_e = rms_s = max_s = np.inf
+        return ConvergenceParams(
+            abs_d_e=abs_d_e, rms_g=rms_g, max_g=max_g, rms_s=rms_s, max_s=max_s
+        )
+
 
 class ExternalOptimiser(BaseOptimiser, ABC):
     @property
diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index de9df67..0d30a76 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -235,16 +235,6 @@ class CRFOptimiser(RFOptimiser):
             f" Current trust radius = {self.alpha:.3f}"
         )
 
-    @property
-    def _g_norm(self) -> GradientRMS:
-        """Calculate the norm of the gradient in the active subspace"""
-
-        if self._coords is None or self._coords.g is None:
-            return super()._g_norm
-
-        gradient = self._coords.g[self._coords.active_indexes]
-        return GradientRMS(np.sqrt(np.mean(np.square(gradient))))
-
     def _initialise_run(self) -> None:
         """Initialise the optimisation"""
         logger.info("Initialising optimisation")
diff --git a/autode/opt/optimisers/prfo.py b/autode/opt/optimisers/prfo.py
index 3fb78a8..6a5a741 100644
--- a/autode/opt/optimisers/prfo.py
+++ b/autode/opt/optimisers/prfo.py
@@ -24,6 +24,7 @@ class PRFOptimiser(CRFOptimiser):
         all others to locate a transition state (TS)
 
         -----------------------------------------------------------------------
+
         Arguments:
             init_alpha: Maximum step size (default Å if unit not given)
 
diff --git a/autode/opt/optimisers/steepest_descent.py b/autode/opt/optimisers/steepest_descent.py
index 8edd3c6..e67980e 100644
--- a/autode/opt/optimisers/steepest_descent.py
+++ b/autode/opt/optimisers/steepest_descent.py
@@ -4,7 +4,7 @@ from autode.opt.optimisers.base import NDOptimiser
 
 
 class SteepestDescent(NDOptimiser, ABC):
-    def __init__(self, maxiter, gtol, etol, step_size=0.2, **kwargs):
+    def __init__(self, maxiter, conv_tol, step_size=0.2, **kwargs):
         """
         Steepest decent optimiser
 
@@ -16,7 +16,7 @@ class SteepestDescent(NDOptimiser, ABC):
 
             :py:meth:`NDOptimiser <autode.opt.optimisers.base.NDOptimiser.__init__>`
         """
-        super().__init__(maxiter=maxiter, gtol=gtol, etol=etol, **kwargs)
+        super().__init__(maxiter=maxiter, conv_tol=conv_tol, **kwargs)
 
         self.alpha = step_size
 
diff --git a/doc/changelog.rst b/doc/changelog.rst
index b66d774..fffecaf 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -16,6 +16,9 @@ Bug Fixes
 *********
 - DIC to Cartesian transform will now always use :code:`PIC.close_to()` to ensure steps along dihedral have the smallest change, even after back-transform is complete
 
+Usability improvements/Changes
+******************************
+- Optimiser convergence criteria have been improved to consider energy change, RMS and max. gradient and step sizes.
 
 1.4.3
 ------


Test patch:
diff --git a/tests/test_bracket/test_dhs.py b/tests/test_bracket/test_dhs.py
index 1c23940..21091e8 100644
--- a/tests/test_bracket/test_dhs.py
+++ b/tests/test_bracket/test_dhs.py
@@ -41,8 +41,7 @@ def test_distance_constrained_optimiser():
         pivot_point=rct_coords,
         maxiter=1,  # just one step
         init_trust=0.2,
-        gtol=1e-3,
-        etol=1e-4,
+        conv_tol="loose",
     )
     opt.run(product, method=XTB())
     assert not opt.converged
@@ -58,10 +57,7 @@ def test_distance_constrained_optimiser():
     assert step_size <= 0.2 + fp_err * 0.2  # for floating point error
 
     opt = DistanceConstrainedOptimiser(
-        pivot_point=rct_coords,
-        maxiter=50,
-        gtol=1e-3,
-        etol=1e-4,
+        pivot_point=rct_coords, maxiter=50, conv_tol="loose"
     )
     opt.run(product, method=XTB())
     assert opt.converged
@@ -77,13 +73,9 @@ def test_dist_constr_optimiser_sd_fallback():
     coords1.update_h_from_cart_h(np.loadtxt("conopt_last_h.txt"))
     pivot = CartesianCoordinates(np.loadtxt("conopt_pivot.txt"))
 
-    # lagrangian step may fail at certain point
+    # lagrangian step may fail at certain points
     opt = DistanceConstrainedOptimiser(
-        pivot_point=pivot,
-        maxiter=2,
-        init_trust=0.2,
-        gtol=1e-3,
-        etol=1e-4,
+        pivot_point=pivot, maxiter=2, init_trust=0.2, conv_tol="loose"
     )
     opt._target_dist = 2.6869833732268
     opt._history.open("test_trj")
@@ -105,10 +97,7 @@ def test_dist_constr_optimiser_energy_rising():
     coords2.e = PotentialEnergy(0.02, "Ha")
     assert (coords2.e - coords1.e) > PotentialEnergy(5, "kcalmol")
     opt = DistanceConstrainedOptimiser(
-        pivot_point=coords1,
-        maxiter=2,
-        gtol=1e-3,
-        etol=1e-4,
+        pivot_point=coords1, maxiter=2, conv_tol="loose"
     )
     opt._history.open("test_trj")
     opt._history.add(coords1)
@@ -257,9 +246,8 @@ def test_dhs_diels_alder():
         initial_species=reactant,
         final_species=product,
         small_step=0.2,
-        maxiter=200,
+        maxiter=100,
         dist_tol=set_dist_tol,
-        gtol=1.0e-3,
     )
 
     dhs.calculate(method=XTB(), n_cores=Config.n_cores)
@@ -318,11 +306,11 @@ def test_dhs_jumping_over_barrier(caplog):
     dhs = DHS(
         initial_species=reactant,
         final_species=product,
-        maxiter=200,
+        maxiter=50,
         large_step=0.5,
         small_step=0.5,
         dist_tol=0.3,  # smaller dist_tol also to make one side jump
-        gtol=5.0e-4,
+        conv_tol="loose",
         barrier_check=True,
         cineb_at_conv=True,
     )
@@ -346,11 +334,11 @@ def test_dhs_stops_if_microiter_exceeded(caplog):
     dhs = DHS(
         initial_species=reactant,
         final_species=product,
-        maxiter=10,
+        maxiter=5,
         large_step=0.2,
         small_step=0.1,
         dist_tol=1.0,
-        gtol=5.0e-4,
+        conv_tol="loose",
         barrier_check=True,
     )
     with caplog.at_level("WARNING"):
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index a761ea9..9e4ee3a 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -168,6 +168,25 @@ def test_hessian_inv():
     assert np.allclose(coords.h_inv, expected_h_inv, atol=1e-10)
 
 
+def test_h_update():
+    coords1 = CartesianCoordinates(np.array([1.0, 2.0]))
+    coords1.g = np.array([0.1, 0.2])
+    coords1.h = 2.0 * np.eye(2)
+
+    coords2 = CartesianCoordinates(np.array([1.1, 2.1]))
+    coords2.g = np.array([0.01, 0.02])
+    assert coords2.h is None
+
+    # must define a valid update type
+    with pytest.raises(RuntimeError):
+        coords2.update_h_from_old_h(coords1, hessian_update_types=[])
+
+    from autode.opt.optimisers.hessian_update import BFGSSR1Update
+
+    coords2.update_h_from_old_h(coords1, hessian_update_types=[BFGSSR1Update])
+    assert coords2.h is not None
+
+
 def test_cartesian_update_clear():
     arr = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0]])
 
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index 2d32f7f..d7a26ee 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -19,7 +19,7 @@ from ..testutils import requires_working_xtb_install
 
 
 def crfo_coords(molecule):
-    optimiser = CRFOptimiser(maxiter=1, gtol=1e-5, etol=1e-5)
+    optimiser = CRFOptimiser(maxiter=1, conv_tol="normal")
     optimiser._species = molecule
     optimiser._build_internal_coordinates()
 
@@ -49,7 +49,7 @@ def test_coordinate_setup():
     dist_consts = mol.constraints.distance
     assert (0, 1) in dist_consts and (1, 0) in dist_consts
 
-    opt = CRFOptimiser(maxiter=1, gtol=1e-5, etol=1e-5)
+    opt = CRFOptimiser(maxiter=1, conv_tol="normal")
 
     with pytest.raises(RuntimeError):
         # Cannot set coordinates without a species
@@ -111,7 +111,7 @@ def test_simple_hessian_update():
 
 
 def test_primitive_projection_discard():
-    optimiser = CRFOptimiser(etol=1, gtol=1, maxiter=1)
+    optimiser = CRFOptimiser(maxiter=1, conv_tol="loose")
     optimiser._species = water_molecule()
 
     # Current distance that will be constrained
@@ -132,15 +132,10 @@ def test_primitive_projection_discard():
     assert np.isclose(r(s.to("cartesian")), r_initial, atol=1e-10)
 
 
-def test_init_g_norm_is_none():
-    optimiser = CRFOptimiser(etol=1, gtol=1, maxiter=1)
-    assert optimiser._g_norm > 0
-
-
 def test_sanitised_zero_length_step():
     """Should be able to update with a null step"""
 
-    optimiser = CRFOptimiser(etol=1, gtol=1, maxiter=1)
+    optimiser = CRFOptimiser(conv_tol="loose", maxiter=1)
     optimiser._coords = CartesianCoordinates(np.array([]))
     optimiser._take_step_within_trust_radius(np.array([]))
 
@@ -163,7 +158,7 @@ def test_xtb_opt_with_distance_constraint():
 
     assert np.isclose(water.distance(0, 1), 0.99, atol=0.01)
 
-    CRFOptimiser.optimise(species=water, method=XTB(), etol=1e-6)
+    CRFOptimiser.optimise(species=water, method=XTB(), conv_tol="loose")
 
     # Optimisation should generate an O-H distance *very* close to 1.1 Å
     assert np.isclose(water.distance(0, 1).to("Å"), 1.1, atol=1e-4)
@@ -259,7 +254,7 @@ def test_xtb_opt_with_two_distance_constraint():
         (0, 2): val.Distance(1.0, units="Å"),
     }
 
-    opt = CRFOptimiser(maxiter=1, gtol=1e-5, etol=1e-5, init_alpha=0.1)
+    opt = CRFOptimiser(maxiter=1, conv_tol="loose", init_alpha=0.1)
     opt._species = water
     opt._method = XTB()
     opt._initialise_run()
@@ -279,7 +274,7 @@ def test_xtb_opt_with_two_distance_constraint():
             water.distance(i, j) - 1.0
         )
 
-    opt = CRFOptimiser(maxiter=10, gtol=1e-5, etol=1e-5, init_alpha=0.1)
+    opt = CRFOptimiser(maxiter=10, conv_tol="loose", init_alpha=0.1)
     opt.run(species=water, method=XTB())
     assert opt.converged
 
@@ -351,7 +346,7 @@ def test_optimise_linear_molecule():
     # the two H-C-C angles are almost linear
     assert val.Angle(170, "deg") < mol.angle(0, 1, 3) < val.Angle(176, "deg")
     assert val.Angle(170, "deg") < mol.angle(2, 0, 1) < val.Angle(176, "deg")
-    opt = CRFOptimiser(maxiter=10, gtol=1e-4, etol=1e-5)
+    opt = CRFOptimiser(maxiter=10, conv_tol="loose")
     opt.run(mol, XTB())
     assert opt.converged
     assert mol.angle(0, 1, 3) > val.Angle(179, "deg")
@@ -365,7 +360,7 @@ def test_optimise_linear_bend_with_ref():
     # the Fe-C-O angle are manually deviated
     assert val.Angle(170, "deg") < mol.angle(2, 3, 4) < val.Angle(176, "deg")
     # large molecule so allow few iters, no need to converge fully
-    opt = CRFOptimiser(maxiter=15, gtol=1e-5, etol=1e-6)
+    opt = CRFOptimiser(maxiter=10, conv_tol="loose")
     opt.run(mol, XTB())
     assert mol.angle(2, 3, 4) > val.Angle(178, "deg")
 
@@ -375,7 +370,7 @@ def test_optimise_linear_bend_with_ref():
 def test_optimise_chain_dihedrals():
     mol = cumulene_mol()
     assert abs(mol.dihedral(6, 3, 4, 8)) < val.Angle(40, "deg")
-    opt = CRFOptimiser(maxiter=20, gtol=1e-4, etol=1e-4)
+    opt = CRFOptimiser(maxiter=20, conv_tol="loose")
     opt.run(mol, XTB())
     # 5-C chain, should be close to 90 degrees
     assert abs(mol.dihedral(6, 3, 4, 8)) > val.Angle(85, "deg")
@@ -413,7 +408,7 @@ def test_composite_bond_constraint():
     # current sum of distances ~ 4.7 A, difference ~ 0.07
     assert np.isclose(mol.distance(1, 11), 2.385, rtol=1e-4)
     assert np.isclose(mol.distance(2, 7), 2.315, rtol=1e-4)
-    opt = CRFOptimiser(maxiter=10, gtol=1e-3, etol=1e-4, extra_prims=[constr])
+    opt = CRFOptimiser(maxiter=10, conv_tol="loose", extra_prims=[constr])
     opt.run(mol, method=XTB())
     # sum of distances should be ~4.6 A
     assert np.isclose(mol.distance(1, 11) + mol.distance(2, 7), 4.6, rtol=1e-4)
@@ -425,12 +420,10 @@ def test_trust_radius_limits():
     import autode.opt.optimisers.crfo
 
     max_lim = autode.opt.optimisers.crfo.MAX_TRUST
-    opt = CRFOptimiser(
-        maxiter=10, gtol=1e-3, etol=1e-4, init_trust=max_lim + 0.1
-    )
+    opt = CRFOptimiser(maxiter=10, conv_tol="loose", init_trust=max_lim + 0.1)
     assert np.isclose(opt.alpha, max_lim)
     min_lim = autode.opt.optimisers.crfo.MIN_TRUST
     opt = CRFOptimiser(
-        maxiter=10, gtol=1e-3, etol=1e-4, init_trust=min_lim - 0.001
+        maxiter=10, conv_tol="loose", init_trust=min_lim - 0.001
     )
     assert np.isclose(opt.alpha, min_lim)
diff --git a/tests/test_opt/test_opt.py b/tests/test_opt/test_opt.py
index 5f53833..3cbea6d 100644
--- a/tests/test_opt/test_opt.py
+++ b/tests/test_opt/test_opt.py
@@ -16,7 +16,11 @@ from .molecules import h2, methane_mol, h_atom
 from .setup import Method
 from autode.utils import NumericStringDict
 from autode.opt.coordinates import CartesianCoordinates
-from autode.opt.optimisers.base import OptimiserHistory, NullOptimiser
+from autode.opt.optimisers.base import (
+    OptimiserHistory,
+    NullOptimiser,
+    ConvergenceParams,
+)
 from autode.opt.optimisers.steepest_descent import (
     CartesianSDOptimiser,
     DIC_SD_Optimiser,
@@ -25,7 +29,7 @@ from autode.opt.optimisers.steepest_descent import (
 
 def sample_cartesian_optimiser():
     return CartesianSDOptimiser(
-        maxiter=1, gtol=GradientRMS(0.1), etol=PotentialEnergy(0.1)
+        maxiter=1, conv_tol=ConvergenceParams(abs_d_e=0.1, rms_g=0.1)
     )
 
 
@@ -39,27 +43,85 @@ def test_optimiser_construct():
         sample_cartesian_optimiser().run(species=methane_mol(), method=None)
 
     # Optimiser needs valid arguments
-    with pytest.raises(ValueError):
+    with pytest.raises(
+        ValueError, match="must be able to run at least one step"
+    ):
+        _ = CartesianSDOptimiser(maxiter=0, conv_tol="normal")
+
+    with pytest.raises(
+        ValueError, match="Value of abs_d_e should be positive"
+    ):
         _ = CartesianSDOptimiser(
-            maxiter=0, gtol=GradientRMS(0.1), etol=PotentialEnergy(0.1)
+            maxiter=1, conv_tol=ConvergenceParams(abs_d_e=-0.1, rms_g=0.1)
         )
 
-    with pytest.raises(ValueError):
+    with pytest.raises(ValueError, match="Value of rms_g should be positive"):
         _ = CartesianSDOptimiser(
-            maxiter=1, gtol=GradientRMS(-0.1), etol=PotentialEnergy(0.1)
+            maxiter=1, conv_tol=ConvergenceParams(abs_d_e=0.1, rms_g=-0.1)
         )
 
-    with pytest.raises(ValueError):
-        _ = CartesianSDOptimiser(
-            maxiter=1, gtol=GradientRMS(-0.1), etol=PotentialEnergy(0.1)
-        )
+    with pytest.raises(ValueError, match="Unknown preset convergence"):
+        _ = CartesianSDOptimiser(maxiter=1, conv_tol="unknown")
 
-    with pytest.raises(ValueError):
+    # should be able to set convergence through setter
+    opt = CartesianSDOptimiser(maxiter=1, conv_tol="loose")
+    opt.conv_tol = "normal"
+    with pytest.raises(ValueError, match="Unknown preset convergence"):
+        opt.conv_tol = "unknown"
+
+    # at least RMS g convergence criteria has to be defined
+    with pytest.raises(
+        ValueError, match="RMS gradient criteria has to be defined"
+    ):
         _ = CartesianSDOptimiser(
-            maxiter=1, gtol=GradientRMS(0.1), etol=PotentialEnergy(-0.1)
+            maxiter=1, conv_tol=ConvergenceParams(abs_d_e=0.1)
         )
 
 
+def test_optimiser_convergence(caplog):
+    opt = CartesianSDOptimiser(
+        maxiter=10,
+        conv_tol=ConvergenceParams(
+            abs_d_e=0.01, rms_g=0.01, max_g=0.01, rms_s=0.01, max_s=0.01
+        ),
+    )
+    coords1 = CartesianCoordinates(np.arange(6, dtype=float))
+    opt._species = Molecule(smiles="N#N")
+    opt._coords = coords1
+    opt._coords.g = np.random.random(6)
+    opt._coords.e = PotentialEnergy(0.1, "Ha")
+
+    # grad + energy < 1/2 + step is < * 3
+    coords2 = coords1 + 0.02
+    coords2.g = np.array([0.004] * 6)
+    coords2.e = PotentialEnergy(0.1 - 0.004, "Ha")
+    opt._coords = coords2
+    with caplog.at_level("WARNING"):
+        assert opt.converged
+    assert "Overachieved gradient and energy" in caplog.text
+    assert "reasonable convergence on step size" in caplog.text
+    caplog.clear()
+    # grad ~ 1/10, dE < *1.5, step < * 2
+    coords2 = coords1 + 0.014
+    coords2.g = np.array([0.0009] * 6)
+    coords2.e = PotentialEnergy(0.1 - 0.015)
+    opt._history._memory[-1] = coords2
+    with caplog.at_level("WARNING"):
+        assert opt.converged
+    assert "Gradient is one order of magnitude below" in caplog.text
+    assert "other parameter(s) are almost converged"
+    caplog.clear()
+    # step achieved, grad ~ 0.7,  dE < * 3
+    coords2 = coords1 + 0.009
+    coords2.g = np.array([0.006] * 6)
+    coords2.e = PotentialEnergy(0.1 - 0.025)
+    opt._history._memory[-1] = coords2
+    with caplog.at_level("WARNING"):
+        assert opt.converged
+    assert "Everything except energy has been converged" in caplog.text
+    assert "Reasonable convergence on energy" in caplog.text
+
+
 def test_initialise_species_and_method():
     optimiser = sample_cartesian_optimiser()
 
@@ -79,63 +141,6 @@ def test_coords_set():
         optimiser._coords = "a"
 
 
-def test_abs_diff_e():
-    # Define a intermediate optimiser state with two sets of coordinates
-    optimiser = sample_cartesian_optimiser()
-    optimiser._history.add(CartesianCoordinates([0.0, 1.0]))
-    optimiser._history.add(CartesianCoordinates([0.0, 1.1]))
-
-    # 2nd iteration for a history of two, indexed from 0
-    assert optimiser.iteration == 1
-
-    # without defined energies |E_0 - E_1| cannot be calculated
-    assert optimiser._abs_delta_e
-
-    # but can be if both structures have a potential energy
-    optimiser._history.final.e = PotentialEnergy(-1.0)
-    optimiser._history.penultimate.e = PotentialEnergy(-1.1)
-
-    diff_e = optimiser._abs_delta_e
-    assert isinstance(diff_e, PotentialEnergy)
-
-    assert np.isclose(diff_e, 0.1, atol=1e-6)
-
-
-def test_g_norm():
-    optimiser = sample_cartesian_optimiser()
-
-    # With no coordinates the norm of the gradient is infinity
-    assert optimiser._coords is None
-    assert not np.isfinite(optimiser._g_norm)
-
-    # Likewise if the gradient is unset
-    optimiser._coords = CartesianCoordinates([1.0, 0.0, 0.0])
-    assert optimiser._coords.g is None
-    assert not np.isfinite(optimiser._g_norm)
-
-
-def test_optimiser_h_update():
-    optimiser = sample_cartesian_optimiser()
-
-    # Remove any possible updater type
-    optimiser._hessian_update_types = []
-
-    c1 = CartesianCoordinates([1.0, 0.0, 0.0])
-    c1.h = np.eye(3)
-
-    optimiser._history.add(c1)
-
-    c2 = CartesianCoordinates([1.1, 0.0, 0.0])
-    c2.h = np.eye(3)
-
-    optimiser._history.add(c2)
-
-    # and try and update the (inverse) hessian, which is impossible without
-    # an updater
-    with pytest.raises(Exception):
-        _ = optimiser._updated_h_inv()
-
-
 def test_history():
     optimiser = sample_cartesian_optimiser()
     assert optimiser.iteration < 1
@@ -162,14 +167,14 @@ def test_xtb_h2_cart_opt():
 
 @work_in_tmp_dir()
 @requires_working_xtb_install
-def test_xtb_h2_cart_opt():
+def test_xtb_h2_cart_opt_2():
     optimiser = CartesianSDOptimiser(
-        maxiter=2,
-        gtol=GradientRMS(0.01),
-        etol=PotentialEnergy(1e-3),
+        maxiter=2, conv_tol=ConvergenceParams(abs_d_e=1e-3, rms_g=0.01)
     )
     optimiser._species = h2()
     optimiser._coords = CartesianCoordinates(optimiser._species.coordinates)
+    optimiser._species.single_point(XTB())
+    optimiser._coords.e = optimiser._species.energy
 
     assert not optimiser.converged
 
@@ -190,8 +195,7 @@ def test_xtb_h2_dic_opt():
     optimiser = DIC_SD_Optimiser(
         step_size=2.5,
         maxiter=10,
-        gtol=GradientRMS(0.01),
-        etol=PotentialEnergy(0.0001),
+        conv_tol=ConvergenceParams(abs_d_e=1e-4, rms_g=0.01),
     )
 
     mol = h2()
@@ -223,8 +227,7 @@ def test_callback_function():
         maxiter=1,
         callback=func,
         callback_kwargs={"m": mol},
-        gtol=GradientRMS(0.1),
-        etol=PotentialEnergy(0.1),
+        conv_tol=ConvergenceParams(rms_g=0.1, abs_d_e=0.1),
     )
 
     optimiser.run(species=mol, method=Method())
@@ -234,7 +237,7 @@ def test_callback_function():
 def test_last_energy_change_with_no_steps():
     mol = h2()
     optimiser = HarmonicPotentialOptimiser(
-        maxiter=2, gtol=GradientRMS(999), etol=PotentialEnergy(999)
+        maxiter=2, conv_tol=ConvergenceParams(abs_d_e=999, rms_g=999)
     )
 
     optimiser.run(mol, method=Method())
@@ -270,7 +273,7 @@ class UnconvergedHarmonicPotentialOptimiser(CartesianSDOptimiser):
 
 def test_last_energy_change_less_than_two_steps():
     optimiser = ConvergedHarmonicPotentialOptimiser(
-        maxiter=2, gtol=GradientRMS(999), etol=PotentialEnergy(999)
+        maxiter=2, conv_tol=ConvergenceParams(abs_d_e=999, rms_g=999)
     )
 
     coords = CartesianCoordinates(np.zeros(1))
@@ -296,11 +299,7 @@ def test_hessian_is_not_recalculated_if_present():
     mol = h2()
     xtb = XTB()
 
-    optimiser = CartesianSDOptimiser(
-        maxiter=1,
-        gtol=GradientRMS(0.01),
-        etol=PotentialEnergy(1e-3),
-    )
+    optimiser = CartesianSDOptimiser(maxiter=1, conv_tol="loose")
     optimiser.run(species=mol, method=xtb, n_cores=1)
 
     mol.calc_hessian(method=xtb)
@@ -314,11 +313,7 @@ def test_hessian_is_not_recalculated_if_present():
 @work_in_tmp_dir()
 @requires_working_xtb_install
 def test_multiple_optimiser_saves_overrides_not_append():
-    optimiser = CartesianSDOptimiser(
-        maxiter=2,
-        gtol=GradientRMS(0.01),
-        etol=PotentialEnergy(1e-3),
-    )
+    optimiser = CartesianSDOptimiser(maxiter=2, conv_tol="loose")
     optimiser.run(method=XTB(), species=h2(), name="tmp.zip")
 
     assert os.path.isfile("tmp.zip")
@@ -327,11 +322,7 @@ def test_multiple_optimiser_saves_overrides_not_append():
 
     old_n_coords = sum([1 for name in names if name.startswith("coords_")])
 
-    optimiser = CartesianSDOptimiser(
-        maxiter=2,
-        gtol=GradientRMS(0.01),
-        etol=PotentialEnergy(1e-3),
-    )
+    optimiser = CartesianSDOptimiser(maxiter=2, conv_tol="loose")
     optimiser.run(method=XTB(), species=h2(), name="tmp.zip")
     # the file "tmp.zip" should be overwritten by new optimiser
     with zipfile.ZipFile("tmp.zip") as file:
@@ -344,7 +335,7 @@ def test_multiple_optimiser_saves_overrides_not_append():
 @work_in_tmp_dir()
 def test_optimiser_plotting_sanity_checks(caplog):
     mol = Molecule(smiles="N#N")
-    opt = CartesianSDOptimiser(maxiter=10, gtol=1e-3, etol=1e-3)
+    opt = CartesianSDOptimiser(maxiter=10, conv_tol="loose")
     coords1 = CartesianCoordinates(mol.coordinates)
     coords1.e = PotentialEnergy(0.1, "Ha")
     coords1.update_g_from_cart_g(
@@ -378,7 +369,7 @@ def test_optimiser_plotting_sanity_checks(caplog):
 def test_optimiser_print_geometries(caplog):
     mol = Molecule(smiles="C=C", name="mymolecule")
     coords1 = CartesianCoordinates(mol.coordinates)
-    opt = CartesianSDOptimiser(maxiter=20, gtol=1e-3, etol=1e-4)
+    opt = CartesianSDOptimiser(maxiter=20, conv_tol="loose")
     opt._coords = coords1
     # cannot print geom without species
     with pytest.raises(AssertionError):
diff --git a/tests/test_opt/test_prfo.py b/tests/test_opt/test_prfo.py
index 09cbc2f..7a915a4 100644
--- a/tests/test_opt/test_prfo.py
+++ b/tests/test_opt/test_prfo.py
@@ -90,7 +90,7 @@ def test_mode_following():
             Atom("H", -2.10961, 4.17548, 1.25945),
         ],
     )
-    opt = PRFOptimiser(maxiter=10, etol=1e-4, gtol=1e-3, imag_mode_idx=0)
+    opt = PRFOptimiser(maxiter=10, conv_tol="normal", imag_mode_idx=0)
     opt._species = mol
     opt._method = xtb
     opt._n_cores = Config.n_cores
diff --git a/tests/test_opt/test_qa.py b/tests/test_opt/test_qa.py
index 2afe12e..fa6770f 100644
--- a/tests/test_opt/test_qa.py
+++ b/tests/test_opt/test_qa.py
@@ -16,7 +16,7 @@ datazip = os.path.join(here, "data", "opt.zip")
 @work_in_zipped_dir(datazip)
 def test_trm_step():
     mol = Molecule("opt-test.xyz")
-    opt = QAOptimiser(maxiter=10, gtol=0.001, etol=0.001, init_trust=0.1)
+    opt = QAOptimiser(maxiter=10, conv_tol="loose", init_trust=0.1)
     opt._species = mol
     opt._build_internal_coordinates()
     assert isinstance(opt._coords, DICWithConstraints)
@@ -42,7 +42,7 @@ def test_trust_update():
         Atom("H", 0.8261, -0.1812, 0.0000),
     ]
     water = Molecule(atoms=water_atoms)
-    opt = QAOptimiser(maxiter=10, gtol=1e-3, etol=1e-4, init_trust=init_trust)
+    opt = QAOptimiser(maxiter=10, conv_tol="loose", init_trust=init_trust)
 
     opt._species = water.copy()
     opt._method = XTB()
diff --git a/tests/test_opt/test_rfo.py b/tests/test_opt/test_rfo.py
index d7907a3..ced3d13 100644
--- a/tests/test_opt/test_rfo.py
+++ b/tests/test_opt/test_rfo.py
@@ -2,6 +2,7 @@ import numpy as np
 from autode.species.molecule import Molecule
 from autode.wrappers.XTB import XTB
 from autode.utils import work_in_tmp_dir
+from autode.opt.optimisers import ConvergenceParams
 from autode.opt.optimisers.rfo import RFOptimiser
 from autode.opt.coordinates import CartesianCoordinates
 from ..testutils import requires_working_xtb_install
@@ -20,11 +21,10 @@ class TestRFOOptimiser2D(RFOptimiser):
         init_x,
         init_y,
         maxiter=30,
-        etol=1e-4,
-        gtol=1e-3,
+        conv_tol=ConvergenceParams(abs_d_e=1e-4, rms_g=1e-3),
         **kwargs,
     ):
-        super().__init__(maxiter=maxiter, etol=etol, gtol=gtol, **kwargs)
+        super().__init__(maxiter=maxiter, conv_tol=conv_tol, **kwargs)
 
         init_arr = np.array([init_x, init_y])
         self._coords = CartesianCoordinates(init_arr)