More comprehensive optimiser convergence checks (#351)

* optimiser criteria

* minor edits

* simplify convergence criteria

* unfinished update

* convergence criteria set

* minor update

* unfinished update

* new update

* unfinished updates

* replace tol in DHS

* improve sign usage

* unfinished bugfix

* refactor class

* update dhs class

* minor update

* fix some tests

* minor update

* minor update

* allow None values in convergence criteria

* refactor convergence params

* fix some tests

* update some tests

* fixed crfo tests

* fixed prfo tests

* more test fixes

* black fix

* fix nwchem issues

* fix xtb issues

* updated parameters

* add tests, remove redundant code

* fix bug

* bug fix 2

* pr suggestion 1

* pr suggestion 2

* pr suggestion 3

* pr suggestion 3

* pr suggestion 4

* update changelog

autode/bracket/dhs.py

@@ -15,13 +15,14 @@
 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 @@ def __init__(
             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 @@ def _initialise_run(self) -> None:
     @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 @@ def rms_tangent_grad(self) -> GradientRMS:
         # 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 @@ def __init__(
         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 @@ def __init__(
             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 @@ def __init__(
                       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 @@ def __init__(
         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 @@ def _step(self) -> None:
 
         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 @@ def _step(self) -> None:
         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"
         )
 

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 @@ def run(self) -> None:
         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)

autode/opt/coordinates/base.py

@@ -343,6 +343,16 @@ def to(self, *args, **kwargs) -> "OptCoordinates":
     def iadd(self, value: np.ndarray) -> "OptCoordinates":
         """Inplace addition of some coordinates"""
 
+    @property
+    @abstractmethod
+    def n_constraints(self) -> int:
+        """Number of constraints in these coordinates"""
+
+    @property
+    @abstractmethod
+    def n_satisfied_constraints(self) -> int:
+        """Number of constraints that are satisfied in these coordinates"""
+
     @property
     @abstractmethod
     def active_indexes(self) -> List[int]:
@@ -358,6 +368,13 @@ def active_mol_indexes(self) -> List[int]:
     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

autode/opt/coordinates/cartesian.py

@@ -9,7 +9,6 @@
 if TYPE_CHECKING:
     from autode.values import Gradient
     from autode.hessians import Hessian
-    from autode.opt.coordinates.primitives import ConstrainedPrimitive
 
 
 class CartesianCoordinates(OptCoordinates):
@@ -61,6 +60,14 @@ def _update_h_from_cart_h(self, arr: Optional["Hessian"]) -> None:
         assert self.h_or_h_inv_has_correct_shape(arr)
         self._h = None if arr is None else np.array(arr)
 
+    @property
+    def n_constraints(self) -> int:
+        return 0
+
+    @property
+    def n_satisfied_constraints(self) -> int:
+        return 0
+
     @property
     def active_indexes(self) -> List[int]:
         return list(range(len(self)))
@@ -106,6 +113,10 @@ def to(self, value: str) -> 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"""

autode/opt/coordinates/dic.py

@@ -139,6 +139,10 @@ def from_cartesian(
         logger.info(f"Transformed in      ...{time() - start_time:.4f} s")
         return dic
 
+    @property
+    def cart_proj_g(self) -> Optional[np.ndarray]:
+        return self.to("cart").g
+
     def _update_g_from_cart_g(self, arr: Optional["Gradient"]) -> None:
         """
         Updates the gradient from a calculated Cartesian gradient
@@ -422,6 +426,20 @@ def iadd(self, value: np.ndarray) -> "OptCoordinates":
 
         return super().iadd(delta_s)
 
+    @property
+    def cart_proj_g(self) -> Optional[np.ndarray]:
+        """Obtain Cartesian gradient with constraints projected out"""
+        if self.g is None:
+            return None
+        # constrained gradient with inactive terms set to zero
+        g_s = self.g
+        g_s[self.inactive_indexes] = 0.0
+        g_s = g_s[: len(self)]
+        # back to Cartesian
+        g_x = np.matmul(self.B.T, g_s)
+        assert len(g_x) == len(self.to("cart"))
+        return g_x
+
     @property
     def g(self):
         """

autode/opt/optimisers/__init__.py

@@ -1,17 +1,20 @@
-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,
 )
 
 __all__ = [
     "NDOptimiser",
+    "ConvergenceParams",
     "RFOptimiser",
     "PRFOptimiser",
     "CRFOptimiser",
+    "QAOptimiser",
     "CartesianSDOptimiser",
     "DIC_SD_Optimiser",
 ]