Repository: duartegroup/autodE
Pull request #354 (2024-09-26 08:47:14)

V1.4.4
@shoubhikraj did you want to add anything else into this release?


Patch:
diff --git a/autode/atoms.py b/autode/atoms.py
index 6e88d8f..675f7a9 100644
--- a/autode/atoms.py
+++ b/autode/atoms.py
@@ -804,7 +804,7 @@ class Atoms(list):
         Raises:
             (IndexError): If i or j are not present
         """
-        return self[j].coord - self[i].coord
+        return np.asarray(self[j].coord - self[i].coord)
 
     def nvector(self, i: int, j: int) -> np.ndarray:
         """
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/hessians.py b/autode/hessians.py
index bbe00fb..c8ef05b 100644
--- a/autode/hessians.py
+++ b/autode/hessians.py
@@ -353,7 +353,7 @@ class Hessian(ValueArray):
             (list(autode.values.Frequency)):
         """
 
-        nus = np.sqrt(np.complex_(lambdas)) / (
+        nus = np.sqrt(np.complex128(lambdas)) / (
             2.0 * np.pi * Constants.ang_to_m * Constants.c_in_cm
         )
         nus *= self._freq_scale_factor
diff --git a/autode/opt/coordinates/base.py b/autode/opt/coordinates/base.py
index af28433..8b26086 100644
--- a/autode/opt/coordinates/base.py
+++ b/autode/opt/coordinates/base.py
@@ -228,13 +228,14 @@ class OptCoordinates(ValueArray, ABC):
         assert isinstance(old_coords, OptCoordinates), "Wrong type!"
         assert old_coords._h is not None
         assert old_coords._g is not None
+        idxs = self.active_mol_indexes
 
         for update_type in hessian_update_types:
             updater = update_type(
                 h=old_coords._h,
-                s=self.raw - old_coords.raw,
+                s=np.array(self) - np.array(old_coords),
                 y=self._g - old_coords._g,
-                subspace_idxs=old_coords.indexes,
+                subspace_idxs=idxs,
             )
 
             if not updater.conditions_met:
@@ -250,6 +251,79 @@ class OptCoordinates(ValueArray, ABC):
             "Could not update the Hessian - no suitable update strategies"
         )
 
+    @property
+    def rfo_shift(self) -> float:
+        """
+        Get the RFO diagonal shift factor λ for the molecular Hessian that
+        can be applied (H - λI) to obtain the RFO downhill step. The shift
+        is only calculated in active subspace
+
+        Returns:
+            (float): The shift parameter
+        """
+        assert self._h is not None
+        # ignore constraint modes
+        n, _ = self._h.shape
+        idxs = self.active_mol_indexes
+        hess = self._h[:, idxs][idxs, :]
+        grad = self._g[idxs]
+
+        h_n, _ = hess.shape
+        # form the augmented Hessian in active subspace
+        aug_h = np.zeros(shape=(h_n + 1, h_n + 1))
+
+        aug_h[:h_n, :h_n] = hess
+        aug_h[-1, :h_n] = grad
+        aug_h[:h_n, -1] = grad
+
+        # first non-zero eigenvalue
+        aug_h_lmda = np.linalg.eigvalsh(aug_h)
+        rfo_lmda = aug_h_lmda[0]
+        assert abs(rfo_lmda) > 1.0e-10
+        return rfo_lmda
+
+    @property
+    def min_eigval(self) -> float:
+        """
+        Obtain the minimum eigenvalue of the molecular Hessian in
+        the active space
+
+        Returns:
+            (float): The minimum eigenvalue
+        """
+        assert self._h is not None
+        n, _ = self._h.shape
+        idxs = self.active_mol_indexes
+        hess = self._h[:, idxs][idxs, :]
+
+        eigvals = np.linalg.eigvalsh(hess)
+        assert abs(eigvals[0]) > 1.0e-10
+        return eigvals[0]
+
+    def pred_quad_delta_e(self, new_coords: np.ndarray) -> float:
+        """
+        Calculate the estimated change in energy at the new coordinates
+        based on the quadratic model (i.e. second order Taylor expansion)
+
+        Args:
+            new_coords(np.ndarray): The new coordinates
+
+        Returns:
+            (float): The predicted change in energy
+        """
+        assert self._g is not None and self._h is not None
+
+        step = np.array(new_coords) - np.array(self)
+
+        idxs = self.active_mol_indexes
+        step = step[idxs]
+        grad = self._g[idxs]
+        hess = self._h[:, idxs][idxs, :]
+
+        pred_delta = np.dot(grad, step)
+        pred_delta += 0.5 * np.linalg.multi_dot((step, hess, step))
+        return pred_delta
+
     def make_hessian_positive_definite(self) -> None:
         """
         Make the Hessian matrix positive definite by shifting eigenvalues
@@ -269,6 +343,38 @@ 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]:
+        """A list of indexes which are active in this coordinate set"""
+
+    @property
+    def active_mol_indexes(self) -> List[int]:
+        """Active indexes that are actually atomic coordinates in the molecule"""
+        return [i for i in self.active_indexes if i < len(self)]
+
+    @property
+    @abstractmethod
+    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 423539f..7b9c00f 100644
--- a/autode/opt/coordinates/cartesian.py
+++ b/autode/opt/coordinates/cartesian.py
@@ -1,5 +1,5 @@
 import numpy as np
-from typing import Optional, TYPE_CHECKING
+from typing import Optional, List, TYPE_CHECKING
 
 from autode.log import logger
 from autode.values import ValueArray
@@ -45,7 +45,7 @@ class CartesianCoordinates(OptCoordinates):
         Arguments:
             arr: Gradient array
         """
-        self.g = None if arr is None else np.array(arr).flatten()
+        self._g = None if arr is None else np.array(arr).flatten()
 
     def _update_h_from_cart_h(self, arr: Optional["Hessian"]) -> None:
         """
@@ -57,7 +57,24 @@ class CartesianCoordinates(OptCoordinates):
         Arguments:
             arr: Hessian matrix
         """
-        self.h = None if arr is None else np.array(arr)
+        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)))
+
+    @property
+    def inactive_indexes(self) -> List[int]:
+        return []
 
     def iadd(self, value: np.ndarray) -> OptCoordinates:
         return np.ndarray.__iadd__(self, value)
@@ -96,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 2d52778..e0a39d9 100644
--- a/autode/opt/coordinates/dic.py
+++ b/autode/opt/coordinates/dic.py
@@ -34,7 +34,7 @@ if TYPE_CHECKING:
     from autode.hessians import Hessian
 
 
-_max_back_transform_iterations = 20
+MAX_BACK_TRANSFORM_ITERS = 50
 
 
 class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
@@ -70,8 +70,10 @@ class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
         Returns:
             (np.ndarray): U
         """
-
-        lambd, u = np.linalg.eigh(primitives.G)
+        # calculate spectroscopic G matrix
+        B = primitives.get_B(x)
+        G = np.dot(B, B.T)
+        lambd, u = np.linalg.eigh(G)
 
         # Form a transform matrix from the primitive internals by removing the
         # redundant subspace comprised of small eigenvalues. This forms a set
@@ -124,8 +126,9 @@ class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
 
         dic.U = U  # Transform matrix primitives -> non-redundant
 
-        dic.B = np.matmul(U.T, primitives.B)
+        dic.B = np.matmul(U.T, primitives.get_B(x))
         dic.B_T_inv = np.linalg.pinv(dic.B)
+        dic._q = q.copy()
         dic._x = x.copy()
         dic.primitives = primitives
 
@@ -136,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
@@ -145,11 +152,11 @@ class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
             arr: Cartesian gradient array
         """
         if arr is None:
-            self._x.g, self.g = None, None
+            self._x.g, self._g = None, None
 
         else:
             self._x.g = arr.flatten()
-            self.g = np.matmul(self.B_T_inv.T, self._x.g)
+            self._g = np.matmul(self.B_T_inv.T, self._x.g)
 
         return None
 
@@ -162,14 +169,16 @@ class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
             arr: Cartesian Hessian matrix
         """
         if arr is None:
-            self._x.h, self.h = None, None
+            self._x.h, self._h = None, None
 
         else:
             self._x.h = arr
 
             # NOTE: This is not the full transformation as noted in
             # 10.1063/1.471864 only an approximate Hessian is required(?)
-            self.h = np.linalg.multi_dot((self.B_T_inv.T, arr, self.B_T_inv))
+            hess = np.linalg.multi_dot((self.B_T_inv.T, arr, self.B_T_inv))
+            assert self.h_or_h_inv_has_correct_shape(hess)
+            self._h = hess
 
         return None
 
@@ -213,22 +222,46 @@ class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
 
         # Initialise
         s_k, x_k = np.array(self, copy=True), self.to("cartesian").copy()
-        q_init = self.primitives(x_k)
+        q_init = self._q
         x_1 = self.to("cartesian") + np.matmul(self.B_T_inv, value)
 
         success = False
-        for i in range(1, _max_back_transform_iterations + 1):
+        rms_s = np.inf
+        # NOTE: J. Comput. Chem., 2013, 34, 1842 suggests if step size
+        # is larger than 0.5 bohr (= 0.2 Å), internal step can be halved
+        # for easier convergence (i.e. damp = 1/2)
+        if np.linalg.norm(value) > 0.2:
+            damp = 0.5
+        else:
+            damp = 1.0
+
+        # hybrid SIBT/IBT algorithm
+        for i in range(1, MAX_BACK_TRANSFORM_ITERS + 1):
             try:
-                x_k = x_k + np.matmul(self.B_T_inv, (s_new - s_k))
+                x_k = x_k + np.matmul(self.B_T_inv, damp * (s_new - s_k))
 
-                # Rebuild the primitives & DIC from the back-transformed Cartesians
+                # Rebuild the DIC from back-transformed Cartesians
                 q_k = self.primitives.close_to(x_k, q_init)
                 s_k = np.matmul(self.U.T, q_k)
-                self.B = np.matmul(self.U.T, self.primitives.B)
-                self.B_T_inv = np.linalg.pinv(self.B)
 
+                # Rebuild the B matrix every 10 steps
+                if i % 10 == 0:
+                    self.B = np.matmul(self.U.T, self.primitives.get_B(x_k))
+                    self.B_T_inv = np.linalg.pinv(self.B)
+
+                rms_s_old = rms_s
                 rms_s = np.sqrt(np.mean(np.square(s_k - s_new)))
 
+                # almost converged, turn off damping
+                if rms_s < 1e-6:
+                    damp = 1.0
+                # RMS going down, reduce damping
+                elif rms_s < rms_s_old and i > 1:
+                    damp = min(1.2 * damp, 1.0)
+                # RMS going up, increase damping
+                elif rms_s > rms_s_old:
+                    damp = max(0.7 * damp, 0.1)
+
             # for ill-conditioned primitives, there might be math error
             except ArithmeticError:
                 break
@@ -245,7 +278,7 @@ class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
         else:
             logger.warning(
                 f"Failed to transform in {i} cycles. "
-                f"Final RMS(s) = {rms_s:.8f}"
+                + f"Final RMS(s) = {rms_s:.8f}"
             )
             x_k = x_1
             if not self.allow_unconverged_back_transform:
@@ -253,11 +286,13 @@ class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
                     "DIC->Cart iterative back-transform did not converge"
                 )
 
-        s_k = np.matmul(self.U.T, self.primitives(x_k))
-        self.B = np.matmul(self.U.T, self.primitives.B)
+        q_k = self.primitives.close_to(x_k, q_init)
+        s_k = np.matmul(self.U.T, q_k)
+        self.B = np.matmul(self.U.T, self.primitives.get_B(x_k))
         self.B_T_inv = np.linalg.pinv(self.B)
 
         self[:] = s_k
+        self._q = q_k
         self._x = x_k
 
         return self
@@ -292,7 +327,7 @@ class DICWithConstraints(DIC):
 
         arr = super().__new__(cls, input_array)
 
-        arr._lambda = None  # Additional lagrangian multipliers
+        arr._lambda = np.array([])  # Additional lagrangian multipliers
         return arr
 
     def __array_finalize__(self, obj: "OptCoordinates") -> None:
@@ -301,11 +336,28 @@ class DICWithConstraints(DIC):
         self._lambda = getattr(obj, "_lambda", None)
         return
 
-    def zero_lagrangian_multipliers(self) -> None:
-        r"""Zero all \lambda_i"""
+    @classmethod
+    def from_cartesian(
+        cls,
+        x: "CartesianCoordinates",
+        primitives: Optional[PIC] = None,
+    ) -> "DICWithConstraints":
+        """
+        Generate delocalised internal coordinates with constraints
+        with the Lagrangian multipliers initialised as zeroes
 
-        self._lambda = np.zeros(shape=(self.n_constraints,))
-        return None
+        Args:
+            x: Cartesian coordinates
+
+            primitives: Primitive internal coordinates. If undefined then use
+                        all pairwise inverse distances
+
+        Returns:
+            (DICWithConstraints): DIC with constraints
+        """
+        dic = super().from_cartesian(x=x, primitives=primitives)
+        dic._lambda = np.zeros(shape=(dic.n_constraints,))
+        return dic
 
     @property
     def raw(self) -> np.ndarray:
@@ -355,88 +407,98 @@ class DICWithConstraints(DIC):
 
         return [i for i in range(n + m) if i not in self.inactive_indexes]
 
-    def _update_g_from_cart_g(self, arr: Optional["Gradient"]) -> None:
-        r"""
-        Updates the gradient from a calculated Cartesian gradient, where
-        the gradient is that of the Lagrangian. Includes dL/d_λi terms where
-        λi is the i-th lagrangian multiplier.
+    def iadd(self, value: np.ndarray) -> "OptCoordinates":
+        """
+        Add a step in internal coordinates (along with Lagrange multipliers)
+        to this set of coordinates, and update the Cartesian coordinates
 
-        -----------------------------------------------------------------------
-        Arguments:
-            arr: Cartesian gradient array
+        Args:
+            value: Difference between current and new DICs, and the multipliers
         """
+        assert len(value) == len(self) + self.n_constraints
+        # separate the coordinates and the lagrange multipliers
+        if self.n_constraints > 0:
+            delta_lambda = value[-self.n_constraints :]
+            self._lambda += delta_lambda
+            delta_s = value[: -self.n_constraints]
+        else:
+            delta_s = value
 
-        if arr is None:
-            self._x.g, self.g = None, None
-            return
+        return super().iadd(delta_s)
 
-        assert self._lambda is not None, "Must have λ defined"
-
-        self._x.g = arr.flatten()
-        n = len(self)
-        m = self.n_constraints
+    @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
 
-        self.g = np.zeros(shape=(n + m,))
+    @property
+    def g(self):
+        """
+        Gradient of the energy, contains the Lagrangian dL/d_λi terms where
+        λi is the i-th lagrangian multiplier.
+        """
+        if self._g is None:
+            return None
 
-        # Set the first part dL/ds_i
-        self.g[:n] = np.matmul(self.B_T_inv.T, self._x.g)
+        n, m = len(self), self.n_constraints
+        arr = np.zeros(shape=(n + m,))
+        arr[:n] = self._g
 
+        # constrained gradient terms
         for i in range(m):
-            self.g[n - m + i] -= self._lambda[i] * 1  # λ dC_i/ds_i
+            arr[n - m + i] -= self._lambda[i] * 1  # λ dC_i/ds_i
 
-        # and the final dL/dλ_i
+        # final dL/dλ_i
         c = self.constrained_primitives
         for i in range(m):
-            self.g[n + i] = -c[i].delta(self._x)  # C_i(x) = Z - Z_ideal
+            arr[n + i] = -c[i].delta(self._x)  # C_i(x) = Z - Z_ideal
 
-        return None
+        return arr
 
-    def _update_h_from_cart_h(self, arr: Optional["Hessian"]) -> None:
-        """
-        Update the DIC Hessian matrix from a Cartesian one
+    @g.setter
+    def g(self, value):
+        """Setting g is not allowed with constraints"""
+        raise RuntimeError("Cannot set gradient with constraints enabled")
 
-        -----------------------------------------------------------------------
-        Arguments:
-            arr: Cartesian Hessian matrix
+    @property
+    def h(self):
         """
+        The Hessian matrix, containing Lagrangian constraint terms
 
-        if arr is None:
-            self._x.h, self.h = None, None
-
-        else:
-            self._x.h = arr
-            n = len(self)
-            m = self.n_constraints
-
-            self.h = np.zeros(shape=(n + m, n + m))
-
-            # Fill in the upper left corner with d^2L/ds_i ds_j
-            # where the second derivative of the constraint is zero
-            self.h[:n, :n] = np.linalg.multi_dot(
-                (self.B_T_inv.T, arr, self.B_T_inv)
-            )
+        Returns:
+            (np.ndarray):
+        """
+        if self._h is None:
+            return None
 
-            # and the d^2L/ds_i dλ_i = -dC_i/ds_i = -1
-            #         d^2L/dλ_i dλ_j = 0
+        n, m = len(self), self.n_constraints
+        arr = np.zeros(shape=(n + m, n + m))
 
-            for i in range(m):
-                self.h[n + i, :] = self.h[:, n + i] = 0.0
+        # Upper left corner is d^2L/ds_i ds_j
+        arr[:n, :n] = self._h
 
-            for i in range(m):
-                self.h[n - m + i, n + i] = self.h[n + i, n - m + i] = -1.0
+        # and the d^2L/ds_i dλ_i = -dC_i/ds_i = -1
+        #         d^2L/dλ_i dλ_j = 0
+        for i in range(m):
+            arr[n + i, :] = arr[:, n + i] = 0.0
 
-        return None
+        for i in range(m):
+            arr[n - m + i, n + i] = arr[n + i, n - m + i] = -1.0
 
-    def h_or_h_inv_has_correct_shape(self, arr: Optional[np.ndarray]):
-        """Does a Hessian or its inverse have the correct shape?"""
-        if arr is None:
-            return True  # None is always valid
+        return arr
 
-        n_rows, n_cols = arr.shape
-        return (
-            arr.ndim == 2
-            and n_rows == n_cols == len(self) + self.n_constraints
-        )
+    @h.setter
+    def h(self, value):
+        raise RuntimeError("Cannot set hessian when constraints are enabled")
 
     def update_lagrange_multipliers(self, arr: np.ndarray) -> None:
         """Update the lagrange multipliers by adding a set of values"""
diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 3de919b..d4b608b 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -24,6 +24,7 @@ from autode.opt.coordinates.primitives import (
     PrimitiveDummyLinearAngle,
     PrimitiveLinearAngle,
     PrimitiveDihedralAngle,
+    PrimitiveImproperDihedral,
     LinearBendType,
 )
 
@@ -47,9 +48,10 @@ class InternalCoordinates(OptCoordinates, ABC):  # lgtm [py/missing-equals]
         arr = super().__new__(cls, input_array, units="Å")
 
         arr._x = None
+        arr._q = None
         arr.primitives = None
 
-        for attr in ("_x", "primitives"):
+        for attr in ("_x", "primitives", "_q"):
             setattr(arr, attr, getattr(input_array, attr, None))
 
         return arr
@@ -58,7 +60,7 @@ class InternalCoordinates(OptCoordinates, ABC):  # lgtm [py/missing-equals]
         """See https://numpy.org/doc/stable/user/basics.subclassing.html"""
         OptCoordinates.__array_finalize__(self, obj)
 
-        for attr in ("_x", "primitives"):
+        for attr in ("_x", "primitives", "_q"):
             setattr(self, attr, getattr(obj, attr, None))
 
         return
@@ -90,8 +92,6 @@ class PIC(list, ABC):
         """
         super().__init__(args)
 
-        self._B: Optional[np.ndarray] = None
-
         if not self._are_all_primitive_coordinates(args):
             raise ValueError(
                 "Cannot construct primitive internal coordinates "
@@ -106,28 +106,11 @@ class PIC(list, ABC):
             super().append(item)
 
     def append(self, item: Primitive) -> None:
-        """Append an item to this set of primitives"""
+        """Appending directly is not allowed, use add() instead"""
         raise NotImplementedError(
             "Please use PIC.add() to add new primitives to the set"
         )
 
-    @property
-    def B(self) -> np.ndarray:
-        """Wilson B matrix"""
-
-        if self._B is None:
-            raise AttributeError(
-                f"{self} had no B matrix. Please calculate "
-                f"the value of the primitives to determine B"
-            )
-
-        return self._B
-
-    @property
-    def G(self) -> np.ndarray:
-        """Spectroscopic G matrix as the symmetrised Wilson B matrix"""
-        return np.dot(self.B, self.B.T)
-
     @classmethod
     def from_cartesian(
         cls,
@@ -145,7 +128,6 @@ class PIC(list, ABC):
         """Populate Primitive-s used in the construction of set"""
 
         q = self._calc_q(x)
-        self._calc_B(x)
 
         return q
 
@@ -159,7 +141,6 @@ class PIC(list, ABC):
         assert len(self) == len(other) and isinstance(other, np.ndarray)
 
         q = self._calc_q(x)
-        self._calc_B(x)
 
         for i, primitive in enumerate(self):
             if isinstance(primitive, PrimitiveDihedralAngle):
@@ -193,7 +174,7 @@ class PIC(list, ABC):
     def _populate_all(self, x: np.ndarray) -> None:
         """Populate primitives from an array of cartesian coordinates"""
 
-    def _calc_B(self, x: np.ndarray) -> None:
+    def get_B(self, x: np.ndarray) -> np.ndarray:
         """Calculate the Wilson B matrix"""
 
         if len(self) == 0:
@@ -209,8 +190,7 @@ class PIC(list, ABC):
         for i, primitive in enumerate(self):
             B[i] = primitive.derivative(x=cart_coords)
 
-        self._B = B
-        return None
+        return B
 
     @staticmethod
     def _are_all_primitive_coordinates(args: tuple) -> bool:
@@ -397,7 +377,7 @@ class AnyPIC(PIC):
                         mol, m, o, n, bonded=True
                     )
                     if r is not None:
-                        self.add(PrimitiveDihedralAngle(m, r, o, n))
+                        self.add(PrimitiveImproperDihedral(m, r, o, n))
                         continue
 
                     # Otherwise, we use a nearby (< 4.0 A) reference atom to
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 52efca8..530e3e7 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -1,5 +1,5 @@
 import numpy as np
-
+import itertools
 from abc import ABC, abstractmethod
 from enum import Enum
 from typing import Tuple, TYPE_CHECKING, List, Optional
@@ -31,7 +31,7 @@ def _get_3d_vecs_from_atom_idxs(
         deriv_order: Order of derivatives for initialising variables
 
     Returns:
-        (list[VectorHyperDual]): A list of differentiable variables
+        (list[DifferentiableVector3D]): A list of differentiable variables
     """
     assert all(isinstance(idx, int) and idx >= 0 for idx in args)
     # get positions in the flat Cartesian array
@@ -65,7 +65,7 @@ class Primitive(ABC):
     @abstractmethod
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
-    ):
+    ) -> VectorHyperDual:
         """
         The function that performs the main evaluation of the PIC,
         and optionally returns derivative or second derivatives.
@@ -234,7 +234,7 @@ class PrimitiveInverseDistance(_DistanceFunction):
 
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
-    ) -> "VectorHyperDual":
+    ) -> VectorHyperDual:
         """1 / |x_i - x_j|"""
         vec_i, vec_j = _get_3d_vecs_from_atom_idxs(
             self.i, self.j, x=x, deriv_order=deriv_order
@@ -256,7 +256,7 @@ class PrimitiveDistance(_DistanceFunction):
 
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
-    ) -> "VectorHyperDual":
+    ) -> VectorHyperDual:
         """|x_i - x_j|"""
         vec_i, vec_j = _get_3d_vecs_from_atom_idxs(
             self.i, self.j, x=x, deriv_order=deriv_order
@@ -318,14 +318,16 @@ class PrimitiveBondAngle(Primitive):
 
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
-    ):
+    ) -> VectorHyperDual:
         """m - o - n angle"""
         vec_m, vec_o, vec_n = _get_3d_vecs_from_atom_idxs(
             self.m, self.o, self.n, x=x, deriv_order=deriv_order
         )
         u = vec_m - vec_o
         v = vec_n - vec_o
-        return DifferentiableMath.acos(u.dot(v) / (u.norm() * v.norm()))
+        res = DifferentiableMath.acos(u.dot(v) / (u.norm() * v.norm()))
+        assert isinstance(res, VectorHyperDual)
+        return res
 
     def __repr__(self):
         return f"Angle({self.m}-{self.o}-{self.n})"
@@ -385,7 +387,7 @@ class PrimitiveDihedralAngle(Primitive):
 
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
-    ) -> "VectorHyperDual":
+    ) -> VectorHyperDual:
         """Dihedral m-o-p-n"""
         # https://en.wikipedia.org/wiki/Dihedral_angle#In_polymer_physics
         _x = x.ravel()
@@ -408,6 +410,13 @@ class PrimitiveDihedralAngle(Primitive):
         return f"Dihedral({self.m}-{self.o}-{self.p}-{self.n})"
 
 
+class PrimitiveImproperDihedral(PrimitiveDihedralAngle):
+    """Out-of-Plan (improper) dihedral angles"""
+
+    def __repr__(self):
+        return f"ImproperDihedral({self.m}-{self.o}-{self.p}-{self.n})"
+
+
 class LinearBendType(Enum):
     """For linear angles, there are two orthogonal directions"""
 
@@ -440,7 +449,7 @@ class LinearAngleBase(Primitive, ABC):
         o_vec: DifferentiableVector3D,
         n_vec: DifferentiableVector3D,
         r_vec: DifferentiableVector3D,
-    ):
+    ) -> VectorHyperDual:
         """
         Evaluate the linear bend from the vector positions of the
         atoms involved in the angle m, o, n, and the reference
@@ -466,11 +475,13 @@ class LinearAngleBase(Primitive, ABC):
 
         # eq. (46) and (47) p 1074
         if self.axis == LinearBendType.BEND:
-            return u.dot(o_n) / o_n.norm()
+            res = u.dot(o_n) / o_n.norm()
         elif self.axis == LinearBendType.COMPLEMENT:
-            return u.dot(o_n.cross(o_m)) / (o_n.norm() * o_m.norm())
+            res = u.dot(o_n.cross(o_m)) / (o_n.norm() * o_m.norm())
         else:
             raise ValueError("Unknown axis for linear bend")
+        assert isinstance(res, VectorHyperDual)
+        return res
 
 
 class PrimitiveLinearAngle(LinearAngleBase):
@@ -478,7 +489,7 @@ class PrimitiveLinearAngle(LinearAngleBase):
 
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
-    ):
+    ) -> VectorHyperDual:
         """Linear Bend angle m-o-n against reference atom r"""
 
         _x = x.ravel()
@@ -529,6 +540,7 @@ class PrimitiveDummyLinearAngle(LinearAngleBase):
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ):
+        """Linear bend m-o-n against a dummy atom"""
         if self._vec_r is None:
             self._vec_r = self._get_dummy_atom(x)
 
@@ -542,3 +554,89 @@ class PrimitiveDummyLinearAngle(LinearAngleBase):
     def __repr__(self):
         axis_str = "B" if self.axis == LinearBendType.BEND else "C"
         return f"LinearBend{axis_str}({self.m}-{self.o}-{self.n}, D)"
+
+
+class CompositeBonds(Primitive):
+    """Linear Combination of several bond distances"""
+
+    def __init__(self, bonds: List[Tuple[int, int]], coeffs: List[float]):
+        """
+        Linear combination of a list of bonds and the corresponding
+        coefficients given as a list of real numbers
+
+        Args:
+            bonds: A list of tuples (i, j) representing bonds
+            coeffs: A list of floating point coefficients in order
+        """
+        super().__init__()
+        assert len(bonds) == len(coeffs), "Number of bonds != coefficients"
+        assert all(isinstance(bond, tuple) for bond in bonds)
+        assert all(len(bond) == 2 for bond in bonds)
+        assert all(
+            isinstance(bond[0], int) and isinstance(bond[1], int)
+            for bond in bonds
+        )
+        assert len(set(bonds)) == len(bonds)
+
+        self._bonds = list(bonds)
+        self._coeffs = [float(c) for c in coeffs]
+
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ) -> VectorHyperDual:
+        """Linear combination of bonds"""
+        all_idxs = list(itertools.chain(*self._bonds))
+        unique_idxs = list(set(all_idxs))
+        _x = x.ravel()
+        atom_vecs = _get_3d_vecs_from_atom_idxs(
+            *unique_idxs, x=_x, deriv_order=deriv_order
+        )
+
+        bonds_combined = None
+        for idx, (i, j) in enumerate(self._bonds):
+            atom_i = atom_vecs[unique_idxs.index(i)]
+            atom_j = atom_vecs[unique_idxs.index(j)]
+            if bonds_combined is None:
+                bonds_combined = self._coeffs[0] * (atom_i - atom_j).norm()
+            else:
+                bonds_combined += self._coeffs[idx] * (atom_i - atom_j).norm()
+
+        assert isinstance(bonds_combined, VectorHyperDual)
+        return bonds_combined
+
+    def __eq__(self, other):
+        """Equality of two linear combination of bonds"""
+        return (
+            isinstance(other, self.__class__)
+            and set(zip(self._bonds)) == set(zip(other._bonds))
+            and np.allclose(self._coeffs, other._coeffs)
+        )  # fmt: skip
+
+    def __repr__(self):
+        return f"CombinationOfBonds(n={len(self._bonds)})"
+
+
+class ConstrainedCompositeBonds(ConstrainedPrimitive, CompositeBonds):
+    """Constrained linear combindation of bonds"""
+
+    def __init__(
+        self, bonds: List[Tuple[int, int]], coeffs: List[float], value: float
+    ):
+        """
+        Linear combination of a list of bonds and the corresponding
+        coefficients given as a list of real numbers
+
+        Args:
+            bonds: A list of tuples (i, j) representing bonds
+            coeffs: A list of floating point coefficients in order
+            value: The target value for this coordinate
+        """
+        CompositeBonds.__init__(self, bonds=bonds, coeffs=coeffs)
+        self._r0 = value
+
+    @property
+    def _value(self) -> float:
+        return self._r0
+
+    def __repr__(self):
+        return f"ConstrainedCombinationOfBonds(n={len(self._bonds)})"
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 a28aa85..0d30a76 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -6,104 +6,234 @@ Notation follows:
 [2] J. Baker, J. Comput. Chem., 13, 240 Ž1992
 """
 import numpy as np
-from typing import Union
+from typing import Union, Optional, List, TYPE_CHECKING
 
 from autode.log import logger
 from autode.values import GradientRMS, Distance
 from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints
 from autode.opt.coordinates.internals import AnyPIC
 from autode.opt.optimisers.rfo import RFOptimiser
+from autode.exceptions import OptimiserStepError
 from autode.opt.optimisers.hessian_update import (
     BFGSDampedUpdate,
     BFGSSR1Update,
 )
 
+if TYPE_CHECKING:
+    from autode.opt.coordinates.primitives import Primitive
+
+# max and min bounds for the trust radius
+MAX_TRUST = 0.2
+MIN_TRUST = 0.01
+
 
 class CRFOptimiser(RFOptimiser):
+    """Constrained optimisation in delocalised internal coordinates"""
+
     def __init__(
-        self, init_alpha: Union[Distance, float] = 0.05, *args, **kwargs
+        self,
+        init_trust: float = 0.1,
+        *args,
+        extra_prims: Optional[List["Primitive"]] = None,
+        trust_update: bool = True,
+        max_move: Union[Distance, float] = Distance(0.12, "ang"),
+        **kwargs,
     ):
         """
         Constrained rational function optimisation
 
         -----------------------------------------------------------------------
         Arguments:
-            init_alpha: Maximum step size, assumed Angstrom if units
-                        not given
+            init_alpha: Initial value of the trust radius
+
+        Keyword Args:
+            extra_prims: A list of aditional coordinates (or constraints) to
+                        add to the DIC optimisation space (optional)
+            max_move: The maximum distance an atom can move in Cartesian
+                    coordinates in a step (assumed units of Å if not given)
+            trust_update: Whether to update the trust radius
 
         See Also:
-            :py:meth:`RFOOptimiser <RFOOptimiser.__init__>`
+            :py:meth:`RFOptimiser <RFOptimiser.__init__>`
         """
         super().__init__(*args, **kwargs)
 
-        self.alpha = Distance(init_alpha, units="ang")
-        assert self.alpha > 0
+        if not (MIN_TRUST < init_trust < MAX_TRUST):
+            init_trust = min(max(init_trust, MIN_TRUST), MAX_TRUST)
+            logger.warning(f"Setting trust radius to {init_trust:.3f}")
+
+        self.alpha = float(init_trust)
+        self._trust_update = bool(trust_update)
+        self._maxmove = Distance(max_move, units="ang")
+        assert self._maxmove > 0
+        self._extra_prims = [] if extra_prims is None else list(extra_prims)
+
         self._hessian_update_types = [BFGSDampedUpdate, BFGSSR1Update]
 
+    def _log_constrained_opt_progress(self):
+        """Log information about the constraints"""
+        n, m = len(self._coords), self._coords.n_constraints
+        s = self._coords.n_satisfied_constraints
+        logger.info(f"Optimising {n} coordinates and {m} lagrange multipliers")
+
+        idxs = self._coords.active_indexes
+        logger.info(
+            f"Satisfied {s} constraints. Active space"
+            f" is {len(idxs)} dimensional"
+        )
+        d2l_ev = np.linalg.eigvalsh(self._coords.h[:, idxs][idxs, :])
+        logger.info(
+            f"Hessian in active space has {sum(k < 0 for k in d2l_ev)} "
+            f"negative eigenvalue(s). Should have {m-s}"
+        )
+        return None
+
     def _step(self) -> None:
         """Partitioned rational function step"""
         assert self._coords is not None, "Must have coords to take a step"
-        assert self._coords.g is not None, "Must have a gradient"
+
         if self.iteration != 0:
             self._coords.update_h_from_old_h(
                 self._history.penultimate, self._hessian_update_types
             )
         assert self._coords.h is not None
 
-        n, m = len(self._coords), self._coords.n_constraints
-        logger.info(f"Optimising {n} coordinates and {m} lagrange multipliers")
+        self._update_trust_radius()
+        self._log_constrained_opt_progress()
 
-        idxs = self._coords.active_indexes
-        n_satisfied_constraints = (n + m - len(idxs)) // 2
-        logger.info(
-            f"Satisfied {n_satisfied_constraints} constraints. "
-            f"Active space is {len(idxs)} dimensional"
-        )
+        # get RFO step
+        delta_s = self._get_rfo_step()
 
-        b, u = np.linalg.eigh(self._coords.h[:, idxs][idxs, :])
-        logger.info(
-            f"∇^2L has {sum(lmda < 0 for lmda in b)} negative "
-            f"eigenvalue(s). Should have {m - n_satisfied_constraints}"
-        )
+        # scale back to trust radius only on non-constraint modes
+        n = len(self._coords)
+        delta_s_q = delta_s[:n]
+        if np.linalg.norm(delta_s_q) > self.alpha:
+            logger.info("Scaling RFO step to trust radius")
+            delta_s = delta_s * self.alpha / np.linalg.norm(delta_s_q)
 
-        logger.info("Calculating transformed gradient vector")
-        f = u.T.dot(self._coords.g[idxs])
+        logger.info("Taking an RFO step")
+        self._take_step_within_max_move(delta_s)
+        return None
 
-        lambda_p = self._lambda_p_from_eigvals_and_gradient(b, f)
-        logger.info(f"Calculated λ_p=+{lambda_p:.8f}")
+    def _get_rfo_step(self):
+        """
+        Calculate the unscaled RFO step, for the correct set of
+        coordinates
 
-        lambda_n = self._lambda_n_from_eigvals_and_gradient(b, f)
-        logger.info(f"Calculated λ_n={lambda_n:.8f}")
+        Returns:
+            (np.ndarray): The RFO step
+        """
+        n, m = len(self._coords), self._coords.n_constraints
+        idxs = self._coords.active_indexes
 
-        # Create the step along the n active DICs and m lagrangian multipliers
-        delta_s_active = np.zeros(shape=(len(idxs),))
+        # only molec. Hessian should be +ve definite
+        lmda = self._coords.rfo_shift
+        hess = self._coords.h - lmda * np.eye(n + m)
+        # no shift on constraints
+        for i in range(m):
+            hess[-m + i, -m + i] = 0.0
+
+        logger.info(f"Calculated RFO λ = {lmda:.4f}")
+        # RFO step in active space
+        hess = hess[:, idxs][idxs, :]
+        grad = self._coords.g[idxs]
+        self._check_shifted_hessian_has_correct_struct(hess)
+        full_step = np.zeros(shape=(n + m))
+        rfo_step = -np.matmul(np.linalg.inv(hess), grad)
+        full_step[idxs] = rfo_step
+
+        return full_step
+
+    def _check_shifted_hessian_has_correct_struct(self, arr) -> None:
+        """
+        Check that the shifted Hessian from RFO or QA has correct
+        eigenvalue structure
 
-        o = m - n_satisfied_constraints
-        logger.info(f"Maximising {o} modes")
+        Args:
+            arr (np.ndarray): Shifted hessian to check
 
-        for i in range(o):
-            delta_s_active -= f[i] * u[:, i] / (b[i] - lambda_p)
+        Raises:
+            (OptimiserStepError): if Hessian does not have correct structure
+        """
+        assert self._coords is not None
+        m = self._coords.n_constraints
+        o = m - self._coords.n_satisfied_constraints
+        ev = np.linalg.eigvalsh(arr)
+        n_negative = sum(k < 0 for k in ev)
+        if not o == n_negative:
+            raise OptimiserStepError(
+                f"Failed to obtain step, shifted Hessian should have {o}"
+                f" negative eigenvalue(s), but has {n_negative}"
+            )
+        return None
 
-        for j in range(n - n_satisfied_constraints):
-            delta_s_active -= f[o + j] * u[:, o + j] / (b[o + j] - lambda_n)
+    def _take_step_within_max_move(self, delta_s: np.ndarray):
+        """
+        Take the step by converting internal coordinates to Cartesian
+        coordinates, and scaling back if the maximum movement of an
+        atom exceeds max_move
 
-        # Set all the non-active components of the step to zero
-        delta_s = np.zeros(shape=(n + m,))
-        delta_s[idxs] = delta_s_active
+        Arguments:
+            delta_s (np.ndarray): The step in internal coordinates
+        """
+        assert self._coords is not None
+
+        self._coords.allow_unconverged_back_transform = True
+        new_coords = self._coords + delta_s
+        cart_delta = new_coords.to("cart") - self._coords.to("cart")
+        cart_displ = np.linalg.norm(cart_delta.reshape((-1, 3)), axis=1)
+        max_displ = np.abs(cart_displ).max()
+
+        self._coords.allow_unconverged_back_transform = False
+        if max_displ > self._maxmove:
+            logger.info(
+                f"Calculated step too large: max. displacement = "
+                f"{max_displ:.3f} Å, scaling down"
+            )
+            # Note because the transformation is not linear this will not
+            # generate a step exactly max(∆x) ≡ α, but is empirically close
+            factor = self._maxmove / max_displ
+            self._coords = self._coords + (factor * delta_s)
+
+        else:
+            self._coords = new_coords
 
-        c = self._take_step_within_trust_radius(delta_s[:n])
-        self._coords.update_lagrange_multipliers(c * delta_s[n:])
         return None
 
-    @property
-    def _g_norm(self) -> GradientRMS:
-        """Calculate the norm of the gradient in the active subspace"""
+    def _update_trust_radius(self):
+        """Updates the trust radius before a geometry step"""
+        assert self._coords is not None, "Must have coordinates!"
 
-        if self._coords is None or self._coords.g is None:
-            return super()._g_norm
+        if self.iteration == 0:
+            return None
 
-        gradient = self._coords.g[self._coords.active_indexes]
-        return GradientRMS(np.sqrt(np.mean(np.square(gradient))))
+        if self._trust_update is False:
+            return None
+
+        coords_l = self._history.penultimate
+        pred_delta_e = coords_l.pred_quad_delta_e(self._coords)
+        trust_ratio = self.last_energy_change / float(pred_delta_e)
+        last_step_size = np.linalg.norm(
+            np.array(coords_l) - np.array(self._coords)
+        )
+
+        if trust_ratio < 0.25:
+            self.alpha = max(0.7 * self.alpha, MIN_TRUST)
+        elif 0.25 < trust_ratio < 0.75:
+            pass
+        elif 0.75 < trust_ratio < 1.25:
+            # increase if step was actually near trust radius
+            if abs(last_step_size - self.alpha) / self.alpha < 0.05:
+                self.alpha = min(1.3 * self.alpha, MAX_TRUST)
+        elif 1.25 < trust_ratio < 1.75:
+            pass
+        elif trust_ratio > 1.75:
+            self.alpha = max(0.7 * self.alpha, MIN_TRUST)
+
+        logger.info(
+            f"Ratio of actual/predicted dE = {trust_ratio:.3f},"
+            f" Current trust radius = {self.alpha:.3f}"
+        )
 
     def _initialise_run(self) -> None:
         """Initialise the optimisation"""
@@ -127,45 +257,10 @@ class CRFOptimiser(RFOptimiser):
 
         cartesian_coords = CartesianCoordinates(self._species.coordinates)
         primitives = AnyPIC.from_species(self._species)
+        for ic in self._extra_prims:
+            primitives.add(ic)
 
         self._coords = DICWithConstraints.from_cartesian(
             x=cartesian_coords, primitives=primitives
         )
-        self._coords.zero_lagrangian_multipliers()
         return None
-
-    def _lambda_p_from_eigvals_and_gradient(
-        self, b: np.ndarray, f: np.ndarray
-    ) -> float:
-        """
-        Calculate the positive eigenvalue of the augmented hessian from
-        the eigenvalues of the full lagrangian Hessian matrix (b) and the
-        gradient in the Hessian eigenbasis
-        """
-        assert self._coords is not None
-
-        m = self._coords.n_constraints - self._coords.n_satisfied_constraints
-
-        aug_h = np.zeros(shape=(m + 1, m + 1))
-        aug_h[:m, :m] = np.diag(b[:m])
-        aug_h[:-1, -1] = aug_h[-1, :-1] = f[:m]
-
-        eigenvalues = np.linalg.eigvalsh(aug_h)
-        return eigenvalues[-1]
-
-    def _lambda_n_from_eigvals_and_gradient(
-        self, b: np.ndarray, f: np.ndarray
-    ) -> float:
-        """Like lambda_p but for the negative component"""
-        assert self._coords is not None
-
-        n_satisfied = self._coords.n_satisfied_constraints
-        m = self._coords.n_constraints - n_satisfied
-        n = len(self._coords) - n_satisfied
-
-        aug_h = np.zeros(shape=(n + 1, n + 1))
-        aug_h[:n, :n] = np.diag(b[m:])
-        aug_h[-1, :n] = aug_h[:n, -1] = f[m:]
-
-        eigenvalues = np.linalg.eigvalsh(aug_h)
-        return eigenvalues[0]
diff --git a/autode/opt/optimisers/prfo.py b/autode/opt/optimisers/prfo.py
index 4b8ae41..6a5a741 100644
--- a/autode/opt/optimisers/prfo.py
+++ b/autode/opt/optimisers/prfo.py
@@ -1,13 +1,12 @@
 """Partitioned rational function optimisation"""
 import numpy as np
-from typing import Union
+from typing import Union, Optional
 
 from autode.log import logger
 from autode.values import Distance
 from autode.opt.optimisers.crfo import CRFOptimiser
 from autode.opt.optimisers.hessian_update import BofillUpdate
 from autode.opt.coordinates.cartesian import CartesianCoordinates
-from autode.exceptions import CalculationException
 
 
 class PRFOptimiser(CRFOptimiser):
@@ -15,6 +14,7 @@ class PRFOptimiser(CRFOptimiser):
         self,
         init_alpha: Union[Distance, float] = 0.05,
         recalc_hessian_every: int = 10,
+        imag_mode_idx: int = 0,
         *args,
         **kwargs,
     ):
@@ -24,12 +24,12 @@ class PRFOptimiser(CRFOptimiser):
         all others to locate a transition state (TS)
 
         -----------------------------------------------------------------------
+
         Arguments:
             init_alpha: Maximum step size (default Å if unit not given)
 
-            imag_mode_idx: Index of the imaginary mode to follow. Default = 0,
-                           the most imaginary mode (i.e. most negative
-                           eigenvalue)
+            imag_mode_idx: Index of the imaginary mode to follow. Default
+                        is 0th mode i.e. the most negative mode
 
         See Also:
             :py:meth:`RFOOptimiser <RFOOptimiser.__init__>`
@@ -39,6 +39,8 @@ class PRFOptimiser(CRFOptimiser):
         self.alpha = Distance(init_alpha, units="ang")
         assert self.alpha > 0
         self.recalc_hessian_every = int(recalc_hessian_every)
+        self._mode_idx = imag_mode_idx
+        self._last_eigvec: Optional[np.ndarray] = None  # store last mode
         self._hessian_update_types = [BofillUpdate]
 
     def _step(self) -> None:
@@ -52,33 +54,75 @@ class PRFOptimiser(CRFOptimiser):
                 self._history.penultimate, self._hessian_update_types
             )
 
-        assert self._coords.h is not None  # _update_hessian must set .h
-        idxs = list(range(len(self._coords)))
+        assert self._coords.h is not None  # must set .h
 
-        b, u = np.linalg.eigh(self._coords.h[:, idxs][idxs, :])
+        b, u = np.linalg.eigh(self._coords.h)
+        f = u.T.dot(self._coords.g)
         n_negative_eigenvalues = sum(lmda < 0 for lmda in b)
         logger.info(
             f"∇^2E has {n_negative_eigenvalues} negative "
             f"eigenvalue(s). Should have 1"
         )
 
-        if n_negative_eigenvalues < 1:
-            raise CalculationException("Lost imaginary (TS) mode")
+        imag_idx = self._get_imag_mode_idx(u)
+        logger.info(f"Following mode {imag_idx} uphill")
+
+        b_max = b[imag_idx]
+        u_max = u[:, imag_idx]
+        f_max = f[imag_idx]
+
+        b_min = np.delete(b, imag_idx)
+        u_min = np.delete(u, imag_idx, axis=1)
+        f_min = np.delete(f, imag_idx)
 
-        f = u.T.dot(self._coords.g[idxs])
-        lambda_p = self._lambda_p_from_eigvals_and_gradient(b, f)
-        lambda_n = self._lambda_n_from_eigvals_and_gradient(b, f)
-        logger.info(f"Calculated λ_p=+{lambda_p:.8f}, λ_n={lambda_n:.8f}")
+        n = len(b)
+        delta_s = np.zeros(shape=(n,))
+        # downhill step
+        aug_h_min = np.zeros(shape=(n, n))
+        aug_h_min[: n - 1, : n - 1] = np.diag(b_min)
+        aug_h_min[:-1, -1] = aug_h_min[-1, :-1] = f_min
+        lambda_n = np.linalg.eigvalsh(aug_h_min)[0]
+        logger.info(f"Calculated λ_n = {lambda_n:.6f}")
 
-        delta_s = np.zeros(shape=(len(idxs),))
-        delta_s -= f[0] * u[:, 0] / (b[0] - lambda_p)
+        for i in range(n - 1):
+            delta_s -= f_min[i] * u_min[:, i] / (b_min[i] - lambda_n)
 
-        for j in range(1, len(idxs)):
-            delta_s -= f[j] * u[:, j] / (b[j] - lambda_n)
+        # uphill step
+        aug_h_max = np.zeros(shape=(2, 2))
+        aug_h_max[:1, :1] = b_max
+        aug_h_max[:-1, -1] = aug_h_max[-1, :-1] = f_max
+        lambda_p = np.linalg.eigvalsh(aug_h_max)[-1]
+        logger.info(f"Calculated λ_p = {lambda_p:.6f}")
 
-        _ = self._take_step_within_trust_radius(delta_s)
+        delta_s -= f_max * u_max / (b_max - lambda_p)
+
+        self._last_eigvec = u[:, imag_idx].flatten()
+        self._take_step_within_trust_radius(delta_s)
         return None
 
+    def _get_imag_mode_idx(self, u: np.ndarray):
+        """
+        Find the imaginary mode to follow upwards in the current step.
+
+        Args:
+            u (np.ndarray): The Hessian eigenvectors
+
+        Returns:
+            (int): Integer
+        """
+        if self.iteration == 0:
+            return self._mode_idx
+
+        overlaps = []
+        for i in range(u.shape[1]):
+            overlaps.append(
+                np.abs(np.dot(u[:, i].flatten(), self._last_eigvec))
+            )
+
+        mode_idx = np.argmax(overlaps)
+        logger.info(f"Overlap with previous TS mode: {overlaps[mode_idx]:.3f}")
+        return mode_idx
+
     def _initialise_run(self) -> None:
         """
         Initialise running a partitioned rational function optimisation by
@@ -86,7 +130,6 @@ class PRFOptimiser(CRFOptimiser):
         """
         assert self._species is not None, "Must have a species to init"
 
-        # self._build_internal_coordinates()
         self._coords = CartesianCoordinates(self._species.coordinates).to(
             "dic"
         )
diff --git a/autode/opt/optimisers/qa.py b/autode/opt/optimisers/qa.py
new file mode 100644
index 0000000..c1994b9
--- /dev/null
+++ b/autode/opt/optimisers/qa.py
@@ -0,0 +1,122 @@
+"""
+Constrained optimisation with quadratic trust radius model
+
+Also known as Quadratic Approximation (QA) or Trust-Radius Model (TRM)
+
+References:
+[1] P. Culot et al. Theor. Chim. Acta, 82, 1992, 189-205
+[2] T. Helgaker, Chem. Phys. Lett., 182(5), 1991, 503-510
+[3] J. T. Golab et al. Chem. Phys., 78, 1983, 175-199
+[4] R. Fletcher, Practical Methods of Optimization, Wiley, Chichester, 1981
+"""
+import numpy as np
+from scipy.optimize import root_scalar
+
+from autode.log import logger
+from autode.opt.optimisers.crfo import CRFOptimiser
+from autode.exceptions import OptimiserStepError
+
+
+class QAOptimiser(CRFOptimiser):
+    """Quadratic trust-radius optimiser in delocalised internal coordinates"""
+
+    def _step(self) -> None:
+        """Trust radius step"""
+        assert self._coords is not None, "Must have coords to take a step"
+
+        if self.iteration != 0:
+            self._coords.update_h_from_old_h(
+                self._history.penultimate, self._hessian_update_types
+            )
+        assert self._coords.h is not None
+
+        self._update_trust_radius()
+        self._log_constrained_opt_progress()
+
+        n = len(self._coords)
+
+        # Take RFO step if within trust radius
+        delta_s_rfo = self._get_rfo_step()
+        if np.linalg.norm(delta_s_rfo[:n]) <= self.alpha:
+            logger.info("Taking an RFO step")
+            self._take_step_within_max_move(delta_s_rfo)
+            return None
+
+        # otherwise use QA step within trust
+        try:
+            delta_s_qa = self._get_qa_step()
+            logger.info("Taking a QA step within trust radius")
+            self._take_step_within_max_move(delta_s_qa)
+            return None
+
+        # if QA fails, used scaled RFO step
+        except OptimiserStepError as exc:
+            logger.info(f"QA step failed: {str(exc)}, using scaled RFO step")
+            factor = self.alpha / np.linalg.norm(delta_s_rfo[:n])
+            self._take_step_within_max_move(delta_s_rfo * factor)
+            return None
+
+    def _get_qa_step(self):
+        """
+        Calculate the QA step within trust radius for the current
+        set of coordinates
+
+        Returns:
+            (np.ndarray): The trust radius step
+        """
+        n, m = len(self._coords), self._coords.n_constraints
+        idxs = self._coords.active_indexes
+
+        def shifted_newton_step(hess, grad, lmda, check=False):
+            """
+            Level-shifted Newton step (H-λI)^-1 . g
+            optional check of Hessian eigenvalue structure
+            """
+            hess = hess - lmda * np.eye(hess.shape[0])
+            # no shift on constraints
+            for i in range(m):
+                hess[-m + i, -m + i] = 0.0
+            full_step = np.zeros_like(grad)
+            hess = hess[:, idxs][idxs, :]
+            grad = grad[idxs]
+            if check:
+                self._check_shifted_hessian_has_correct_struct(hess)
+            qa_step = -np.matmul(np.linalg.inv(hess), grad)
+            full_step[idxs] = qa_step
+            return full_step
+
+        def qa_step_error(lmda):
+            """Error in step size"""
+            ds = shifted_newton_step(self._coords.h, self._coords.g, lmda)
+            ds_atoms = ds[:n]
+            return np.linalg.norm(ds_atoms) - self.alpha
+
+        # if molar Hessian +ve definite & step within trust use simple qN
+        min_b = self._coords.min_eigval
+        if min_b > 0 and qa_step_error(0.0) <= 0.0:
+            return shifted_newton_step(
+                self._coords.h, self._coords.g, 0.0, True
+            )
+
+        # Find λ in range (-inf, b)
+        for k in range(500):
+            right_bound = min_b - 0.5**k
+            if qa_step_error(right_bound) > 0:
+                break
+        assert qa_step_error(right_bound) > 0
+
+        for k in range(-6, 10):
+            left_bound = right_bound - 2**k
+            if qa_step_error(left_bound) < 0:
+                break
+        if not qa_step_error(left_bound) < 0:
+            raise OptimiserStepError("Unable to find bounds for root search")
+
+        res = root_scalar(f=qa_step_error, bracket=[left_bound, right_bound])
+        if (not res.converged) or (res.root >= min_b):
+            raise OptimiserStepError("QA root search failed")
+
+        logger.info(f"Calculated QA λ = {res.root:.4f}")
+        return shifted_newton_step(
+            self._coords.h, self._coords.g, res.root, True
+        )
diff --git a/autode/opt/optimisers/rfo.py b/autode/opt/optimisers/rfo.py
index dbb0533..25bae00 100644
--- a/autode/opt/optimisers/rfo.py
+++ b/autode/opt/optimisers/rfo.py
@@ -35,7 +35,7 @@ class RFOptimiser(NDOptimiser):
         """
         super().__init__(*args, **kwargs)
 
-        self.alpha = Distance(init_alpha, units="ang")
+        self.alpha = float(Distance(init_alpha, units="ang"))
         assert self.alpha > 0
         self._hessian_update_types = [BFGSPDUpdate, NullUpdate]
 
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/autode/values.py b/autode/values.py
index ea771c6..81d2cde 100644
--- a/autode/values.py
+++ b/autode/values.py
@@ -619,9 +619,10 @@ class ValueArray(ABC, np.ndarray):
             other is not None
             and hasattr(other, "shape")
             and other.shape == self.shape
-            and np.allclose(self, other)
+            and np.allclose(
+                np.asarray(self), np.asarray(other), atol=1e-64, rtol=1e-64
+            )
         )
-
         return eq
 
     def __ne__(self, other):
diff --git a/autode/wrappers/XTB.py b/autode/wrappers/XTB.py
index 338b05e..a685f02 100644
--- a/autode/wrappers/XTB.py
+++ b/autode/wrappers/XTB.py
@@ -109,6 +109,12 @@ class XTB(autode.wrappers.methods.ExternalMethodOEG):
     def print_xcontrol_file(self, calc: "CalculationExecutor", molecule):
         """Print an XTB input file with constraints and point charges"""
 
+        # do not print constraints if not optimisation
+        if not isinstance(calc.input.keywords, OptKeywords) and (
+            calc.input.point_charges is None
+        ):
+            return
+
         xcontrol_filename = f"xcontrol_{calc.name}"
         with open(xcontrol_filename, "w") as xcontrol_file:
             self.print_distance_constraints(xcontrol_file, molecule)
diff --git a/doc/changelog.rst b/doc/changelog.rst
index d70b8d3..fffecaf 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -1,6 +1,25 @@
 Changelog
 =========
 
+1.4.4
+------
+-------
+
+Functionality improvements
+**************************
+- Improved constrained optimisation (:code:`CRFOptimiser`) and handling of multiple constraints
+- Adds compatability with numpy v2.0
+- Improved implementation of the RFO-TRM (:code:`QAOptimiser`) optimiser that can handle constraints
+- Added static internal back-transform and damping for faster and easier DIC to Cartesian coordinate transformation
+
+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
 ------
 -------
diff --git a/setup.py b/setup.py
index 3854df2..8e145d2 100644
--- a/setup.py
+++ b/setup.py
@@ -32,7 +32,7 @@ extensions = [
 
 setup(
     name="autode",
-    version="1.4.3",
+    version="1.4.4",
     python_requires=">3.7",
     packages=[
         "autode",


Test patch:
diff --git a/tests/test_atoms.py b/tests/test_atoms.py
index 4c115c0..0832656 100644
--- a/tests/test_atoms.py
+++ b/tests/test_atoms.py
@@ -41,10 +41,12 @@ def test_atoms():
 
     h_atoms = Atoms([Atom("H"), Atom("H", x=1.0)])
     assert isinstance(h_atoms.com, Coordinate)
-    assert np.allclose(h_atoms.com, np.array([0.5, 0.0, 0.0]))
+    assert np.allclose(np.asarray(h_atoms.com), np.array([0.5, 0.0, 0.0]))
     assert not h_atoms.contain_metals
 
-    assert h_atoms.vector(0, 1) == np.array([1.0, 0.0, 0.0])
+    v = h_atoms.vector(0, 1)
+    assert isinstance(v, np.ndarray)
+    assert np.allclose(v, np.array([1.0, 0.0, 0.0]))
 
     # Moment of inertia
     assert np.sum(np.diag(h_atoms.moi)) > 0.0
diff --git a/tests/test_bracket/test_dhs.py b/tests/test_bracket/test_dhs.py
index e7e39b4..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)
@@ -157,7 +146,9 @@ def test_dhs_image_pair_ts_guess(caplog):
 
     assert np.allclose(peak.coordinates.flatten(), imgpair.left_coords)
     assert np.isclose(peak.energy, -0.144)
-    assert np.allclose(peak.gradient.flatten(), imgpair.left_coords.g)
+    assert np.allclose(
+        np.asarray(peak.gradient.flatten()), np.asarray(imgpair.left_coords.g)
+    )
 
 
 @requires_working_xtb_install
@@ -255,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)
@@ -316,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,
     )
@@ -344,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 eb8e210..9e4ee3a 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -22,7 +22,7 @@ from autode.opt.coordinates.internals import (
     AnyPIC,
 )
 from autode.opt.coordinates.cartesian import CartesianCoordinates
-from autode.opt.coordinates.dic import DIC
+from autode.opt.coordinates.dic import DIC, DICWithConstraints
 from autode.opt.coordinates.primitives import (
     PrimitiveInverseDistance,
     PrimitiveDistance,
@@ -30,9 +30,12 @@ from autode.opt.coordinates.primitives import (
     PrimitiveBondAngle,
     ConstrainedPrimitiveBondAngle,
     PrimitiveDihedralAngle,
+    PrimitiveImproperDihedral,
     PrimitiveLinearAngle,
     PrimitiveDummyLinearAngle,
     LinearBendType,
+    CompositeBonds,
+    ConstrainedCompositeBonds,
 )
 
 
@@ -111,10 +114,10 @@ def test_cartesian_coordinates_hessian_update():
     # Simple coordinates with 2 atoms in 3 D
     coords = CartesianCoordinates(np.arange(0, 6).reshape((2, 3)))
 
-    with pytest.raises(ValueError):
+    with pytest.raises(AssertionError):
         coords.update_h_from_cart_h(arr=np.array([]))
 
-    with pytest.raises(ValueError):
+    with pytest.raises(AssertionError):
         coords.update_h_from_cart_h(arr=np.array([1.0]))
 
     # Hessian needs to be 6x6
@@ -165,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]])
 
@@ -199,6 +221,29 @@ def test_basic_dic_properties():
         _ = x.to("unknown coordinates")
 
 
+def test_dic_constraints():
+    mol = water_mol()
+    mol.constraints.distance = {(0, 1): 1.5}
+
+    pic = AnyPIC.from_species(mol)
+    x = CartesianCoordinates(mol.coordinates)
+    q = DICWithConstraints.from_cartesian(x, pic)
+    assert q.n_constraints == 1
+    assert q.g is None and q.h is None
+
+    # with constraints grad or hessian cannot be set directly
+    with pytest.raises(RuntimeError):
+        q.g = np.arange(3)
+    with pytest.raises(RuntimeError):
+        q.h = np.arange(3)
+
+    q.update_g_from_cart_g(np.random.rand(9))
+    q.update_h_from_cart_h(np.random.rand(9, 9))
+    # one extra dimension from Lagrange multiplier
+    assert q.g.shape == (4,)
+    assert q.h.shape == (4, 4)
+
+
 def test_invalid_pic_construction():
     # Cannot construct some primitives e.g. PrimitiveInverseDistances from non Primitive
     # internal coordinates
@@ -236,6 +281,13 @@ def test_cart_to_dic():
     x -= 0.1
 
 
+def test_cartesian_all_indexes_active():
+    arr = np.arange(6)
+    x = CartesianCoordinates(arr)
+    assert x.active_indexes == list(range(6))
+    assert x.inactive_indexes == list()
+
+
 def test_simple_dic_to_cart():
     arr = np.array([[0.0, 0.0, 0.0], [2.0, 0.0, 0.0]])
 
@@ -638,6 +690,16 @@ def test_dihedral_equality():
     )
 
 
+def test_composite_bonds_equality():
+    a = CompositeBonds(bonds=[(1, 2), (2, 3)], coeffs=[0.5, 1.2])
+    b = CompositeBonds(bonds=[(1, 2), (2, 3)], coeffs=[0.5, 1.2])
+    c = CompositeBonds(bonds=[(0, 5), (2, 4)], coeffs=[0.5, 1.2])
+    d = CompositeBonds(bonds=[(1, 2), (2, 3)], coeffs=[0.1, 1.2])
+    assert a == b
+    assert a != c  # different bonds
+    assert a != d  # different coefficient
+
+
 def test_linear_angle():
     acetylene = Molecule(
         atoms=[
@@ -654,8 +716,8 @@ def test_linear_angle():
     assert angle._vec_r is not None
     old_r_vec = angle._vec_r
     # the dummy atom should not change after the first call
-    _ = angle(x)
-    _ = angle(x)
+    _ = angle(x + 0.05)
+    _ = angle(x - 0.07)
     assert angle._vec_r is old_r_vec
 
     axis_vec = np.array(np.array(angle._vec_r._data) - x.reshape(-1, 3)[1])
@@ -684,6 +746,9 @@ def test_primitives_consistent_with_mol_values():
     assert np.isclose(ang(coords), h2o2.angle(0, 2, 1), rtol=1e-8)
     dihedral = PrimitiveDihedralAngle(2, 0, 1, 3)
     assert np.isclose(dihedral(coords), h2o2.dihedral(2, 0, 1, 3), rtol=1e-8)
+    ic = CompositeBonds([(0, 1), (0, 2)], [0.3, 0.7])
+    mol_val = 0.3 * h2o2.distance(0, 1) + 0.7 * h2o2.distance(0, 2)
+    assert np.isclose(mol_val, ic(coords))
 
 
 # fmt: off
@@ -713,6 +778,7 @@ extra_mols = [
         ]
     ),
     feco5_mol(), # for testing linear angles
+    h2o2_mol(),
 ]
 
 test_mols = [
@@ -726,6 +792,7 @@ test_prims = [
     PrimitiveDihedralAngle(2, 0, 1, 3),
     PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND),
     PrimitiveLinearAngle(2, 3, 4, 8, LinearBendType.BEND),
+    CompositeBonds([(0, 1), (0, 2)], [1, 1]),
 ]
 # fmt: on
 
@@ -783,6 +850,9 @@ def test_repr():
         PrimitiveLinearAngle(0, 1, 2, 3, LinearBendType.BEND),
         PrimitiveLinearAngle(0, 1, 2, 3, LinearBendType.COMPLEMENT),
         PrimitiveDummyLinearAngle(0, 1, 2, LinearBendType.BEND),
+        PrimitiveImproperDihedral(0, 1, 2, 3),
+        CompositeBonds(bonds=[(1, 2), (2, 3)], coeffs=[1, 1]),
+        ConstrainedCompositeBonds([(1, 2), (2, 3)], [1, 1], 0.2),
     ]
 
     for p in prims:
@@ -844,11 +914,11 @@ def test_pic_generation_linear_angle_ref():
     assert not any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic)
     assert PrimitiveLinearAngle(4, 3, 2, 8, LinearBendType.BEND) in pic
     # for C-Fe-C, only one out-of-plane dihedral should be present
-    assert PrimitiveDihedralAngle(3, 5, 2, 1) in pic
+    assert PrimitiveImproperDihedral(3, 5, 2, 1) in pic
     assert sum(isinstance(ic, PrimitiveDihedralAngle) for ic in pic) == 1
     # check degrees of freedom = 3N - 6
-    _ = pic(m.coordinates.flatten())
-    assert np.linalg.matrix_rank(pic.B) == 3 * m.n_atoms - 6
+    x = m.coordinates.flatten()
+    assert np.linalg.matrix_rank(pic.get_B(x)) == 3 * m.n_atoms - 6
 
 
 def test_pic_generation_linear_angle_dummy():
@@ -864,8 +934,8 @@ def test_pic_generation_linear_angle_dummy():
     assert any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic)
 
     # degrees of freedom = 3N - 5 for linear molecules
-    _ = pic(mol.coordinates.flatten())
-    assert np.linalg.matrix_rank(pic.B) == 3 * mol.n_atoms - 5
+    x = mol.coordinates.flatten()
+    assert np.linalg.matrix_rank(pic.get_B(x)) == 3 * mol.n_atoms - 5
 
 
 @work_in_tmp_dir()
@@ -907,8 +977,8 @@ def test_pic_generation_disjoint_graph():
     assert PrimitiveDistance(2, 3) not in pic
     assert PrimitiveBondAngle(1, 2, 3) not in pic
     # check degrees of freedom = 3N - 6
-    _ = pic(mol.coordinates.flatten())
-    assert np.linalg.matrix_rank(pic.B) == 3 * mol.n_atoms - 6
+    x = mol.coordinates.flatten()
+    assert np.linalg.matrix_rank(pic.get_B(x)) == 3 * mol.n_atoms - 6
 
     # if the bond between 2, 3 is made into a constraint, it will generate angles
     mol.constraints.distance = {(2, 3): mol.distance(2, 3)}
@@ -928,8 +998,8 @@ def test_pic_generation_chain_dihedrals():
     assert PrimitiveDihedralAngle(7, 4, 3, 6) in pic
 
     # check that the 3N-6 degrees of freedom are maintained
-    _ = pic(cumulene.coordinates.flatten())
-    assert np.linalg.matrix_rank(pic.B) == 3 * cumulene.n_atoms - 6
+    x = cumulene.coordinates.flatten()
+    assert np.linalg.matrix_rank(pic.get_B(x)) == 3 * cumulene.n_atoms - 6
 
 
 def test_pic_generation_square_planar():
@@ -947,5 +1017,5 @@ def test_pic_generation_square_planar():
     # for sq planar, out-of-plane dihedrals are needed to have
     # all degrees of freedom
     pic = AnyPIC.from_species(ptcl4)
-    _ = pic(ptcl4.coordinates.flatten())
-    assert np.linalg.matrix_rank(pic.B) == 3 * ptcl4.n_atoms - 6
+    x = ptcl4.coordinates.flatten()
+    assert np.linalg.matrix_rank(pic.get_B(x)) == 3 * ptcl4.n_atoms - 6
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index f040278..d7a26ee 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -9,14 +9,17 @@ from autode.methods import XTB
 from autode.opt.coordinates.internals import PIC
 from autode.opt.optimisers.crfo import CRFOptimiser
 from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints
-from autode.opt.coordinates.primitives import PrimitiveDihedralAngle
+from autode.opt.coordinates.primitives import (
+    PrimitiveDihedralAngle,
+    ConstrainedCompositeBonds,
+)
 from autode.utils import work_in_tmp_dir
 from .molecules import h2o2_mol, acetylene_mol, feco5_mol, cumulene_mol
 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()
 
@@ -46,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
@@ -62,7 +65,6 @@ def test_coordinate_setup():
 
     # Initial lagrangian multiplier is close to zero, which is the last
     # component in the optimisation space
-    opt._coords.zero_lagrangian_multipliers()
     assert np.isclose(opt._coords._lambda[0], 0.0)
 
 
@@ -70,14 +72,14 @@ def crfo_water_coords():
     return crfo_coords(molecule=water_molecule())
 
 
-def test_setting_invalid_lagrange_multipliers():
+def test_adding_invalid_step():
     s = crfo_water_coords()
-    # the vector of lagrange multipliers must be of length one to be
-    # appropriate for a water molecule with a single constrained distance
-    invalid_multipliers = np.ones(shape=(4,))
+    # the added step must be the length of the coordinates plus
+    # one lagragne multiplier
+    invalid_step = np.ones(shape=(5,))
 
-    with pytest.raises(ValueError):
-        s.update_lagrange_multipliers(invalid_multipliers)
+    with pytest.raises(AssertionError):
+        s += invalid_step
 
 
 def test_simple_gradient_update():
@@ -109,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
@@ -121,7 +123,7 @@ def test_primitive_projection_discard():
     assert len(s) == 3
 
     # Shift on the first couple of DIC but nothing on the final one
-    s += np.array([0.03, -0.07, 0.0])
+    s += np.array([0.03, -0.07, 0.0, 0.0])
 
     def r(_x):
         return np.linalg.norm(_x[:3] - _x[3:6])
@@ -130,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([]))
 
@@ -161,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)
@@ -182,7 +179,7 @@ def test_step_c2h3():
     coords = crfo_coords(m)
 
     # Should be able to add an arbitrary vector to the coordinates
-    coords += np.random.uniform(-0.1, 0.1, size=coords.shape)
+    coords += np.random.uniform(-0.1, 0.1, size=coords.raw.shape)
 
 
 def test_baker1997_example():
@@ -257,13 +254,13 @@ 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()
 
     # Moving the angle should preserve the distances
-    s = opt._coords + np.array([0.1, 0.0, 0.0])
+    s = opt._coords + np.array([0.1, 0.0, 0.0, 0.0, 0.0])  # append multiplier
     x = s.to("cart").reshape((3, 3))
     for i, j in ((0, 1), (0, 2)):
         assert np.isclose(
@@ -277,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
 
@@ -349,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")
@@ -363,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")
 
@@ -373,8 +370,60 @@ 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=15, 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")
     assert abs(mol.dihedral(6, 3, 4, 8)) < val.Angle(95, "deg")
+
+
+@requires_working_xtb_install
+@work_in_tmp_dir()
+def test_composite_bond_constraint():
+    atoms = [
+        Atom("C", -2.6862, 1.0780, -0.1640),
+        Atom("C", -2.1836, -0.0798, -0.5820),
+        Atom("C", -0.5315, 2.2123, -0.0294),
+        Atom("H", -1.1416, -0.3055, -0.5232),
+        Atom("C", -1.8773, 2.2239, 0.0534),
+        Atom("H", 0.0331, 1.3028, -0.0181),
+        Atom("H", -0.4693, 2.6248, -1.9811),
+        Atom("C", -0.2061, 1.6165, -2.2424),
+        Atom("H", -2.8152, -0.9212, -0.7994),
+        Atom("H", -2.3965, 3.1639, 0.1692),
+        Atom("H", 0.8400, 1.4387, -2.4301),
+        Atom("C", -1.1022, 0.6635, -2.5738),
+        Atom("H", -3.7553, 1.2122, -0.0856),
+        Atom("H", 0.0433, 3.1121, 0.1002),
+        Atom("H", -2.1548, 0.8243, -2.5510),
+        Atom("C", -0.6762, -0.6020, -3.1972),
+        Atom("O", 0.4285, -0.9110, -3.5228),
+        Atom("H", -1.5372, -1.3068, -3.3795),
+    ]
+    # asymmetric Diels-Alder (butadiene + acrolein)
+    mol = Molecule(atoms=atoms)
+    constr = ConstrainedCompositeBonds(
+        bonds=[(1, 11), (2, 7)], coeffs=[1, 1], value=4.6
+    )
+    # 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, 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)
+    # difference should be much higher, as TS is asymmetric
+    assert mol.distance(1, 11) - mol.distance(2, 7) > 0.1
+
+
+def test_trust_radius_limits():
+    import autode.opt.optimisers.crfo
+
+    max_lim = autode.opt.optimisers.crfo.MAX_TRUST
+    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, 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 af216de..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:
@@ -341,11 +332,44 @@ def test_multiple_optimiser_saves_overrides_not_append():
     assert old_n_coords == n_coords
 
 
+@work_in_tmp_dir()
+def test_optimiser_plotting_sanity_checks(caplog):
+    mol = Molecule(smiles="N#N")
+    opt = CartesianSDOptimiser(maxiter=10, conv_tol="loose")
+    coords1 = CartesianCoordinates(mol.coordinates)
+    coords1.e = PotentialEnergy(0.1, "Ha")
+    coords1.update_g_from_cart_g(
+        np.array([0.01, 0.02, 0.05, 0.06, 0.03, 0.07])
+    )
+    opt._coords = coords1
+    opt._species = mol
+    assert opt.iteration == 0
+    assert not opt.converged
+    # plotting does not work if less than 2 points
+    with caplog.at_level("WARNING"):
+        opt.plot_optimisation(filename="test-plot.pdf")
+    assert not os.path.isfile("test-plot.pdf")
+    assert "Less than 2 points, cannot draw optimisation" in caplog.text
+
+    opt._coords = coords1.copy()
+    opt._coords.e = PotentialEnergy(0.0, "Ha")
+    assert not opt.converged
+    # either rms_g or energy plot has to be requested
+    with caplog.at_level("ERROR"):
+        opt.plot_optimisation("test-plot.pdf", False, False)
+    assert not os.path.isfile("test-plot.pdf")
+    assert "Must plot either energies or RMS gradients" in caplog.text
+    with caplog.at_level("WARNING"):
+        opt.plot_optimisation("test-plot.pdf", plot_energy=True)
+    assert os.path.isfile("test-plot.pdf")
+    assert "Optimisation is not converged, drawing a plot" in caplog.text
+
+
 @work_in_tmp_dir()
 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 d48815e..7a915a4 100644
--- a/tests/test_opt/test_prfo.py
+++ b/tests/test_opt/test_prfo.py
@@ -3,6 +3,7 @@ from autode.species.molecule import Molecule
 from autode.atoms import Atom
 from autode.methods import XTB
 from autode.opt.optimisers import PRFOptimiser
+from autode.config import Config
 from autode.utils import work_in_tmp_dir
 from ..testutils import requires_working_xtb_install
 
@@ -67,6 +68,42 @@ def test_diels_alder_ts_opt():
         print(xyz_file_string, file=file)
 
     mol = Molecule("init.xyz")
-    PRFOptimiser.optimise(mol, method=xtb, maxiter=50, init_alpha=0.05)
+    PRFOptimiser.optimise(mol, method=xtb, maxiter=25, init_alpha=0.05)
     assert has_single_imag_freq_at_xtb_level(mol)
-    # print(mol.imaginary_frequencies)  # should be ~600 cm-1
+    freq = mol.imaginary_frequencies[0]
+    assert np.isclose(freq, -600, atol=30)  # should be ~600 cm-1
+
+
+@requires_working_xtb_install
+@work_in_tmp_dir()
+def test_mode_following():
+    mol = Molecule(
+        name="sn2_ts",
+        charge=-1,
+        solvent_name="water",
+        atoms=[
+            Atom("F", -4.17085, 3.55524, 1.59944),
+            Atom("Cl", -0.75962, 3.53830, -0.72354),
+            Atom("C", -2.51988, 3.54681, 0.47836),
+            Atom("H", -3.15836, 3.99230, -0.27495),
+            Atom("H", -2.54985, 2.47411, 0.62732),
+            Atom("H", -2.10961, 4.17548, 1.25945),
+        ],
+    )
+    opt = PRFOptimiser(maxiter=10, conv_tol="normal", imag_mode_idx=0)
+    opt._species = mol
+    opt._method = xtb
+    opt._n_cores = Config.n_cores
+    opt._initialise_run()
+    # take a step
+    opt._step()
+    opt._update_gradient_and_energy()
+    opt._update_hessian()
+    # shift the Hessian modes by exchanging eigenvalues of first two modes
+    b, u = np.linalg.eigh(opt._coords.h)
+    b[0], b[1] = b[1], b[0]
+    new_h = np.linalg.multi_dot((u, np.diag(b), u.T)).real
+    b, u = np.linalg.eigh(new_h)
+    new_idx = opt._get_imag_mode_idx(u)
+    # the chosen index should be 1, based on overlap
+    assert new_idx == 1
diff --git a/tests/test_opt/test_qa.py b/tests/test_opt/test_qa.py
new file mode 100644
index 0000000..fa6770f
--- /dev/null
+++ b/tests/test_opt/test_qa.py
@@ -0,0 +1,99 @@
+import os
+import numpy as np
+import pytest
+from autode import Molecule, Atom, Config
+from autode.methods import XTB
+from autode.opt.optimisers.qa import QAOptimiser
+from autode.opt.coordinates import DICWithConstraints
+from autode.utils import work_in_tmp_dir
+from ..testutils import work_in_zipped_dir, requires_working_xtb_install
+
+
+here = os.path.dirname(os.path.abspath(__file__))
+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, conv_tol="loose", init_trust=0.1)
+    opt._species = mol
+    opt._build_internal_coordinates()
+    assert isinstance(opt._coords, DICWithConstraints)
+
+    grad = np.loadtxt("opt-test_grad.txt")
+    hess = np.loadtxt("opt-test_hess.txt")
+    opt._coords.update_g_from_cart_g(grad)
+    opt._coords.update_h_from_cart_h(hess)
+
+    opt._step()
+    step = np.array(opt._history.final) - np.array(opt._history.penultimate)
+    step_size = np.linalg.norm(step)
+    assert np.isclose(step_size, 0.1)
+
+
+@work_in_tmp_dir()
+@requires_working_xtb_install
+def test_trust_update():
+    init_trust = 0.05
+    water_atoms = [
+        Atom("O", -0.0011, 0.3631, -0.0000),
+        Atom("H", -0.8250, -0.1819, -0.0000),
+        Atom("H", 0.8261, -0.1812, 0.0000),
+    ]
+    water = Molecule(atoms=water_atoms)
+    opt = QAOptimiser(maxiter=10, conv_tol="loose", init_trust=init_trust)
+
+    opt._species = water.copy()
+    opt._method = XTB()
+    opt._n_cores = Config.n_cores
+    opt._initialise_run()
+    # store last grad
+    last_g = opt._coords.g.copy()
+    last_h = opt._coords.h.copy()
+
+    opt._step()
+    opt._update_gradient_and_energy()
+    last_step = np.array(opt._coords) - np.array(opt._history[-2])
+    pred_delta_e = float(np.dot(last_g, last_step))
+    pred_delta_e += 0.5 * np.linalg.multi_dot((last_step, last_h, last_step))
+    # pred_dE should be around -0.002544605 Ha (depends on xTB version)
+
+    def simulate_energy_change_ratio_update_trust(ratio):
+        opt.alpha = init_trust
+        opt._history.final.e = (
+            opt._history.penultimate.e + ratio * pred_delta_e
+        )
+        opt._update_trust_radius()
+
+    # should not update if trust update turned off
+    opt._trust_update = False
+    simulate_energy_change_ratio_update_trust(0.2)
+    assert np.isclose(opt.alpha, init_trust)
+
+    opt._trust_update = True
+    simulate_energy_change_ratio_update_trust(0.2)
+    assert np.isclose(opt.alpha, 0.7 * init_trust)
+
+    simulate_energy_change_ratio_update_trust(0.5)
+    assert np.isclose(opt.alpha, init_trust)
+
+    simulate_energy_change_ratio_update_trust(1.0)
+    assert (np.linalg.norm(last_step) - init_trust) / init_trust < 0.05
+    assert np.isclose(opt.alpha, 1.3 * init_trust)
+
+    simulate_energy_change_ratio_update_trust(1.3)
+    assert np.isclose(opt.alpha, init_trust)
+
+    simulate_energy_change_ratio_update_trust(1.8)
+    assert np.isclose(opt.alpha, 0.7 * init_trust)
+
+
+@work_in_tmp_dir()
+@requires_working_xtb_install
+def test_molecular_opt_qa():
+    mol = Molecule(smiles="CCO")
+    constr_distance = mol.distance(1, 3) + 0.1
+    mol.constraints.distance = {(1, 3): constr_distance}
+    QAOptimiser.optimise(mol, method=XTB(), maxiter=10)
+    assert np.isclose(mol.distance(1, 3), constr_distance, 1e-6)
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)
diff --git a/tests/test_values.py b/tests/test_values.py
index 3955c01..e2c48e4 100644
--- a/tests/test_values.py
+++ b/tests/test_values.py
@@ -49,8 +49,7 @@ def test_coordinate():
     assert "coord" in repr(coord).lower()
 
     assert coord is not None
-    # Equality defaults to np.allclose
-    assert coord == np.zeros(3)
+    assert np.allclose(coord, np.zeros(3))
 
 
 def test_coordinates():
@@ -114,8 +113,7 @@ def test_inplace_modification():
 def test_copy_conversion():
     x = Gradient([[1.0, 1.0, 1.0]], units="Ha / Å")
     y = x.to("eV / Å")
-
-    assert not np.allclose(x, y)
+    assert not np.allclose(np.asarray(x), np.asarray(y))
     assert np.allclose(x, np.ones(shape=(1, 3)))