DHS improvements (duartegroup#338)

* recursion error bugfix

* add steepest descent step

* add test

* add test

* add test

* add damped step

* add logging

* add tests

* pr suggestion

autode/bracket/dhs.py

@@ -148,8 +148,7 @@ 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.last_energy_change < self.etol
+            self.rms_tangent_grad < self.gtol and self._abs_delta_e < self.etol
         )
 
     @property
@@ -195,23 +194,65 @@ def _step(self) -> None:
         the pivot point. A line search is done if it is not the first
         iteration (and it has not been turned off), and then a
         quasi-Newton step with a Lagrangian constraint for the distance
-        is taken
+        is taken (falls back to steepest descent with projected gradient
+        if this fails).
         """
         assert self._coords is not None, "Must have set coordinates"
 
+        # if energy is rising, interpolate halfway between last step
+        if self.iteration >= 1 and (
+            self.last_energy_change > PotentialEnergy(5, "kcalmol")
+        ):
+            logger.warning("Energy rising, going back half a step")
+            half_interp = (self._coords + self._history.penultimate) / 2
+            self._coords = half_interp
+            return None
+
         if self.iteration >= 1 and self._do_line_search:
             coords, grad = self._line_search_on_sphere()
         else:
             coords, grad = self._coords, self._coords.g
 
-        step = self._get_lagrangian_step(coords, grad)
-
-        step_size = np.linalg.norm(step)
-        logger.info(f"Taking a quasi-Newton step: {step_size:.3f} Å")
+        try:
+            step = self._get_lagrangian_step(coords, grad)
+            logger.info(
+                f"Taking a quasi-Newton step: {np.linalg.norm(step):.3f} Å"
+            )
+        except OptimiserStepError:
+            step = self._get_sd_step(coords, grad)
+            logger.warning(
+                f"Failed to take quasi-Newton step, taking steepest "
+                f"descent step instead: {np.linalg.norm(step):.3f} Å"
+            )
 
         # the step is on the interpolated coordinates (if done)
         actual_step = (coords + step) - self._coords
         self._coords = self._coords + actual_step
+        return None
+
+    def _get_sd_step(self, coords, grad) -> np.ndarray:
+        """
+        Obtain a steepest descent step minimising the tangential
+        gradient. This step cannot perfectly maintain the same
+        distance from pivot point. The step size is at most half
+        of the trust radius.
+
+        Args:
+            coords: Previous coordinates
+            grad: Previous gradient
+
+        Returns:
+            (np.ndarray): Step in Cartesian coordinates
+        """
+        dist_vec = coords - self._pivot
+        dist_hat = dist_vec / np.linalg.norm(dist_vec)
+        perp_grad = grad - np.dot(grad, dist_hat) * dist_hat
+
+        sd_step = -perp_grad
+        if np.linalg.norm(sd_step) > self.alpha / 2:
+            sd_step *= (self.alpha / 2) / np.linalg.norm(sd_step)
+
+        return sd_step
 
     def _get_lagrangian_step(self, coords, grad) -> np.ndarray:
         """

tests/test_bracket/test_dhs.py

@@ -14,6 +14,7 @@
     TruncatedTaylor,
     DHSImagePair,
     ImageSide,
+    OptimiserStepError,
 )
 from autode.opt.coordinates import CartesianCoordinates
 from autode import Config
@@ -97,6 +98,53 @@ def test_distance_constrained_optimiser():
     assert np.isclose(new_distance, distance)
 
 
+@work_in_zipped_dir(datazip)
+def test_dist_constr_optimiser_sd_fallback():
+    coords1 = CartesianCoordinates(np.loadtxt("conopt_last.txt"))
+    coords1.update_g_from_cart_g(np.loadtxt("conopt_last_g.txt"))
+    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
+    opt = DistanceConstrainedOptimiser(
+        pivot_point=pivot,
+        maxiter=2,
+        init_trust=0.2,
+        gtol=1e-3,
+        etol=1e-4,
+    )
+    opt._target_dist = 2.6869833732268
+    opt._history.open("test_trj")
+    opt._history.add(coords1)
+    with pytest.raises(OptimiserStepError):
+        opt._get_lagrangian_step(coords1, coords1.g)
+    # however, steepest descent step should work
+    sd_step = opt._get_sd_step(coords1, coords1.g)
+    opt._step()
+    assert np.allclose(opt._coords, coords1 + sd_step)
+
+
+@work_in_tmp_dir()
+def test_dist_constr_optimiser_energy_rising():
+    coords1 = CartesianCoordinates(np.arange(6, dtype=float))
+    coords1.e = PotentialEnergy(0.01, "Ha")
+    step = np.random.random(6)
+    coords2 = coords1 + step
+    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,
+    )
+    opt._history.open("test_trj")
+    opt._history.add(coords1)
+    opt._history.add(coords2)
+    opt._step()
+    assert np.allclose(opt._coords, coords1 + step * 0.5)
+
+
 def test_dhs_image_pair():
     mol1 = Molecule(smiles="CCO")
     mol2 = mol1.new_species()