Repository: duartegroup/autodE
Pull request #342 (2024-05-24 19:46:48)

DHS variable step size
(1) Uses multiple step sizes depending on the distance thresholds
(2) Replaces only gradient based peak detection to cubic polynomial based peak detection (both energies and gradients)


***

## Checklist

* [ ] The changes include an associated explanation of how/why
* [ ] Test pass
* [ ] Documentation has been updated
* [ ] Changelog has been updated


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/duartegroup/autodE/pull/342?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) Report
Attention: Patch coverage is `98.15951%` with `12 lines` in your changes are missing coverage. Please review.
> Project coverage is 97.34%. Comparing base [(`034743c`)](https://app.codecov.io/gh/duartegroup/autodE/commit/034743cb1b679711fb3801ebbd97f02a4b3c1acf?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) to head [(`98ac75f`)](https://app.codecov.io/gh/duartegroup/autodE/commit/98ac75fd099db5a9deb06ad62b8eabb63d3a85a5?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup).

| [Files](https://app.codecov.io/gh/duartegroup/autodE/pull/342?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) | Patch % | Lines |
|---|---|---|
| [autode/opt/optimisers/utils.py](https://app.codecov.io/gh/duartegroup/autodE/pull/342?src=pr&el=tree&filepath=autode%2Fopt%2Foptimisers%2Futils.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL29wdC9vcHRpbWlzZXJzL3V0aWxzLnB5) | 88.13% | [7 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/342?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [autode/bracket/ieip.py](https://app.codecov.io/gh/duartegroup/autodE/pull/342?src=pr&el=tree&filepath=autode%2Fbracket%2Fieip.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL2JyYWNrZXQvaWVpcC5weQ==) | 91.42% | [3 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/342?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [autode/bracket/dhs.py](https://app.codecov.io/gh/duartegroup/autodE/pull/342?src=pr&el=tree&filepath=autode%2Fbracket%2Fdhs.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL2JyYWNrZXQvZGhzLnB5) | 96.07% | [2 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/342?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |

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


```diff
@@            Coverage Diff             @@
##           master     #342      +/-   ##
==========================================
- Coverage   97.44%   97.34%   -0.10%     
==========================================
  Files         209      201       -8     
  Lines       23747    23160     -587     
==========================================
- Hits        23141    22546     -595     
- Misses        606      614       +8     
```

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

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


</details>

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

shoubhikraj: @t-young31 Could you please check one last time before merge?

Patch:
diff --git a/autode/bracket/dhs.py b/autode/bracket/dhs.py
index 9503815..0a9c39b 100644
--- a/autode/bracket/dhs.py
+++ b/autode/bracket/dhs.py
@@ -11,7 +11,8 @@ from enum import Enum
 
 from autode.values import Distance, Angle, GradientRMS, PotentialEnergy
 from autode.bracket.imagepair import EuclideanImagePair
-from autode.opt.coordinates import OptCoordinates, CartesianCoordinates
+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
@@ -23,50 +24,6 @@ if TYPE_CHECKING:
     from autode.wrappers.methods import Method
 
 
-class TruncatedTaylor:
-    """The truncated taylor surface from current grad and hessian"""
-
-    def __init__(
-        self,
-        centre: Union[OptCoordinates, np.ndarray],
-        grad: np.ndarray,
-        hess: np.ndarray,
-    ):
-        """
-        Second-order Taylor expansion around a point
-
-        Args:
-            centre (OptCoordinates|np.ndarray): The coordinate point
-            grad (np.ndarray): Gradient at that point
-            hess (np.ndarray): Hessian at that point
-        """
-        self.centre = centre
-        if hasattr(centre, "e") and centre.e is not None:
-            self.e = centre.e
-        else:
-            # the energy can be relative and need not be absolute
-            self.e = 0.0
-        self.grad = grad
-        self.hess = hess
-        n_atoms = grad.shape[0]
-        assert hess.shape == (n_atoms, n_atoms)
-
-    def value(self, coords: np.ndarray) -> float:
-        """Energy (or relative energy if point did not have energy)"""
-        # E = E(0) + g^T . dx + 0.5 * dx^T. H. dx
-        dx = (coords - self.centre).flatten()
-        new_e = self.e + np.dot(self.grad, dx)
-        new_e += 0.5 * np.linalg.multi_dot((dx, self.hess, dx))
-        return new_e
-
-    def gradient(self, coords: np.ndarray) -> np.ndarray:
-        """Gradient at supplied coordinate"""
-        # g = g(0) + H . dx
-        dx = (coords - self.centre).flatten()
-        new_g = self.grad + np.matmul(self.hess, dx)
-        return new_g
-
-
 class DistanceConstrainedOptimiser(RFOptimiser):
     """
     Constrained optimisation of a molecule, with the Euclidean
@@ -500,8 +457,8 @@ class DHSImagePair(EuclideanImagePair):
         self, side: ImageSide, step_size: float
     ) -> np.ndarray:
         """
-        Obtain the DHS step on the specified side, with the specified step
-        size
+        Obtain the DHS extrapolation step on the specified side,
+        with the specified step size
 
         Args:
             side (ImageSide): left or right
@@ -531,33 +488,37 @@ class DHS(BaseBracketMethod):
         self,
         initial_species: "Species",
         final_species: "Species",
-        step_size: Union[Distance, float] = Distance(0.1, "ang"),
+        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"),
         **kwargs,
     ):
         """
-        Dewar-Healy-Stewart method to find transition states.
-
-        1) The step size is 0.1 which is quite conservative, so may want
-        to increase that if convergence is slow; 2) The distance tolerance
-        should not be lowered any more than 1.0 Angstrom as DHS is unstable
-        when the distance is low, and there is a tendency for one image to
-        jump over the barrier
+        Dewar-Healy-Stewart method to find transition states. The distance
+        tolerance convergence criteria should not be much lower than 0.5 Angstrom
+        as DHS is unstable when the distance is low, and there is a tendency for
+        one image to jumpy over the barrier.
 
         Args:
             initial_species: The "reactant" species
 
             final_species: The "product" species
 
-            step_size: The size of the DHS step taken along
-                        the linear path between reactant and
-                        product in Angstroms
+            large_step: The size of the DHS step when distance between the
+                        images is larger than switch_thresh (Angstrom)
+
+            small_step: The size of the DHS step when distance between the
+                        images is smaller than swtich_thresh (Angstrom)
+
+            switch_thresh: When distance between the two images is less than
+                        this cutoff, smaller DHS extrapolation steps are taken
 
         Keyword Args:
 
             maxiter: Maximum number of en/grad evaluations
 
             dist_tol: The distance tolerance at which DHS will
-                      stop, values less than 1.0 Angstrom are not
+                      stop, values less than 0.5 Angstrom are not
                       recommended.
 
             gtol: Gradient tolerance for the optimiser micro-iterations
@@ -577,18 +538,22 @@ class DHS(BaseBracketMethod):
         self._method: Optional[Method] = None
         self._n_cores: Optional[int] = None
 
-        self._step_size = Distance(abs(step_size), "ang")
-        if self._step_size > self.imgpair.dist:
+        self._large_step = Distance(abs(large_step), "ang")
+        self._small_step = Distance(abs(small_step), "ang")
+        self._sw_thresh = Distance(abs(switch_thresh), "ang")
+        assert self._small_step < self._large_step
+
+        self._step_size: Optional[Distance] = None
+        if self._large_step > self.imgpair.dist:
             logger.warning(
-                f"Step size ({self._step_size:.3f} Å) for {self._name}"
+                f"Step size ({self._large_step:.3f} Å) for {self._name}"
                 f" is larger than the starting Euclidean distance between"
                 f" images ({self.imgpair.dist:.3f} Å). This calculation"
                 f" will likely run into errors."
             )
 
-        # NOTE: In DHS the micro-iterations are done separately, in
-        # an optimiser, so to keep track of the actual number of
-        # en/grad calls, this local variable is used
+        # NOTE: In DHS the micro-iterations are done separately in
+        # an optimiser, so keep track with local variable
         self._current_microiters: int = 0
 
     def _initialise_run(self) -> None:
@@ -607,6 +572,12 @@ class DHS(BaseBracketMethod):
         assert self._method is not None, "Must have a set method"
         assert self.imgpair.left_coords.e and self.imgpair.right_coords.e
 
+        if self.imgpair.dist > self._sw_thresh:
+            self._step_size = self._large_step
+        else:
+            self._step_size = self._small_step
+        opt_trust = min(self._step_size, Distance(0.1, "ang"))
+
         if self.imgpair.left_coords.e < self.imgpair.right_coords.e:
             side = ImageSide.left
             pivot = self.imgpair.right_coords
@@ -628,6 +599,7 @@ class DHS(BaseBracketMethod):
             maxiter=curr_maxiter,
             gtol=self._gtol,
             etol=1.0e-3,  # seems like a reasonable etol
+            init_trust=opt_trust,
             pivot_point=pivot,
             old_coords_read_hess=old_coords,
         )
@@ -648,6 +620,7 @@ class DHS(BaseBracketMethod):
 
         # put results back into imagepair
         self.imgpair.put_coord_by_side(opt.final_coordinates, side)  # type: ignore
+        opt.clean_up()
         return None
 
     def _calculate(
@@ -699,6 +672,7 @@ class DHS(BaseBracketMethod):
         Returns:
             (CartesianCoordinates): New predicted coordinates for that side
         """
+        assert self._step_size is not None
         # take a DHS step of the size given
         dhs_step = self.imgpair.get_dhs_step_by_side(side, self._step_size)
 
@@ -707,7 +681,7 @@ class DHS(BaseBracketMethod):
 
         logger.info(
             f"DHS step on {side} image: taking a step of"
-            f" size {self._step_size:.4f}"
+            f" size {np.linalg.norm(dhs_step):.4f} Å"
         )
         return new_coord
 
@@ -754,6 +728,7 @@ class DHSGS(DHS):
             (CartesianCoordinates): New predicted coordinates for that side
         """
         assert self.imgpair is not None, "Must have an image pair"
+        assert self._step_size is not None
 
         dhs_step = self.imgpair.get_dhs_step_by_side(side, self._step_size)
         gs_step = self.imgpair.get_last_step_by_side(side)
@@ -764,7 +739,7 @@ class DHSGS(DHS):
             dhs_step = dhs_step / (1 - self._gs_mix)
         else:
             # rescale GS step as well so that one vector doesn't dominate
-            gs_step *= self._step_size / np.linalg.norm(gs_step)
+            gs_step *= np.linalg.norm(dhs_step) / np.linalg.norm(gs_step)
 
         old_coord = self.imgpair.get_coord_by_side(side)
         new_coord = (
diff --git a/autode/bracket/imagepair.py b/autode/bracket/imagepair.py
index 6a1ad00..8023edc 100644
--- a/autode/bracket/imagepair.py
+++ b/autode/bracket/imagepair.py
@@ -5,7 +5,7 @@ that require a pair of images
 import itertools
 import numpy as np
 from abc import ABC, abstractmethod
-from typing import Optional, Tuple, TYPE_CHECKING, List, Iterator
+from typing import Optional, Tuple, TYPE_CHECKING, Union, Iterator
 from enum import Enum
 
 from autode.values import Distance, PotentialEnergy, Gradient
@@ -14,6 +14,7 @@ from autode.methods import get_lmethod
 from autode.neb import CINEB
 from autode.opt.coordinates import CartesianCoordinates
 from autode.opt.optimisers.hessian_update import BofillUpdate
+from autode.opt.optimisers.utils import Polynomial2PointFit
 from autode.opt.optimisers.base import OptimiserHistory, print_geometries_from
 from autode.plotting import plot_bracket_method_energy_profile
 from autode.utils import work_in_tmp_dir, ProcessPool
@@ -450,56 +451,27 @@ class EuclideanImagePair(BaseImagePair, ABC):
     @property
     def has_jumped_over_barrier(self) -> bool:
         """
-        A quick test of whether the images are still separated by a barrier
-        is to check whether the gradient vectors are pointing outwards
-        compared to the linear path connecting the two images. In case
-        there are multiple barriers in the way, a distance threshold is
-        also used to guess if it is likely that one image has jumped over.
-
-        This is a slightly modified version of the method proposed in ref:
-        Y. Liu, H. Qui, M. Lei, J. Chem. Theory. Comput., 2023
-        https://doi.org/10.1021/acs.jctc.3c00151
+        A quick test of whether the images are still separated by a barrier,
+        implemented via fitting a cubic polynomial along the linear path
+        connecting the two images and checking for a peak. This is only an
+        approximation.
         """
+        assert self.left_coords is not None and self.right_coords is not None
+        assert self.left_coords.e and self.right_coords.e
         assert (
             self.left_coords.g is not None and self.right_coords.g is not None
         )
+        cubic_poly = Polynomial2PointFit.cubic_fit(
+            self.left_coords, self.right_coords
+        )
 
-        def cos_angle(vec1, vec2) -> float:
-            """Returns the cos(theta) between two vectors (1D arrays)"""
-            dot = float(np.dot(vec1, vec2)) / (
-                np.linalg.norm(vec1) * np.linalg.norm(vec2)
-            )
-            return dot
-
-        # NOTE: The angle between the force vector on right image
-        # and the distance vector must be more than 90 degrees. Similarly
-        # for left image it must be less than 90 degrees. This would mean
-        # that the parallel component of the forces on each image are
-        # pointing away from each other. (force is negative gradient).
-
-        left_cos_theta = cos_angle(-self.left_coords.g, self.dist_vec)
-        right_cos_theta = cos_angle(-self.right_coords.g, self.dist_vec)
-
-        assert -1.0 < left_cos_theta < 1.0
-        assert -1.0 < right_cos_theta < 1.0
-
-        # cos(theta) < 0.0 means angle > 90 degrees and vice versa
-        if right_cos_theta < 0.0 < left_cos_theta:
+        # NOTE: Interpolation seems reasonable upto ~1.2 Angstrom. If distance
+        # is larger, detecting peak is impossible without calculating energies
+        # so we assume there is a barrier between the images
+        if self.dist > Distance(1.2, "ang"):
             return False
-
-        # However, if there are multiple barriers in the path (i.e. multi-step
-        # reaction), it will identify as having jumped over, even if it didn't.
-        # The distance between the images would be high if there are multiple
-        # barriers. A threshold of 1 angstrom is used as it seems the risk of
-        # jumping over is high (if there is only one barrier) below this.
-        # (according to Kilmes et al., J. Phys.: Condens. Matter, 22 2010, 074203)
-        # This is of course, somewhat arbitrary, and will not work if really
-        # large steps are taken OR two barriers are very close in distance
-
-        if self.dist <= Distance(1.0, "ang"):
-            return True
         else:
-            return False
+            return cubic_poly.get_extremum(0.0, 1.0, get_max=True) is None
 
     def run_cineb_from_end_points(self) -> None:
         """
diff --git a/autode/opt/optimisers/utils.py b/autode/opt/optimisers/utils.py
new file mode 100644
index 0000000..936d7a2
--- /dev/null
+++ b/autode/opt/optimisers/utils.py
@@ -0,0 +1,159 @@
+"""
+Various operations for coordinates and optimisers
+"""
+import numpy as np
+from numpy.polynomial import Polynomial
+from typing import Union, TYPE_CHECKING
+
+if TYPE_CHECKING:
+    from autode.opt.coordinates.base import OptCoordinates
+
+
+class TruncatedTaylor:
+    """The truncated taylor surface from current grad and hessian"""
+
+    def __init__(
+        self,
+        centre: Union["OptCoordinates", np.ndarray],
+        grad: np.ndarray,
+        hess: np.ndarray,
+    ):
+        """
+        Second-order Taylor expansion around a point
+
+        Args:
+            centre (OptCoordinates|np.ndarray): The coordinate point
+            grad (np.ndarray): Gradient at that point
+            hess (np.ndarray): Hessian at that point
+        """
+        self.centre = centre
+        if hasattr(centre, "e") and centre.e is not None:
+            self.e = centre.e
+        else:
+            # the energy can be relative and need not be absolute
+            self.e = 0.0
+        self.grad = grad
+        self.hess = hess
+        n_atoms = grad.shape[0]
+        assert hess.shape == (n_atoms, n_atoms)
+
+    def value(self, coords: np.ndarray) -> float:
+        """Energy (or relative energy if point did not have energy)"""
+        # E = E(0) + g^T . dx + 0.5 * dx^T. H. dx
+        dx = (coords - self.centre).flatten()
+        new_e = self.e + np.dot(self.grad, dx)
+        new_e += 0.5 * np.linalg.multi_dot((dx, self.hess, dx))
+        return new_e
+
+    def gradient(self, coords: np.ndarray) -> np.ndarray:
+        """Gradient at supplied coordinate"""
+        # g = g(0) + H . dx
+        dx = (coords - self.centre).flatten()
+        new_g = self.grad + np.matmul(self.hess, dx)
+        return new_g
+
+
+def _get_energies_proj_gradients(
+    coords0: "OptCoordinates", coords1: "OptCoordinates"
+):
+    """
+    Get energies and projected gradients from two set of
+    coordinates
+    """
+    assert coords0.e and coords1.e
+    assert coords0.g is not None and coords1.g is not None
+    dist_vec = coords1.raw - coords0.raw
+    e0 = float(coords0.e)
+    g0 = np.dot(coords0.g, dist_vec)
+    e1 = float(coords1.e)
+    g1 = np.dot(coords1.g, dist_vec)
+    return e0, e1, g0, g1
+
+
+class Polynomial2PointFit(Polynomial):
+    """
+    1D polynomial along the line connecting two coordinates
+    """
+
+    @classmethod
+    def cubic_fit(cls, coords0: "OptCoordinates", coords1: "OptCoordinates"):
+        """
+        Obtain a cubic polynomial from the energies and gradients
+        at two sets of coordinates: f(x) = d + cx + bx**2 + ax**3
+
+        Args:
+            coords0 (OptCoordinates):
+            coords1 (OptCoordinates):
+
+        Returns:
+            (Polynomial2PointFit): The fitted polynomial (normalised)
+        """
+        e0, e1, g0, g1 = _get_energies_proj_gradients(coords0, coords1)
+        # f(0) = d; f(1) = a + b + c + d
+        d = e0
+        # f'(0) = c => a + b = f(1) - c - d
+        c = g0
+        a_b = e1 - c - d
+        # f'(1) = 3a + 2b + c => 3a + 2b = f'(1) - c
+        a3_2b = g1 - c
+        a = a3_2b - 2 * a_b
+        b = a_b - a
+        return cls([d, c, b, a])
+
+    def get_extremum(
+        self, l_bound: float = 0.0, u_bound: float = 1.0, get_max=False
+    ) -> Union[float, None]:
+        """
+        Obtain the maximum/minimum of a polynomial f(x), within
+        two bounds. If there are multiple, return the highest/lowest
+        respectively.
+
+        Args:
+            l_bound (float):
+            u_bound (float):
+            get_max (bool): Maximum or minimum requested
+
+        Returns:
+            (float|None): The x value at min or max f(x), None
+                        if not found or not within bounds
+        """
+        # points with derivative 0 are critical points
+        crit_points = self.deriv().roots()
+
+        if l_bound > u_bound:
+            u_bound, l_bound = l_bound, u_bound
+        crit_points = crit_points[crit_points < u_bound]
+        crit_points = crit_points[crit_points > l_bound]
+
+        if len(crit_points) == 0:
+            return None
+
+        maxima = []
+        minima = []
+        for point in crit_points:
+            for i in range(2, 6):
+                ith_deriv = self.deriv(i)(point)
+                # if zero, move up an order of derivative
+                if -1.0e-14 < ith_deriv < 1.0e-14:
+                    continue
+                # derivative > 0 and i is even => max
+                elif ith_deriv < 0 and i % 2 == 0:
+                    maxima.append(point)
+                # derivative > 0 and i is even => min
+                elif ith_deriv > 0 and i % 2 == 0:
+                    minima.append(point)
+                # otherwise inflection point
+                else:
+                    break
+
+        if get_max:
+            if len(maxima) == 0:
+                return None
+            max_vals = [self(x) for x in maxima]
+            return maxima[np.argmax(max_vals)]
+
+        else:
+            if len(minima) == 0:
+                return None
+            min_vals = [self(x) for x in minima]
+            return minima[np.argmin(min_vals)]
diff --git a/doc/changelog.rst b/doc/changelog.rst
index 6e7c2e4..d70b8d3 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -7,7 +7,8 @@ Changelog
 
 Functionality improvements
 **************************
-- ...
+- DHS and DHS-GS can now switch between two step sizes
+- Peak detection in bracket methods now uses cubic polynomial fit with energies and gradients instead of only projecting gradients
 
 Bug Fixes
 *********


Test patch:
diff --git a/tests/test_bracket/test_dhs.py b/tests/test_bracket/test_dhs.py
index bcd3c3d..e7e39b4 100644
--- a/tests/test_bracket/test_dhs.py
+++ b/tests/test_bracket/test_dhs.py
@@ -11,7 +11,6 @@ from autode.bracket.dhs import (
     DHS,
     DHSGS,
     DistanceConstrainedOptimiser,
-    TruncatedTaylor,
     DHSImagePair,
     ImageSide,
     OptimiserStepError,
@@ -24,33 +23,6 @@ here = os.path.dirname(os.path.abspath(__file__))
 datazip = os.path.join(here, "data", "geometries.zip")
 
 
-@requires_working_xtb_install
-@work_in_tmp_dir()
-def test_truncated_taylor_surface():
-    mol = Molecule(smiles="CCO")
-    mol.calc_hessian(method=XTB())
-    coords = CartesianCoordinates(mol.coordinates)
-    coords.update_g_from_cart_g(mol.gradient)
-    coords.update_h_from_cart_h(mol.hessian)
-    coords.make_hessian_positive_definite()
-
-    # for positive definite hessian, minimum of taylor surface would
-    # be a simple Newton step
-    minim = coords - (np.linalg.inv(coords.h) @ coords.g)
-
-    # minimizing surface should give the same result
-    surface = TruncatedTaylor(coords, coords.g, coords.h)
-    res = minimize(
-        method="CG",
-        fun=surface.value,
-        x0=np.array(coords),
-        jac=surface.gradient,
-    )
-
-    assert res.success
-    assert np.allclose(res.x, minim, rtol=1e-4)
-
-
 @requires_working_xtb_install
 @work_in_zipped_dir(datazip)
 def test_distance_constrained_optimiser():
@@ -199,7 +171,8 @@ def test_dhs_single_step():
         initial_species=reactant,
         final_species=product,
         maxiter=200,
-        step_size=step_size,
+        large_step=0.2,
+        switch_thresh=1.5,
         dist_tol=1.0,
     )
 
@@ -213,13 +186,14 @@ def test_dhs_single_step():
     old_dist = imgpair.dist
     assert imgpair.left_coords.e > imgpair.right_coords.e
 
+    assert dhs.imgpair.dist > 1.5
     # take a single step
     dhs._step()
     # step should be on lower energy image
     assert len(imgpair._left_history) == 1
     assert len(imgpair._right_history) == 2
     new_dist = imgpair.dist
-    # image should move exactly by step_size
+    # image should move exactly by large step
     assert np.isclose(old_dist - new_dist, step_size)
 
 
@@ -235,7 +209,8 @@ def test_dhs_gs_single_step(caplog):
         initial_species=reactant,
         final_species=product,
         maxiter=200,
-        step_size=step_size,
+        large_step=step_size,
+        switch_thresh=1.5,
         dist_tol=1.0,
         gs_mix=0.5,
     )
@@ -243,6 +218,7 @@ def test_dhs_gs_single_step(caplog):
     dhs_gs.imgpair.set_method_and_n_cores(method=XTB(), n_cores=1)
     dhs_gs._method = XTB()
     dhs_gs._initialise_run()
+    assert dhs_gs.imgpair.dist > 1.5
     # take one step
     with caplog.at_level("INFO"):
         dhs_gs._step()
@@ -278,8 +254,8 @@ def test_dhs_diels_alder():
     dhs = DHS(
         initial_species=reactant,
         final_species=product,
+        small_step=0.2,
         maxiter=200,
-        step_size=0.2,
         dist_tol=set_dist_tol,
         gtol=1.0e-3,
     )
@@ -341,7 +317,8 @@ def test_dhs_jumping_over_barrier(caplog):
         initial_species=reactant,
         final_species=product,
         maxiter=200,
-        step_size=0.6,
+        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,
         barrier_check=True,
@@ -368,7 +345,8 @@ def test_dhs_stops_if_microiter_exceeded(caplog):
         initial_species=reactant,
         final_species=product,
         maxiter=10,
-        step_size=0.2,
+        large_step=0.2,
+        small_step=0.1,
         dist_tol=1.0,
         gtol=5.0e-4,
         barrier_check=True,
diff --git a/tests/test_opt/test_opt_utils.py b/tests/test_opt/test_opt_utils.py
new file mode 100644
index 0000000..c0507ef
--- /dev/null
+++ b/tests/test_opt/test_opt_utils.py
@@ -0,0 +1,84 @@
+from autode import Molecule
+from autode.utils import work_in_tmp_dir
+from autode.opt.coordinates import CartesianCoordinates
+from autode.bracket.imagepair import _calculate_engrad_for_species
+from autode.methods import XTB
+from autode.opt.optimisers.utils import TruncatedTaylor, Polynomial2PointFit
+from scipy.optimize import minimize
+import numpy as np
+from numpy.polynomial import Polynomial
+from ..testutils import requires_working_xtb_install
+
+
+@requires_working_xtb_install
+@work_in_tmp_dir()
+def test_truncated_taylor_surface():
+    mol = Molecule(smiles="CCO")
+    mol.calc_hessian(method=XTB())
+    coords = CartesianCoordinates(mol.coordinates)
+    coords.update_g_from_cart_g(mol.gradient)
+    coords.update_h_from_cart_h(mol.hessian)
+    coords.make_hessian_positive_definite()
+
+    # for positive definite hessian, minimum of taylor surface would
+    # be a simple Newton step
+    minim = coords - (np.linalg.inv(coords.h) @ coords.g)
+
+    # minimizing surface should give the same result
+    surface = TruncatedTaylor(coords, coords.g, coords.h)
+    res = minimize(
+        method="CG",
+        fun=surface.value,
+        x0=np.array(coords),
+        jac=surface.gradient,
+    )
+
+    assert res.success
+    assert np.allclose(res.x, minim, rtol=1e-4)
+
+
+@requires_working_xtb_install
+@work_in_tmp_dir()
+def test_cubic_interpolation():
+    mol = Molecule(smiles="CCO")
+
+    en1, grad1 = _calculate_engrad_for_species(mol, XTB(), 1)
+    coords1 = CartesianCoordinates(mol.coordinates)
+    coords1.e = en1
+    coords1.update_g_from_cart_g(grad1)
+    # take a steepest descent step
+    step = -coords1.g * 0.1 / np.linalg.norm(coords1.g)
+    coords2 = coords1 + step
+    mol.coordinates = coords2
+    en2, grad2 = _calculate_engrad_for_species(mol, XTB(), 1)
+    coords2.e = en2
+    coords2.update_g_from_cart_g(grad2)
+    # cubic interp
+    cubic_poly = Polynomial2PointFit.cubic_fit(coords1, coords2)
+    # energy at a point between 0 and 1
+    coords3 = coords1 + 0.25 * step
+    mol.coordinates = coords3
+    en3, _ = _calculate_engrad_for_species(mol, XTB(), 1)
+    interp_en = cubic_poly(0.25)
+    assert np.isclose(interp_en, en3, atol=1e-3)
+
+
+def test_polynomial_max_min():
+    # 1 - x + 2x^2 + x^3
+    poly = Polynomial([1, -1, 2, 1])
+    # known minimum = 0.215
+    x_min = Polynomial2PointFit.get_extremum(poly, 0, 1)
+    assert np.isclose(x_min, 0.215, atol=1e-3)
+    # known maximum = -1.549
+    x_max = Polynomial2PointFit.get_extremum(poly, -2, 0, get_max=True)
+    assert np.isclose(x_max, -1.549, atol=1e-3)
+    # order of bounds should not matter
+    x_max = Polynomial2PointFit.get_extremum(poly, 0, -2, get_max=True)
+    assert np.isclose(x_max, -1.549, atol=1e-3)
+    # should not count inflection points
+    poly = Polynomial([0, 0, 0, 1])
+    x_min = Polynomial2PointFit.get_extremum(poly, -np.inf, np.inf)
+    x_max = Polynomial2PointFit.get_extremum(
+        poly, -np.inf, np.inf, get_max=True
+    )
+    assert x_min is None and x_max is None