Merge pull request #201 from duartegroup/v1.3.3

V1.3.3

autode/__init__.py

@@ -38,7 +38,7 @@
   - Merge when tests pass
 """
 
-__version__ = "1.3.2"
+__version__ = "1.3.3"
 
 
 __all__ = [

autode/calculations/executors.py

@@ -360,7 +360,7 @@ def run(self) -> None:
         type_ = PRFOptimiser if self._calc_is_ts_opt else CRFOptimiser
 
         self.optimiser = type_(
-            init_alpha=0.1,  # Å
+            init_alpha=self._step_size,
             maxiter=self._max_opt_cycles,
             etol=self.etol,
             gtol=self.gtol,
@@ -433,7 +433,11 @@ def _max_opt_cycles(self) -> int:
                 if isinstance(kwd, kws.MaxOptCycles)
             )
         except StopIteration:
-            return 30
+            return 50
+
+    @property
+    def _step_size(self) -> float:
+        return 0.05 if self._calc_is_ts_opt else 0.1
 
     @property
     def _opt_trajectory_name(self) -> str:

autode/config.py

@@ -305,7 +305,7 @@ class NWChem:
             grad=[def2svp, pbe0, "task dft gradient"],
             low_sp=[def2svp, pbe0, "task dft energy"],
             opt=[def2svp, pbe0, MaxOptCycles(100), "task dft gradient"],
-            opt_ts=[def2svp, pbe0, "task dft gradient"],
+            opt_ts=[def2svp, pbe0, MaxOptCycles(50), "task dft gradient"],
             hess=[def2svp, pbe0, "task dft freq"],
             sp=[def2tzvp, pbe0, "task dft energy"],
             ecp=def2ecp,

autode/opt/coordinates/cartesian.py

@@ -89,3 +89,9 @@ def to(self, value: str) -> OptCoordinates:
             raise ValueError(
                 f"Cannot convert Cartesian coordinates to {value}"
             )
+
+    @property
+    def expected_number_of_dof(self) -> int:
+        """Expected number of degrees of freedom for the system"""
+        n_atoms = len(self.flatten()) // 3
+        return 3 * n_atoms - 6

autode/opt/coordinates/dic.py

@@ -69,6 +69,13 @@ def _calc_U(primitives: PIC, x: "CartesianCoordinates") -> np.ndarray:
         # of 3N - 6 non-redundant internals for a system of N atoms
         idxs = np.where(np.abs(lambd) > 1e-10)[0]
 
+        if len(idxs) < x.expected_number_of_dof:
+            raise RuntimeError(
+                "Failed to create a complete set of delocalised internal "
+                f"coordinates. {len(idxs)} < 3 N_atoms - 6. Likely due to "
+                f"missing primitives"
+            )
+
         logger.info(f"Removed {len(lambd) - len(idxs)} redundant vectors")
         return u[:, idxs]
 
@@ -248,6 +255,16 @@ def iadd(
     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"""
+        return list(range(len(self)))
+
+    @property
+    def inactive_indexes(self) -> List[int]:
+        """A list of indexes for the non-active modes in this coordinate set"""
+        return []
+
 
 class DICWithConstraints(DIC):
     r"""
@@ -434,7 +451,7 @@ def _schmidt_orthogonalise(arr: np.ndarray, *indexes: int) -> np.ndarray:
     provide pure primitive coordinates, which can then be constrained simply
     """
     logger.info(
-        f"Schmidt-orthogonalizing. Using {indexes} as " f"orthonormal vectors"
+        f"Schmidt-orthogonalizing. Using {indexes} as orthonormal vectors"
     )
 
     u = np.zeros_like(arr)

autode/opt/optimisers/base.py

@@ -1,13 +1,15 @@
+import os.path
 import numpy as np
 
 from abc import ABC, abstractmethod
 from typing import Union, Optional, Callable, Any
 from autode.log import logger
 from autode.utils import NumericStringDict
 from autode.config import Config
-from autode.values import GradientRMS, PotentialEnergy
+from autode.values import GradientRMS, PotentialEnergy, method_string
 from autode.opt.coordinates.base import OptCoordinates
 from autode.opt.optimisers.hessian_update import NullUpdate
+from autode.exceptions import CalculationException
 
 
 class BaseOptimiser(ABC):
@@ -213,7 +215,7 @@ def _update_gradient_and_energy(self) -> None:
         grad.clean_up(force=True, everything=True)
 
         if self._species.gradient is None:
-            raise RuntimeError(
+            raise CalculationException(
                 "Calculation failed to calculate a gradient. "
                 "Cannot continue!"
             )
@@ -227,16 +229,43 @@ def _update_hessian_gradient_and_energy(self) -> None:
         Update the energy, gradient and Hessian using the method. Will
         transform from the current coordinates type to Cartesian coordinates
         to perform the calculation, then back. Uses a numerical Hessian if
-        analytic Hessians are not implemented for this method
+        analytic Hessians are not implemented for this method. Does not
+        perform a Hessian evaluation if the molecule's energy is evaluated
+        at the same level of theory that would be used for the Hessian
+        evaluation.
 
         -----------------------------------------------------------------------
         Raises:
             (autode.exceptions.CalculationException):
         """
+        should_calc_hessian = True
+
+        if (
+            _energy_method_string(self._species)
+            == method_string(self._method, self._method.keywords.hess)
+            and self._species.hessian is not None
+        ):
+            logger.info(
+                "Have a calculated the energy at the same level of "
+                "theory as this optimisation and a present Hessian. "
+                "Not calculating a new Hessian"
+            )
+            should_calc_hessian = False
+
         self._update_gradient_and_energy()
 
+        if should_calc_hessian:
+            self._update_hessian()
+        else:
+            self._coords.update_h_from_cart_h(
+                self._species.hessian.to("Ha Å^-2")
+            )
+        return None
+
+    def _update_hessian(self) -> None:
+        """Update the Hessian of a species"""
         species = self._species.new_species(
-            name=f"{self._species.name}" f"_opt_{self.iteration}"
+            name=f"{self._species.name}_opt_{self.iteration}"
         )
         species.coordinates = self._coords.to("cartesian")
 
@@ -247,8 +276,7 @@ def _update_hessian_gradient_and_energy(self) -> None:
         )
 
         self._species.hessian = species.hessian.copy()
-        self._coords.update_h_from_cart_h(self._species.hessian)
-        return None
+        self._coords.update_h_from_cart_h(self._species.hessian.to("Ha Å^-2"))
 
     @property
     def _space_has_degrees_of_freedom(self) -> bool:
@@ -547,6 +575,10 @@ def save(self, filename: str) -> None:
             f" maxiter = {self._maxiter}"
         )
 
+        if os.path.exists(filename):
+            logger.warning(f"FIle {filename} existed. Overwriting")
+            open(filename, "w").close()
+
         for i, coordinates in enumerate(self._history):
 
             energy = coordinates.e
@@ -839,3 +871,7 @@ def __call__(self, coordinates: OptCoordinates) -> Any:
 
         logger.info("Calling callback function")
         return self._f(coordinates, **self._kwargs)
+
+
+def _energy_method_string(species: "Species") -> str:
+    return "" if species.energy is None else species.energy.method_str

autode/opt/optimisers/crfo.py

@@ -117,8 +117,18 @@ def _build_internal_coordinates(self):
             )
 
         cartesian_coords = CartesianCoordinates(self._species.coordinates)
+        primitives = self._primitives
+
+        if len(primitives) < cartesian_coords.expected_number_of_dof:
+            logger.info(
+                "Had an incomplete set of primitives. Adding "
+                "additional distances"
+            )
+            for i, j in combinations(range(self._species.n_atoms), 2):
+                primitives.append(Distance(i, j))
+
         self._coords = DICWithConstraints.from_cartesian(
-            x=cartesian_coords, primitives=self._primitives
+            x=cartesian_coords, primitives=primitives
         )
         self._coords.zero_lagrangian_multipliers()
         return None
@@ -127,7 +137,15 @@ def _build_internal_coordinates(self):
     def _primitives(self) -> PIC:
         """Primitive internal coordinates in this molecule"""
         logger.info("Generating primitive internal coordinates")
-        graph = self._species.graph
+        graph = self._species.graph.copy()
+
+        # Any distance constraints should also count as 'bonds' when forming
+        # the set of primitive internal coordinates, so that there is a
+        # single molecule if those distances are approaching dissociation
+        if self._species.constraints.distance is not None:
+            logger.info("Adding distance constraints as primitives")
+            for (i, j) in self._species.constraints.distance:
+                graph.add_edge(i, j)
 
         pic = PIC()
 

autode/opt/optimisers/prfo.py

@@ -1,14 +1,19 @@
 """Partitioned rational function optimisation"""
 import numpy as np
 from autode.log import logger
-from autode.opt.coordinates import CartesianCoordinates
-from autode.opt.optimisers.rfo import RFOptimiser
+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(RFOptimiser):
+class PRFOptimiser(CRFOptimiser):
     def __init__(
-        self, init_alpha: float = 0.1, imag_mode_idx: int = 0, *args, **kwargs
+        self,
+        init_alpha: float = 0.05,
+        recalc_hessian_every: int = 10,
+        *args,
+        **kwargs,
     ):
         """
         Partitioned rational function optimiser (PRFO) using a maximum step
@@ -17,7 +22,7 @@ def __init__(
 
         -----------------------------------------------------------------------
         Arguments:
-            init_alpha: Maximum step size
+            init_alpha: Maximum step size (Å)
 
             imag_mode_idx: Index of the imaginary mode to follow. Default = 0,
                            the most imaginary mode (i.e. most negative
@@ -29,86 +34,57 @@ def __init__(
         super().__init__(*args, **kwargs)
 
         self.alpha = float(init_alpha)
-        self.imag_mode_idx = int(imag_mode_idx)
-
+        self.recalc_hessian_every = int(recalc_hessian_every)
         self._hessian_update_types = [BofillUpdate]
 
     def _step(self) -> None:
         """Partitioned rational function step"""
 
-        self._coords.h = self._updated_h()
-
-        # Eigenvalues (\tilde{H}_kk in ref [1]) and eigenvectors (V in ref [1])
-        # of the Hessian matrix
-        lmda, v = np.linalg.eigh(self._coords.h)
+        if self.should_calculate_hessian:
+            self._update_hessian()
+        else:
+            self._coords.h = self._updated_h()
 
-        if np.min(lmda) > 0:
-            raise RuntimeError(
-                "Hessian had no negative eigenvalues, cannot " "follow to a TS"
-            )
+        idxs = list(range(len(self._coords)))
 
+        b, u = np.linalg.eigh(self._coords.h[:, idxs][idxs, :])
+        n_negative_eigenvalues = sum(lmda < 0 for lmda in b)
         logger.info(
-            f"Maximising along mode {self.imag_mode_idx} with "
-            f"λ={lmda[self.imag_mode_idx]:.4f}"
-        )
-
-        # Gradient in the eigenbasis of the Hessian. egn 49 in ref. 50
-        g_tilde = np.matmul(v.T, self._coords.g)
-
-        # Initialised step in the Hessian eigenbasis
-        s_tilde = np.zeros_like(g_tilde)
-
-        # For a step in Cartesian coordinates the Hessian will have zero
-        # eigenvalues for translation/rotation - keep track of them
-        non_zero_lmda = np.where(np.abs(lmda) > 1e-8)[0]
-
-        # Augmented Hessian 1 along the imaginary mode to maximise, with the
-        # form (see eqn. 59 in ref [1]):
-        #  (\tilde{H}_11  \tilde{g}_1) (\tilde{s}_1)  =      (\tilde{s}_1)
-        #  (                         ) (           )  =  ν_R (           )
-        #  (\tilde{g}_1        0     ) (    1      )         (     1     )
-        #
-        aug1 = np.array(
-            [
-                [lmda[self.imag_mode_idx], g_tilde[self.imag_mode_idx]],
-                [g_tilde[self.imag_mode_idx], 0.0],
-            ]
+            f"∇^2E has {n_negative_eigenvalues} negative "
+            f"eigenvalue(s). Should have 1"
         )
-        _, aug1_v = np.linalg.eigh(aug1)
 
-        # component of the step along the imaginary mode is the first element
-        # of the eigenvector with the largest eigenvalue (1), scaled by the
-        # final element
-        s_tilde[self.imag_mode_idx] = aug1_v[0, 1] / aug1_v[1, 1]
+        if n_negative_eigenvalues < 1:
+            raise CalculationException("Lost imaginary (TS) mode")
 
-        # Augmented Hessian along all other modes with non-zero eigenvalues,
-        # that are also not the imaginary mode to be followed
-        non_mode_lmda = np.delete(non_zero_lmda, [self.imag_mode_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}")
 
-        # see eqn. 60 in ref. [1] for the structure of this matrix!
-        augn = np.diag(np.concatenate((lmda[non_mode_lmda], np.zeros(1))))
-        augn[:-1, -1] = g_tilde[non_mode_lmda]
-        augn[-1, :-1] = g_tilde[non_mode_lmda]
+        delta_s = np.zeros(shape=(len(idxs),))
+        delta_s -= f[0] * u[:, 0] / (b[0] - lambda_p)
 
-        _, augn_v = np.linalg.eigh(augn)
+        for j in range(1, len(idxs)):
+            delta_s -= f[j] * u[:, j] / (b[j] - lambda_n)
 
-        # The step along all other components is then the all but the final
-        # component of the eigenvector with the smallest eigenvalue (0)
-        s_tilde[non_mode_lmda] = augn_v[:-1, 0] / augn_v[-1, 0]
-
-        # Transform back from the eigenbasis with eqn. 52 in ref [1]
-        delta_s = np.matmul(v, s_tilde)
-
-        self._take_step_within_trust_radius(delta_s)
+        _ = self._take_step_within_trust_radius(delta_s)
         return None
 
     def _initialise_run(self) -> None:
         """
         Initialise running a partitioned rational function optimisation by
         setting the coordinates and Hessian
         """
+        # self._build_internal_coordinates()
         self._coords = CartesianCoordinates(self._species.coordinates).to(
             "dic"
         )
         self._update_hessian_gradient_and_energy()
         return None
+
+    @property
+    def should_calculate_hessian(self) -> bool:
+        """Should an explicit Hessian calculation be performed?"""
+        n = self.iteration
+        return n > 1 and n % self.recalc_hessian_every == 0