Add trust-radius quasi-Newton optimiser (duartegroup#262)

* stub commit

* add qa step

* fix qa step to use cartesian step size

* improvements to coordinates

* damping updates

* TRIM step updates

* hessian update tests

* TRIM tests

* fix qa step to use cartesian step size

* minor updates

* unfinished updates

* minor updates

* updates to QA

* optimiser plot

* optimiser plotting added

* stub for flowchart update

* flowchart update

* trust update

* damping optimiser

* hessian update tests

* DIC flag to disallow unconverged IBT

* improvements to coordinates

* damping updates

* TRIM step updates

* hessian update tests

* TRIM tests

* tests for trimoptimiser

* tests for plotting

* unconverged ibt is allowed by default

* add to init

* reject high energy increase

* print ratio as well

* BFGS-SR1 scheme added

* pr updates

* pr updates duartegroup#2

* rename

* bfgs-sr1-update tests

* cartesian optimiser

* reimplement damping

* minor fixes

* pr updates, test updates

* fixes for tests, change name of RFO trust argument

* revert RFO and change step size control

* test updates

* xyz trajectory writing

* polynomial interpolation: commit

* pr suggestions

* test update

* changes to TRM

* use hessian eigenbasis

* fix 0 eigenvalue error

autode/exceptions.py

@@ -128,3 +128,11 @@ class CouldNotPlotSmoothProfile(AutodeException):
 
 class TemplateLoadingFailed(AutodeException):
     """A template file was not in the correct format"""
+
+
+class OptimiserStepError(AutodeException):
+    """Unable to calculate a valid Optimiser step"""
+
+
+class CoordinateTransformFailed(AutodeException):
+    """Internal coordinate to Cartesian transform failed"""

autode/opt/coordinates/base.py

@@ -39,15 +39,22 @@ def __new__(
         arr.B = None  # Wilson B matrix
         arr.B_T_inv = None  # Generalised inverse of B
         arr.U = np.eye(len(arr))  # Transform matrix
+        arr.allow_unconverged_back_transform = True  # for internal coords
 
         return arr
 
     def __array_finalize__(self, obj: "OptCoordinates") -> None:
         """See https://numpy.org/doc/stable/user/basics.subclassing.html"""
 
-        for attr in ("units", "_e", "_g", "_h", "_h_inv", "U", "B", "B_T_inv"):
+        # fmt: off
+        for attr in (
+            "units",    "_e",   "_g",   "_h",
+            "_h_inv",   "U",    "B",    "B_T_inv",
+            "allow_unconverged_back_transform",
+        ):
             self.__dict__[attr] = getattr(obj, attr, None)
 
+        # fmt: on
         return None
 
     @property

autode/opt/coordinates/dic.py

@@ -24,6 +24,7 @@
     InverseDistances,
     InternalCoordinates,
 )
+from autode.exceptions import CoordinateTransformFailed
 
 _max_back_transform_iterations = 20
 
@@ -239,8 +240,10 @@ def iadd(
                     f"Final RMS(s) = {rms_s:.8f}"
                 )
                 x_k = x_1
-                if self._allow_unconverged_back_transform:
-                    break
+                if not self.allow_unconverged_back_transform:
+                    raise CoordinateTransformFailed(
+                        "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)
@@ -251,10 +254,6 @@ def iadd(
 
         return self
 
-    @property
-    def _allow_unconverged_back_transform(self) -> bool:
-        return True
-
     @property
     def active_indexes(self) -> List[int]:
         """A list of indexes for the active modes in this coordinate set"""

autode/opt/optimisers/__init__.py

@@ -3,6 +3,7 @@
 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.trm import HybridTRMOptimiser
 from autode.opt.optimisers.steepest_descent import (
     CartesianSDOptimiser,
     DIC_SD_Optimiser,
@@ -16,4 +17,5 @@
     "CRFOptimiser",
     "CartesianSDOptimiser",
     "DIC_SD_Optimiser",
+    "HybridTRMOptimiser",
 ]

autode/opt/optimisers/base.py

@@ -2,6 +2,7 @@
 import numpy as np
 
 from abc import ABC, abstractmethod
+from collections import UserList
 from typing import Union, Optional, Callable, Any
 from autode.log import logger
 from autode.utils import NumericStringDict
@@ -10,6 +11,7 @@
 from autode.opt.coordinates.base import OptCoordinates
 from autode.opt.optimisers.hessian_update import NullUpdate
 from autode.exceptions import CalculationException
+from autode.plotting import plot_optimiser_profile
 
 
 class BaseOptimiser(ABC):
@@ -66,7 +68,7 @@ def __init__(
         self._maxiter = int(maxiter)
         self._n_cores: int = Config.n_cores
 
-        self._history = _OptimiserHistory()
+        self._history = OptimiserHistory()
 
         self._coords = coords
         self._species: Optional["autode.species.Species"] = None
@@ -769,8 +771,82 @@ def _best_hessian_updater(self) -> "HessianUpdater":
             "suitable update strategies"
         )
 
+    def plot_optimisation(
+        self,
+        filename: Optional[str] = None,
+        plot_energy: bool = True,
+        plot_rms_grad: bool = True,
+    ) -> None:
+        """
+        Draw the plot of the energies and/or rms_gradients of
+        the optimisation so far
+
+        ----------------------------------------------------------------------
+        Args:
+            filename (str): Name of the file to plot
+            plot_energy (bool): Whether to plot energy
+            plot_rms_grad (bool): Whether to plot RMS of gradient
+        """
+        if self.iteration < 1:
+            logger.warning("Less than 2 points, cannot draw optimisation plot")
+            return None
+
+        if not self.converged:
+            logger.warning(
+                "Optimisation is not converged, drawing a plot "
+                "of optimiser profile until current iteration"
+            )
+
+        filename = (
+            f"{self._species.name}_opt_plot.pdf"
+            if filename is None
+            else str(filename)
+        )
+
+        plot_optimiser_profile(
+            self._history,
+            filename=filename,
+            plot_energy=plot_energy,
+            plot_rms_grad=plot_rms_grad,
+        )
+        return None
+
+    def print_geometries(self, filename: Optional[str] = None) -> None:
+        """
+        Write the trajectory of the optimiser in .xyz format
+
+        Args:
+            filename: Name of the trajectory file (xyz)
+        """
+        assert self._species is not None
+        if self.iteration < 1:
+            logger.warning(
+                "Optimiser did no steps, not saving .xyz trajectory"
+            )
+            return None
+
+        filename = (
+            f"{self._species.name}_opt.trj.xyz"
+            if filename is None
+            else str(filename)
+        )
+
+        if os.path.isfile(filename):
+            logger.warning(f"File {filename} exists, overwriting...")
+            os.remove(filename)
+
+        tmp_spc = self._species.new_species()
+        for coord in self._history:
+            tmp_spc.coordinates = coord.to("cart")
+            tmp_spc.energy = coord.e
+            tmp_spc.print_xyz_file(
+                filename=filename,
+                append=True,
+            )
+        return None
+
 
-class _OptimiserHistory(list):
+class OptimiserHistory(UserList):
     """Sequential history of coordinates"""
 
     @property

autode/opt/optimisers/crfo.py

@@ -99,7 +99,7 @@ def _g_norm(self) -> GradientRMS:
 
     def _initialise_run(self) -> None:
         """Initialise the optimisation"""
-        logger.info("Initialising constrained optimisation")
+        logger.info("Initialising optimisation")
 
         self._build_internal_coordinates()
         self._coords.update_h_from_cart_h(self._low_level_cart_hessian)

autode/opt/optimisers/hessian_update.py

@@ -451,5 +451,130 @@ def conditions_met(self) -> bool:
         return True
 
 
+class FlowchartUpdate(HessianUpdater):
+    """
+    A hybrid update scheme combining BFGS, SR1 and PSB Hessian update
+    formulae. Proposed in  A. B. Birkholz and H. B. Schlegel in
+    Theor. Chem. Acc., 135 (84), 2016. This implementation is slightly
+    modified.
+    """
+
+    def __repr__(self):
+        return "Flowchart"
+
+    @property
+    def _updated_h(self) -> np.ndarray:
+        """
+        Flowchart (or FlowPSB) Hessian update scheme, that dynamically
+        switches between BFGS and SR1 depending on some criteria.
+
+        Alternatively switches to PSB update as a fallback if none of
+        the criteria are satisfied. Notation follows A. B. Birkholz,
+        H. B. Schlegel, Theor. Chem. Acc., 135 (84), 2016.
+
+        Returns:
+            (np.ndarray): H
+        """
+        z = self.y - np.matmul(self.h, self.s)
+        sr1_criteria = np.dot(z, self.s) / (
+            np.linalg.norm(z) * np.linalg.norm(self.s)
+        )
+        bfgs_criteria = np.dot(self.y, self.s) / (
+            np.linalg.norm(self.y) * np.linalg.norm(self.s)
+        )
+        if sr1_criteria < -0.1:
+            h_new = self.h + np.outer(z, z) / np.dot(z, self.s)
+            return h_new
+        elif bfgs_criteria > 0.1:
+            bfgs_delta_h = np.outer(self.y, self.y) / np.dot(self.y, self.s)
+            bfgs_delta_h -= np.linalg.multi_dot(
+                (self.h, self.s.reshape(-1, 1), self.s.reshape(1, -1), self.h)
+            ) / np.linalg.multi_dot(
+                (self.s.flatten(), self.h, self.s.flatten())
+            )
+            h_new = self.h + bfgs_delta_h
+            return h_new
+        else:
+            # Notation copied from Bofill update
+            G_i_1 = self.h  # G_{i-1}
+            dE_i = self.y - np.dot(G_i_1, self.s)  # ΔE_i = Δg_i - G_{i-1}Δx_i
+            dxTdg = np.dot(self.s, self.y)
+            G_i_PSB = (
+                G_i_1
+                + (
+                    (np.outer(dE_i, self.s) + np.outer(self.s, dE_i))
+                    / np.dot(self.s, self.s)
+                )
+                - (
+                    (
+                        (dxTdg - np.linalg.multi_dot((self.s, G_i_1, self.s)))
+                        * np.outer(self.s, self.s)
+                    )
+                    / np.dot(self.s, self.s) ** 2
+                )
+            )
+            return G_i_PSB
+
+    @property
+    def _updated_h_inv(self) -> np.ndarray:
+        """Flowchart update is only available for Hessian"""
+        return np.linalg.inv(self._updated_h)
+
+    @property
+    def conditions_met(self) -> bool:
+        """
+        Flowchart update does not have any conditions, as
+        update scheme is dynamically selected
+        """
+        return True
+
+
+class BFGSSR1Update(HessianUpdater):
+    """
+    Interpolates between BFGS and SR1 update in a fashion similar
+    to Bofill updates, but suited for minimisations. Proposed by
+    Farkas and Schlegel in J. Chem. Phys., 111, 1999, 10806
+    """
+
+    def __repr__(self):
+        return "BFGS-SR1"
+
+    @property
+    def _updated_h(self) -> np.ndarray:
+        """
+        Hybrid BFGS and SR1 update. The mixing parameter phi is defined
+        as the square root of the (1 - phi_Bofill) used in Bofill update.
+
+        Returns:
+            (np.ndarray): The updated hessian
+        """
+        bfgs_delta_h = np.outer(self.y, self.y) / np.dot(self.y, self.s)
+        bfgs_delta_h -= np.linalg.multi_dot(
+            (self.h, self.s.reshape(-1, 1), self.s.reshape(1, -1), self.h)
+        ) / np.linalg.multi_dot((self.s.flatten(), self.h, self.s.flatten()))
+
+        y_hs = self.y - np.matmul(self.h, self.s)
+        sr1_delta_h = np.outer(y_hs, y_hs) / np.dot(y_hs, self.s)
+
+        # definition according to Farkas, Schlegel, J Chem. Phys., 111, 1999
+        # NOTE: this phi is (1 - original_phi_bofill)
+        phi = np.dot(self.s, y_hs) ** 2 / (
+            np.dot(self.s, self.s) * np.dot(y_hs, y_hs)
+        )
+        sqrt_phi = np.sqrt(phi)
+        logger.info(f"BFGS-SR1 update: ϕ = {sqrt_phi:.4f}")
+        return self.h + sqrt_phi * sr1_delta_h + (1 - sqrt_phi) * bfgs_delta_h
+
+    @property
+    def _updated_h_inv(self) -> np.ndarray:
+        """For BFGS-SR1 update, only hessian is available"""
+        return np.linalg.inv(self._updated_h)
+
+    @property
+    def conditions_met(self) -> bool:
+        """No conditions need to be satisfied for BFGS-SR1 update"""
+        return True
+
+
 def _ensure_hermitian(matrix: np.ndarray) -> np.ndarray:
     return (matrix + matrix.T) / 2.0

autode/opt/optimisers/rfo.py

@@ -108,7 +108,9 @@ def _take_step_within_trust_radius(
         if len(delta_s) == 0:  # No need to sanitise a null step
             return 0.0
 
-        new_coords = self._coords + factor * delta_s
+        self._coords.allow_unconverged_back_transform = True
+        step = factor * delta_s
+        new_coords = self._coords + step
         cartesian_delta = new_coords.to("cart") - self._coords.to("cart")
         max_component = np.max(np.abs(cartesian_delta))
 
@@ -121,7 +123,8 @@ def _take_step_within_trust_radius(
             # Note because the transformation is not linear this will not
             # generate a step exactly max(∆x) ≡ α, but is empirically close
             factor = self.alpha / max_component
-            new_coords = self._coords + self.alpha / max_component * delta_s
+            step = factor * delta_s
 
-        self._coords = new_coords
+        self._coords.allow_unconverged_back_transform = False
+        self._coords = self._coords + step
         return factor