diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 5c391641b..eff13511e 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -22,7 +22,7 @@ jobs: fail-fast: true matrix: os: ["ubuntu-latest", "macos-latest", "windows-latest"] - python-version: ["3.8", "3.9", "3.10"] + python-version: ["3.8", "3.12"] defaults: run: @@ -31,7 +31,7 @@ jobs: steps: - uses: actions/checkout@v3 - - uses: conda-incubator/setup-miniconda@v2 + - uses: conda-incubator/setup-miniconda@v3 with: python-version: ${{ matrix.python-version }} activate-environment: test diff --git a/.github/workflows/pytest_cov.yml b/.github/workflows/pytest_cov.yml index 0c0573bca..88d9ca084 100644 --- a/.github/workflows/pytest_cov.yml +++ b/.github/workflows/pytest_cov.yml @@ -23,7 +23,7 @@ jobs: fail-fast: true matrix: os: ["ubuntu-latest", "windows-latest"] - python-version: ["3.10"] + python-version: ["3.11"] defaults: run: @@ -32,7 +32,7 @@ jobs: steps: - uses: actions/checkout@v3 - - uses: conda-incubator/setup-miniconda@v2 + - uses: conda-incubator/setup-miniconda@v3 with: python-version: ${{ matrix.python-version }} activate-environment: test @@ -56,7 +56,8 @@ jobs: run: | pytest --cov=./ --cov-report=xml - - uses: codecov/codecov-action@v3 + - uses: codecov/codecov-action@v4 with: flags: unittests fail_ci_if_error: true + token: ${{ secrets.CODECOV_TOKEN }} # required diff --git a/README.md b/README.md index 3b8e65086..558dbe879 100644 --- a/README.md +++ b/README.md @@ -95,3 +95,5 @@ If **autodE** is used in a publication please consider citing the [paper](https: - Jonathon Vandezande ([@jevandezande](https://github.com/jevandezande)) - Shoubhik Maiti ([@shoubhikraj](https://github.com/shoubhikraj)) - Daniel Hollas ([@danielhollas](https://github.com/danielhollas)) +- Nils Heunemann ([@nilsheunemann](https://github.com/NilsHeunemann)) +- Sijie Fu ([@sijiefu](https://github.com/SijieFu)) diff --git a/autode/bracket/base.py b/autode/bracket/base.py index 11af32579..60a9e9b20 100644 --- a/autode/bracket/base.py +++ b/autode/bracket/base.py @@ -122,7 +122,7 @@ def _exceeded_maximum_iteration(self) -> bool: if self._micro_iter >= self._maxiter: logger.error( f"Reached the maximum number of micro-iterations " - f"*{self._maxiter}" + f"*{self._maxiter}*" ) return True else: @@ -165,6 +165,9 @@ def _calculate( n_cores = Config.n_cores if n_cores is None else int(n_cores) self.imgpair.set_method_and_n_cores(method, n_cores) + self.imgpair.initialise_trj( + f"{self._name}_left_history.zip", f"{self._name}_right_history.zip" + ) self._initialise_run() logger.info(f"Starting {self._name} method to find transition state") @@ -205,7 +208,7 @@ def _calculate( f"iterations (optimiser steps). {self._name} is " f"{'converged' if self.converged else 'not converged'}" ) - + self.imgpair.close_trj() self.print_geometries() self.plot_energies() if self.converged and self.ts_guess is not None: diff --git a/autode/bracket/dhs.py b/autode/bracket/dhs.py index 542396450..0a9c39b5e 100644 --- a/autode/bracket/dhs.py +++ b/autode/bracket/dhs.py @@ -9,9 +9,10 @@ from typing import Tuple, Union, Optional, Any, TYPE_CHECKING from enum import Enum -from autode.values import Distance, Angle, GradientRMS +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 @@ 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 @@ -132,18 +89,12 @@ def _initialise_run(self) -> None: self._target_dist = np.linalg.norm(self.dist_vec) self._update_gradient_and_energy() - # Hack to get the Hessian update from old coordinates + # Update the Hessian from old coordinates, if exists if self._old_coords is not None and self._old_coords.h is not None: assert isinstance(self._old_coords, CartesianCoordinates) - sub_opt = DistanceConstrainedOptimiser( - pivot_point=self._coords, # any dummy coordinate will work - maxiter=20, - gtol=1.0e-4, - etol=1.0e-4, + self._coords.update_h_from_old_h( + self._old_coords, self._hessian_update_types ) - sub_opt._coords = self._old_coords - sub_opt._coords = self._coords - self._coords.h = np.array(sub_opt._updated_h()) else: # no hessian available, use low level method self._coords.update_h_from_cart_h(self._low_level_cart_hessian) @@ -154,8 +105,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 @@ -191,7 +141,9 @@ def _update_gradient_and_energy(self) -> None: super()._update_gradient_and_energy() if self.iteration != 0: assert self._coords is not None, "Must have set coordinates" - self._coords.h = self._updated_h() + self._coords.update_h_from_old_h( + self._history.penultimate, self._hessian_update_types + ) def _step(self) -> None: """ @@ -199,23 +151,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: """ @@ -301,6 +295,7 @@ def _line_search_on_sphere( # Eq (12) to (15) in J. Chem. Phys., 90, 1989, 2154 # Notation follows the publication last_coords = self._history[-2] + assert last_coords is not None p_prime = self.dist_vec g_prime_per = self._coords.g - p_prime * ( @@ -394,18 +389,25 @@ def ts_guess(self) -> Optional["Species"]: # jumping over the barrier. So we iterate through all coords energies = [] + max_e = PotentialEnergy(-np.inf) + peak_coords: Optional[CartesianCoordinates] = None for coord in self._total_history: energies.append(coord.e) - if any(x is None for x in energies): - logger.error( - "Energy values are missing in the trajectory of this" - " image-pair. Unable to obtain transition state guess" - ) - return None - peak_coords = self._total_history[np.argmax(energies)] + if coord.e is None: + logger.error( + "Energy values are missing in the trajectory of this" + " image-pair. Unable to obtain transition state guess" + ) + return None + if coord.e > max_e: + max_e = coord.e + peak_coords = coord + + assert peak_coords is not None tmp_spc.coordinates = peak_coords tmp_spc.energy = peak_coords.e - tmp_spc.gradient = peak_coords.g.reshape(-1, 3).copy() + if peak_coords.g is not None: + tmp_spc.gradient = peak_coords.g.reshape(-1, 3).copy() return tmp_spc def get_coord_by_side(self, side: ImageSide) -> CartesianCoordinates: @@ -449,14 +451,14 @@ def get_last_step_by_side(self, side: ImageSide) -> Optional[np.ndarray]: if len(hist) < 2: return None - return hist[-1] - hist[-2] + return hist.final - hist.penultimate def get_dhs_step_by_side( 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 @@ -486,33 +488,37 @@ def __init__( 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 @@ -528,21 +534,26 @@ def __init__( initial_species, final_species ) - # DHS needs to keep an extra reference to the calculation method + # DHS needs to keep an extra reference method and n_cores self._method: Optional[Method] = None + self._n_cores: Optional[int] = None + + 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 = Distance(abs(step_size), "ang") - if self._step_size > self.imgpair.dist: + 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: @@ -561,6 +572,12 @@ def _step(self) -> None: 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 @@ -582,12 +599,13 @@ def _step(self) -> None: 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, ) tmp_spc = self._species.copy() tmp_spc.coordinates = new_coord - opt.run(tmp_spc, self._method) + opt.run(tmp_spc, self._method, self._n_cores) self._micro_iter = self._micro_iter + opt.iteration # not converged can only happen if exceeded maxiter of optimiser @@ -602,6 +620,7 @@ def _step(self) -> None: # put results back into imagepair self.imgpair.put_coord_by_side(opt.final_coordinates, side) # type: ignore + opt.clean_up() return None def _calculate( @@ -615,6 +634,7 @@ def _calculate( n_cores (int): Number of cores to use for calculation """ self._method = method + self._n_cores = n_cores super()._calculate(method, n_cores) @property @@ -652,6 +672,7 @@ def _get_dhs_step(self, side: ImageSide) -> CartesianCoordinates: 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) @@ -660,7 +681,7 @@ def _get_dhs_step(self, side: ImageSide) -> CartesianCoordinates: 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 @@ -707,6 +728,7 @@ def _get_dhs_step(self, side: ImageSide) -> CartesianCoordinates: (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) @@ -717,7 +739,7 @@ def _get_dhs_step(self, side: ImageSide) -> CartesianCoordinates: 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/ieip.py b/autode/bracket/ieip.py index 8105eeefb..27d334ff5 100644 --- a/autode/bracket/ieip.py +++ b/autode/bracket/ieip.py @@ -224,7 +224,7 @@ def redistribute_imagepair( return None -class IEIPMicroImagePair(EuclideanImagePair): +class IEIPMicroIters: """ Class to carry out the micro-iterations for the i-EIP method @@ -232,32 +232,11 @@ class IEIPMicroImagePair(EuclideanImagePair): def __init__( self, - species: "Species", left_coords: CartesianCoordinates, right_coords: CartesianCoordinates, micro_step_size: Union[Distance, float], target_dist: Union[Distance, float], ): - """ - Initialise an image-pair for i-EIP micro-iterations - - Args: - species (Species): The species object - left_coords (CartesianCoordinates): Coordinates of left image, must - have defined energy, gradient and hessian - right_coords (CartesianCoordinates): Coordinates of right image, must - have defined energy, gradient and hessian - micro_step_size (Distance|float): Step size for each micro-iteration - target_dist (Distance|float): Target distance for image-pair - """ - # dummy species are required for superclass init - left_image, right_image = species.copy(), species.copy() - left_image.coordinates = left_coords - right_image.coordinates = right_coords - assert isinstance(left_coords, CartesianCoordinates) - assert isinstance(right_coords, CartesianCoordinates) - super().__init__(left_image=left_image, right_image=right_image) - # generate the Taylor expansion surface from gradient and hessian assert left_coords.g is not None and left_coords.h is not None assert right_coords.g is not None and right_coords.h is not None @@ -269,11 +248,12 @@ def __init__( ) self._micro_step = float(Distance(micro_step_size, "ang")) self._target_dist = float(Distance(target_dist, "ang")) - - @property - def ts_guess(self) -> Optional["Species"]: - """This class does not implement TS guess""" - raise NotImplementedError + self.n_micro_iters = 0 # counter + self.left_coords = left_coords.copy() + self.right_coords = right_coords.copy() + # keep this in memory to calculate how much the coords have moved + self._start_left_coords = left_coords + self._start_right_coords = right_coords def update_both_img_engrad(self) -> None: """ @@ -292,12 +272,9 @@ def update_both_img_engrad(self) -> None: return None @property - def n_hat(self): - """ - Unit vector pointing from right to left image, the vector is - reverse (negative) of how it is defined in publication [1] - """ - return self.dist_vec / np.linalg.norm(self.dist_vec) + def _n_hat(self): + dist_vec = np.array(self.left_coords - self.right_coords) + return dist_vec / np.linalg.norm(dist_vec) def _get_perpendicular_micro_steps(self) -> List[np.ndarray]: """ @@ -314,10 +291,13 @@ def _get_perpendicular_micro_steps(self) -> List[np.ndarray]: steps = [] for coord in [self.left_coords, self.right_coords]: force = -coord.g # type: ignore - force_parall = self.n_hat * np.dot(force, self.n_hat) + force_parall = self._n_hat * np.dot(force, self._n_hat) force_perp = force - force_parall - delta_x_perp = force_perp / np.linalg.norm(force_perp) - delta_x_perp *= min(np.linalg.norm(force_perp), self._micro_step) + if np.linalg.norm(force_perp) > self._micro_step: + delta_x_perp = force_perp / np.linalg.norm(force_perp) + delta_x_perp *= self._micro_step + else: + delta_x_perp = force_perp steps.append(delta_x_perp) return steps @@ -332,13 +312,17 @@ def _get_energy_micro_steps(self) -> List[np.ndarray]: (list[np.ndarray]): A list of steps for the left and right image, in order """ - # NOTE: The sign is flipped here, because distance vector is + # NOTE: The sign is flipped here, because n_hat is # defined in the opposite direction i.e. left - right assert self.left_coords.e and self.right_coords.e - f_de = (self.left_coords.e - self.right_coords.e) / float(self.dist) - f_de = self.n_hat * float(f_de) - delta_x_e = f_de / np.linalg.norm(f_de) - delta_x_e *= min(np.linalg.norm(f_de), self._micro_step) + dist = np.linalg.norm(self.left_coords - self.right_coords) + f_de = (self.left_coords.e - self.right_coords.e) / float(dist) + f_de = self._n_hat * float(f_de) + if np.linalg.norm(f_de) > self._micro_step: + delta_x_e = f_de / np.linalg.norm(f_de) + delta_x_e *= self._micro_step + else: + delta_x_e = f_de return [delta_x_e, delta_x_e] def _get_distance_micro_steps(self): @@ -354,10 +338,14 @@ def _get_distance_micro_steps(self): # NOTE: The factor k that appears in eqn.(1) of the i-EIP paper # has been absorbed into the term in this function (i.e. the function # returns the displacements with the proper sign) - f_l = 2 * (self.dist - self._target_dist) / self.dist - f_l = -self.dist_vec * f_l - delta_x_l = f_l / np.linalg.norm(f_l) - delta_x_l *= min(np.linalg.norm(f_l), self._micro_step) + dist_vec = np.array(self.left_coords - self.right_coords) + dist = np.linalg.norm(dist_vec) + f_l = -dist_vec * 2 * (dist - self._target_dist) / dist + if np.linalg.norm(f_l) > self._micro_step: + delta_x_l = f_l / np.linalg.norm(f_l) + delta_x_l *= self._micro_step + else: + delta_x_l = f_l return [delta_x_l, -delta_x_l] def take_micro_step(self) -> None: @@ -373,45 +361,25 @@ def take_micro_step(self) -> None: left_step = perp_steps[0] + energy_steps[0] + dist_steps[0] right_step = perp_steps[1] + energy_steps[1] + dist_steps[1] - # scale the steps to have the selected micro step size - left_step /= np.linalg.norm(left_step) - left_step *= min(np.linalg.norm(left_step), self._micro_step) - right_step /= np.linalg.norm(right_step) - right_step *= min(np.linalg.norm(right_step), self._micro_step) + # scale the steps within the microiter step size + if np.linalg.norm(left_step) > self._micro_step: + left_step *= self._micro_step / np.linalg.norm(left_step) + if np.linalg.norm(right_step) > self._micro_step: + right_step *= self._micro_step / np.linalg.norm(right_step) self.left_coords = self.left_coords + left_step self.right_coords = self.right_coords + right_step - self._flush_old_matrices() + self.n_micro_iters += 1 return None - def _flush_old_matrices(self): - """The old gradient matrices should be removed to save memory""" - if self.total_iters / 2 > 2: - self._left_history[-3].g = None - self._right_history[-3].g = None - @property def max_displacement(self) -> float: - """ - The maximum displacement on either of the images compared - to the starting point (for the current macro-iteration) - - Returns: - (float): The max displacement - """ - left_displ = np.linalg.norm(self.left_coords - self._left_history[0]) + left_displ = np.linalg.norm(self.left_coords - self._start_left_coords) right_displ = np.linalg.norm( - self.right_coords - self._right_history[0] + self.right_coords - self._start_right_coords ) return max(left_displ, right_displ) - @property - def n_micro_steps(self): - """ - Total number of micro-iterations performed on this image pair - """ - return int(self.total_iters / 2) - class IEIP(BaseBracketMethod): """ @@ -538,8 +506,8 @@ def _exceeded_maximum_iteration(self) -> bool: """Whether it has exceeded the number of maximum iterations""" if self._macro_iter >= self._maxiter: logger.error( - f"Reached the maximum number of micro-iterations " - f"*{self._maxiter}" + f"Reached the maximum number of macro-iterations " + f"*{self._maxiter}*" ) return True else: @@ -572,8 +540,7 @@ def _step(self) -> None: # Turn off logging for micro-iterations logger.disabled = True assert self._target_dist is not None - micro_imgpair = IEIPMicroImagePair( - species=self._species, + micro_imgpair = IEIPMicroIters( left_coords=self.imgpair.left_coords, right_coords=self.imgpair.right_coords, micro_step_size=self._micro_step_size, @@ -581,14 +548,13 @@ def _step(self) -> None: ) while not ( - micro_imgpair.n_micro_steps >= self._max_micro_per + micro_imgpair.n_micro_iters >= self._max_micro_per or micro_imgpair.max_displacement > self._max_macro_step ): micro_imgpair.update_both_img_engrad() micro_imgpair.take_micro_step() self._micro_iter += 1 logger.disabled = False - self.imgpair.left_coords = micro_imgpair.left_coords self.imgpair.right_coords = micro_imgpair.right_coords self.imgpair.update_both_img_engrad() @@ -596,7 +562,7 @@ def _step(self) -> None: logger.info( f"Completed one i-EIP macro-iteration with " - f"{micro_imgpair.n_micro_steps} micro-iterations; maximum " + f"{micro_imgpair.n_micro_iters} micro-iterations; maximum " f"image displacement = {micro_imgpair.max_displacement:.3f}.\n" f"Left image step: {self.imgpair.last_left_step_size:.3f}, " f"Right image step: {self.imgpair.last_right_step_size:.3f}" diff --git a/autode/bracket/imagepair.py b/autode/bracket/imagepair.py index bff07055d..8023edcb0 100644 --- a/autode/bracket/imagepair.py +++ b/autode/bracket/imagepair.py @@ -2,10 +2,10 @@ Base classes for implementing all bracketing methods 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 +from typing import Optional, Tuple, TYPE_CHECKING, Union, Iterator from enum import Enum from autode.values import Distance, PotentialEnergy, Gradient @@ -14,7 +14,8 @@ from autode.neb import CINEB from autode.opt.coordinates import CartesianCoordinates from autode.opt.optimisers.hessian_update import BofillUpdate -from autode.opt.optimisers.base import OptimiserHistory +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 from autode.log import logger @@ -23,9 +24,6 @@ from autode.species import Species from autode.wrappers.methods import Method from autode.hessians import Hessian - from autode.values import Energy - -_flush_old_hessians = True def _calculate_engrad_for_species( @@ -123,7 +121,7 @@ def __init__( self._method = None self._hess_method = None self._n_cores = None - self._hessian_update_type = BofillUpdate + self._hessian_update_types = [BofillUpdate] self._left_history = OptimiserHistory() self._right_history = OptimiserHistory() @@ -195,6 +193,30 @@ def _align_species(self) -> None: rotated_p_mat = np.dot(rot_mat, p_mat.T).T self._left_image.coordinates = rotated_p_mat + def initialise_trj( + self, + left_history_name: str = "left_history_save.zip", + right_history_name: str = "right_history_save.zip", + ) -> None: + """ + Initialise the trajectory save files (history of coordinates on + left and right images) + + Args: + left_history_name: Name of savefile for left history + right_history_name: Name of savefile for right history + """ + self._left_history.open(left_history_name) + self._right_history.open(right_history_name) + + def close_trj(self): + """ + Put all coordinates in memory onto disk in the trajectory + save files, and close the trajectories + """ + self._left_history.close() + self._right_history.close() + def set_method_and_n_cores( self, method: "Method", @@ -250,6 +272,7 @@ def total_iters(self) -> int: @property def left_coords(self) -> CartesianCoordinates: """The coordinates of the left image""" + assert isinstance(self._left_history[-1], CartesianCoordinates) return self._left_history[-1] @left_coords.setter @@ -269,18 +292,16 @@ def left_coords(self, value: CartesianCoordinates): raise ValueError(f"Must have {self.n_atoms * 3} entries") if isinstance(value, CartesianCoordinates): - self._left_history.append(value.copy()) + self._left_history.add(value.copy()) else: raise TypeError self._left_image.coordinates = self.left_coords - if _flush_old_hessians: - if len(self._left_history) >= 3: - self._left_history[-3].h = None @property def right_coords(self) -> CartesianCoordinates: """The coordinates of the right image""" + assert isinstance(self._right_history[-1], CartesianCoordinates) return self._right_history[-1] @right_coords.setter @@ -300,14 +321,11 @@ def right_coords(self, value: CartesianCoordinates): raise ValueError(f"Must have {self.n_atoms * 3} entries") if isinstance(value, CartesianCoordinates): - self._right_history.append(value.copy()) + self._right_history.add(value.copy()) else: raise TypeError self._right_image.coordinates = self.right_coords - if _flush_old_hessians: - if len(self._right_history) >= 3: - self._right_history[-3].h = None @property @abstractmethod @@ -385,18 +403,10 @@ def update_both_img_hessian_by_formula(self): Update the molecular hessian for both images by update formula """ for history in [self._left_history, self._right_history]: - coords_l, coords_k = history.final, history.penultimate - assert coords_l.g is not None, "Gradient is not set!" - assert coords_l.h is None, "Hessian already exists!" - - updater = self._hessian_update_type( - h=coords_k.h, - s=coords_l.raw - coords_k.raw, - y=coords_l.g - coords_k.g, + history.final.update_h_from_old_h( + history.penultimate, self._hessian_update_types ) - coords_l.h = np.array(updater.updated_h) - return None @@ -441,56 +451,27 @@ def dist(self) -> Distance: @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: """ @@ -523,17 +504,17 @@ def run_cineb_from_end_points(self) -> None: return None @property - def _total_history(self) -> OptimiserHistory: + def _total_history(self) -> Iterator[CartesianCoordinates]: """ The total history of the image-pair, including any CI run from the endpoints """ - history = OptimiserHistory() - history.extend(self._left_history) + cineb_coords = [] if self._cineb_coords is not None: - history.append(self._cineb_coords) - history.extend(self._right_history[::-1]) # reverse order - return history + cineb_coords.append(self._cineb_coords) + return itertools.chain( + self._left_history, cineb_coords, reversed(self._right_history) + ) def print_geometries( self, @@ -550,14 +531,20 @@ def print_geometries( logger.warning("Cannot write trajectory, not enough points") return None - self._left_history.print_geometries( - species=self._left_image, filename=init_trj_filename + print_geometries_from( + self._left_history, + species=self._left_image, + filename=init_trj_filename, ) - self._right_history.print_geometries( - species=self._right_image, filename=final_trj_filename + print_geometries_from( + self._right_history, + species=self._right_image, + filename=final_trj_filename, ) - self._total_history.print_geometries( - species=self._left_image, filename=total_trj_filename + print_geometries_from( + self._total_history, + species=self._left_image, + filename=total_trj_filename, ) return None @@ -592,7 +579,8 @@ class Metrics(Enum): logger.warning("Cannot plot energies, not enough points") return None - if any(coord.e is None for coord in self._total_history): + all_energies = [coord.e for coord in self._total_history] + if any(en is None for en in all_energies): logger.error( "One or more coordinates do not have associated" " energies, unable to produce energy plot!" @@ -602,23 +590,25 @@ class Metrics(Enum): num_left_points = len(self._left_history) num_right_points = len(self._right_history) first_point = self._left_history[0] - points: List[Tuple[int, "Energy"]] = [] # list of tuples + points: list = [] # list of tuples - lowest_en = min(coord.e for coord in self._total_history) + lowest_en = min(all_energies) # type: ignore + last_coord = None for idx, coord in enumerate(self._total_history): - en = coord.e - lowest_en + en = coord.e - lowest_en # type: ignore if metric == Metrics.relative: if idx == 0: x = 0 else: - x = np.linalg.norm(coord - self._total_history[idx - 1]) + x = np.linalg.norm(coord - last_coord) x += points[idx - 1][0] # add previous distance elif metric == Metrics.from_start: x = np.linalg.norm(coord - first_point) else: # metric == Metrics.index: x = idx points.append((x, en)) + last_coord = coord left_points = points[:num_left_points] if self._cineb_coords is not None: diff --git a/autode/calculations/executors.py b/autode/calculations/executors.py index 0f5b4e10f..2e7ebe95e 100644 --- a/autode/calculations/executors.py +++ b/autode/calculations/executors.py @@ -29,6 +29,7 @@ from autode.species.species import Species from autode.wrappers.methods import Method from autode.wrappers.keywords import Keywords + from autode.opt.optimisers.base import NDOptimiser class CalculationExecutor: @@ -363,7 +364,7 @@ def run(self) -> None: type_ = PRFOptimiser if self._calc_is_ts_opt else CRFOptimiser - self.optimiser = type_( + self.optimiser: "NDOptimiser" = type_( init_alpha=self._step_size, maxiter=self._max_opt_cycles, etol=self.etol, @@ -373,10 +374,16 @@ def run(self) -> None: method.keywords.grad = kws.GradientKeywords(self.input.keywords) self.optimiser.run( - species=self.molecule, method=method, n_cores=self.n_cores + species=self.molecule, + method=method, + n_cores=self.n_cores, + name=self._opt_trajectory_name, + ) + self.optimiser.print_geometries( + self._opt_trajectory_name[:-4] + if self._opt_trajectory_name.endswith(".zip") + else self._opt_trajectory_name ) - - self.optimiser.save(self._opt_trajectory_name) if self.molecule.n_atoms == 1: return self._run_single_energy_evaluation() @@ -445,7 +452,7 @@ def _step_size(self) -> float: @property def _opt_trajectory_name(self) -> str: - return f"{self.name}_opt_trj.xyz" + return f"{self.name}_opt_trj.zip" @property def _opt_trajectory_exists(self) -> bool: @@ -462,9 +469,11 @@ def _set_properties_from_optimiser(self) -> None: if final_coords is None: raise ex.CalculationException("Final coordinates undefined") - self.molecule.coordinates = final_coords.reshape((-1, 3)) - if final_coords.g is not None: - self.molecule.gradient = final_coords.g.reshape((-1, 3)) + cart_coords = final_coords.to("cart") + self.molecule.coordinates = cart_coords.reshape((-1, 3)) + if cart_coords.g is not None: + self.molecule.gradient = cart_coords.g.reshape((-1, 3)) + self.molecule.energy = final_coords.e return None diff --git a/autode/input_output.py b/autode/input_output.py index 58c505392..cfca0ba38 100644 --- a/autode/input_output.py +++ b/autode/input_output.py @@ -92,7 +92,9 @@ def atoms_to_xyz_file( for atom in atoms: x, y, z = atom.coord # (Å) - print(f"{atom.label:<3}{x:10.5f}{y:10.5f}{z:10.5f}", file=xyz_file) + print( + f"{atom.label:<3} {x:10.5f} {y:10.5f} {z:10.5f}", file=xyz_file + ) return None diff --git a/autode/opt/coordinates/base.py b/autode/opt/coordinates/base.py index 747680e99..af2843311 100644 --- a/autode/opt/coordinates/base.py +++ b/autode/opt/coordinates/base.py @@ -12,6 +12,8 @@ from autode.units import Unit from autode.values import Gradient from autode.hessians import Hessian + from typing import Type + from autode.opt.optimisers.hessian_update import HessianUpdater class OptCoordinates(ValueArray, ABC): @@ -205,6 +207,49 @@ def update_h_from_cart_h( N atoms, zeroing the second derivatives if required""" return self._update_h_from_cart_h(arr) + def update_h_from_old_h( + self, + old_coords: "OptCoordinates", + hessian_update_types: List["Type[HessianUpdater]"], + ) -> None: + r""" + Update the Hessian :math:`H` from an old Hessian using an update + scheme. Requires the gradient to be set, and the old set of + coordinates with gradient to be available + + Args: + old_coords (OptCoordinates): Old set of coordinates with + gradient and hessian defined + hessian_update_types (list[type[HessianUpdater]]): A list of + hessian updater classes - the first updater that + meets the mathematical conditions will be used + """ + assert self._g is not None + assert isinstance(old_coords, OptCoordinates), "Wrong type!" + assert old_coords._h is not None + assert old_coords._g is not None + + for update_type in hessian_update_types: + updater = update_type( + h=old_coords._h, + s=self.raw - old_coords.raw, + y=self._g - old_coords._g, + subspace_idxs=old_coords.indexes, + ) + + if not updater.conditions_met: + logger.info(f"Conditions for {update_type} not met") + continue + + new_h = updater.updated_h + assert self.h_or_h_inv_has_correct_shape(new_h) + self._h = new_h + return None + + raise RuntimeError( + "Could not update the Hessian - no suitable update strategies" + ) + def make_hessian_positive_definite(self) -> None: """ Make the Hessian matrix positive definite by shifting eigenvalues diff --git a/autode/opt/optimisers/__init__.py b/autode/opt/optimisers/__init__.py index 0e2e37704..3b355dc79 100644 --- a/autode/opt/optimisers/__init__.py +++ b/autode/opt/optimisers/__init__.py @@ -1,9 +1,7 @@ from autode.opt.optimisers.base import NDOptimiser -from autode.opt.optimisers.bfgs import BFGSOptimiser 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, @@ -11,11 +9,9 @@ __all__ = [ "NDOptimiser", - "BFGSOptimiser", "RFOptimiser", "PRFOptimiser", "CRFOptimiser", "CartesianSDOptimiser", "DIC_SD_Optimiser", - "HybridTRMOptimiser", ] diff --git a/autode/opt/optimisers/base.py b/autode/opt/optimisers/base.py index 2a76058d1..0e3468c64 100644 --- a/autode/opt/optimisers/base.py +++ b/autode/opt/optimisers/base.py @@ -1,12 +1,23 @@ -import os.path +import os +import pickle + import numpy as np from abc import ABC, abstractmethod -from collections import UserList -from typing import Type, List, Union, Optional, Callable, Any, TYPE_CHECKING +from zipfile import ZipFile, is_zipfile +from collections import deque +from typing import ( + Type, + List, + Union, + Optional, + Callable, + Any, + TYPE_CHECKING, + Iterator, +) from autode.log import logger -from autode.utils import NumericStringDict from autode.config import Config from autode.values import GradientRMS, PotentialEnergy, method_string from autode.opt.coordinates.base import OptCoordinates @@ -37,9 +48,6 @@ def last_energy_change(self) -> PotentialEnergy: def run(self, *args: Any, **kwargs: Any) -> None: raise NotImplementedError - def save(self, filename: str) -> None: - raise NotImplementedError - @property def final_coordinates(self): raise NotImplementedError @@ -111,11 +119,20 @@ def optimise( >>> Optimiser.optimise(mol,method=ade.methods.ORCA()) """ + @property + def optimiser_params(self) -> dict: + """ + The parameters which are needed to intialise the optimiser and + will be saved in the optimiser trajectory + """ + return {"maxiter": self._maxiter} + def run( self, species: "Species", method: "Method", n_cores: Optional[int] = None, + name: Optional[str] = None, ) -> None: """ Run the optimiser. Updates species.atoms and species.energy @@ -130,6 +147,8 @@ def run( n_cores: Number of cores to use for calculations. If None then use autode.Config.n_cores + + name: The name of the optimisation save file """ self._n_cores = n_cores if n_cores is not None else Config.n_cores self._initialise_species_and_method(species, method) @@ -139,6 +158,11 @@ def run( logger.info("Optimisation is in a 0D space – terminating") return None + if name is None: + name = f"{self._species.name}_opt_trj.zip" + + self._history.open(filename=name) + self._history.save_opt_params(self.optimiser_params) self._initialise_run() logger.info( @@ -158,6 +182,7 @@ def run( break logger.info(f"Converged: {self.converged}, in {self.iteration} cycles") + self._history.close() return None @property @@ -318,7 +343,7 @@ def _coords(self) -> Optional[OptCoordinates]: logger.warning("Optimiser had no history, thus no coordinates") return None - return self._history[-1] + return self._history.final @_coords.setter def _coords(self, value: Optional[OptCoordinates]) -> None: @@ -337,7 +362,7 @@ def _coords(self, value: Optional[OptCoordinates]) -> None: return elif isinstance(value, OptCoordinates): - self._history.append(value.copy()) + self._history.add(value.copy()) else: raise ValueError( @@ -436,9 +461,6 @@ def last_energy_change(self) -> PotentialEnergy: def run(self, **kwargs: Any) -> None: pass - def save(self, filename: str) -> None: - pass - @property def final_coordinates(self): raise RuntimeError("A NullOptimiser has no coordinates") @@ -524,6 +546,11 @@ def etol(self, value: Union[int, float, PotentialEnergy]): self._etol = PotentialEnergy(value) + @property + def optimiser_params(self): + """Optimiser params to save""" + return {"maxiter": self._maxiter, "gtol": self.gtol, "etol": self.etol} + @classmethod def optimise( cls, @@ -593,93 +620,27 @@ def converged(self) -> bool: if self._abs_delta_e < self.etol / 10: logger.warning( f"Energy change is overachieved. " - f'{self.etol.to("kcal")/10:.3E} kcal mol-1. ' + f'{self.etol.to("kcal") / 10:.3E} kcal mol-1. ' f"Signaling convergence" ) return True return self._abs_delta_e < self.etol and self._g_norm < self.gtol - def save(self, filename: str) -> None: + def clean_up(self) -> None: """ - Save the entire state of the optimiser to a file + Clean up by removing the trajectory file on disk """ - assert self._species is not None, "Must have a species to save" - - if len(self._history) == 0: - logger.warning("Optimiser did no steps. Not saving a trajectory") - return None - - atomic_symbols = self._species.atomic_symbols - title_str = ( - f" etol = {self.etol.to('Ha')} Ha" - f" gtol = {self.gtol.to('Ha Å^-1')} Ha Å^-1" - 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 - cart_coords = coordinates.to("cartesian").reshape((-1, 3)) - gradient = cart_coords.g.reshape((-1, 3)) - - n_atoms = len(atomic_symbols) - assert n_atoms == len(cart_coords) == len(gradient) - - with open(filename, "a") as file: - print( - n_atoms, - f"E = {energy} Ha" + (title_str if i == 0 else ""), - sep="\n", - file=file, - ) - - for j, symbol in enumerate(atomic_symbols): - x, y, z = cart_coords[j] - dedx, dedy, dedz = gradient[j] - - print( - f"{symbol:<3}{x:10.5f}{y:10.5f}{z:10.5f}" - f"{dedx:15.5f}{dedy:10.5f}{dedz:10.5f}", - file=file, - ) - return None + self._history.clean_up() @classmethod def from_file(cls, filename: str) -> "NDOptimiser": """ - Create an optimiser from a file i.e. reload a saved state + Create an optimiser from a trajectory file i.e. reload a saved state """ - from autode.opt.coordinates.cartesian import CartesianCoordinates - - lines = open(filename, "r").readlines() - n_atoms = int(lines[0].split()[0]) - - title_line = NumericStringDict(lines[1]) - optimiser = cls( - maxiter=int(title_line["maxiter"]), - gtol=GradientRMS(title_line["gtol"]), - etol=PotentialEnergy(title_line["etol"]), - ) - - for i in range(0, len(lines), n_atoms + 2): - raw_coordinates = np.zeros(shape=(n_atoms, 3)) - gradient = np.zeros(shape=(n_atoms, 3)) - - for j, line in enumerate(lines[i + 2 : i + n_atoms + 2]): - _, x, y, z, dedx, dedy, dedz = line.split() - raw_coordinates[j, :] = [float(x), float(y), float(z)] - gradient[j, :] = [float(dedx), float(dedy), float(dedz)] - - coords = CartesianCoordinates(raw_coordinates) - coords.e = NumericStringDict(lines[i + 1])["E"] - coords.g = gradient.flatten() - - optimiser._history.append(coords) - + hist = OptimiserHistory.load(filename) + optimiser = cls(**hist.get_opt_params()) + optimiser._history = hist return optimiser @property @@ -747,82 +708,6 @@ def _log_convergence(self) -> None: logger.info(log_string) return None - def _updated_h_inv(self) -> np.ndarray: - r""" - Update the inverse of the Hessian matrix :math:`H^{-1}` for the - current set of coordinates. If the first iteration then use the true - inverse of the (estimated) Hessian, otherwise update the inverse - using a viable update strategy. - - - .. math:: - - H_{l - 1}^{-1} \rightarrow H_{l}^{-1} - - """ - assert self._coords is not None, "Must have coordinates to get H" - - if self.iteration == 0: - logger.info("First iteration so using exact inverse, H^-1") - return np.linalg.inv(self._coords.h) - - return self._best_hessian_updater.updated_h_inv - - def _updated_h(self) -> np.ndarray: - r""" - Update the Hessian matrix :math:`H` for the current set of - coordinates. If the first iteration then use the initial Hessian - - .. math:: - - H_{l - 1} \rightarrow H_{l} - - """ - assert self._coords is not None, "Must have coordinates to get H" - - if self.iteration == 0: - logger.info("First iteration so not updating the Hessian") - return self._coords.h - - return self._best_hessian_updater.updated_h - - @property - def _best_hessian_updater(self) -> "HessianUpdater": - """ - Find the best Hessian update strategy by enumerating all the possible - Hessian update types implemented for this optimiser and returning the - first that meets the criteria to be used. - - ----------------------------------------------------------------------- - Returns: - (autode.opt.optimisers.hessian_update.HessianUpdater): - - Raises: - (RuntimeError): If no suitable strategies are found - """ - coords_l, coords_k = self._history.final, self._history.penultimate - assert coords_k.g is not None and coords_l.g is not None - - for update_type in self._hessian_update_types: - updater = update_type( - h=coords_k.h, - h_inv=coords_k.h_inv, - s=coords_l.raw - coords_k.raw, - y=coords_l.g - coords_k.g, - subspace_idxs=coords_l.indexes, - ) - - if not updater.conditions_met: - logger.info(f"Conditions for {update_type} not met") - continue - - return updater - - raise RuntimeError( - "Could not update the inverse Hessian - no " - "suitable update strategies" - ) - def plot_optimisation( self, filename: Optional[str] = None, @@ -885,113 +770,306 @@ def print_geometries(self, filename: Optional[str] = None) -> None: if filename is None else str(filename) ) - - self._history.print_geometries(self._species, filename=filename) + print_geometries_from( + self._history, species=self._species, filename=filename + ) return None -class OptimiserHistory(UserList): - """Sequential history of coordinates""" +class OptimiserHistory: + """ + Sequential trajectory of coordinates with a maximum length for + storage on memory. Shunts data to disk if trajectory file is + opened, otherwise old coordinates more than the maximum number + are lost. + """ + + def __init__(self, maxlen: Optional[int] = 2) -> None: + self._filename: Optional[str] = None # filename with abs. path + self._memory: deque = deque(maxlen=maxlen) # coords in mem + self._maxlen = maxlen if maxlen is not None else float("inf") + self._is_closed = False # whether trajectory is open + self._len = 0 # count of total number of coords @property - def penultimate(self) -> OptCoordinates: + def final(self): """ - Last but one set of coordinates (the penultimate set) + Last set of coordinates in memory ----------------------------------------------------------------------- Returns: (OptCoordinates): """ - if len(self) < 2: + if len(self._memory) < 1: raise IndexError( - "Cannot obtain the penultimate set of " - f"coordinates, only had {len(self)}" + "Cannot obtain the final set of " + f"coordinates, memory is empty" ) - - return self[-2] + return self._memory[-1] @property - def final(self) -> OptCoordinates: + def penultimate(self): """ - Last set of coordinates + Last but one set of coordinates from memory (the penultimate set) ----------------------------------------------------------------------- Returns: (OptCoordinates): """ - if len(self) < 1: + if len(self._memory) < 2: raise IndexError( - "Cannot obtain the final set of coordinates from " - "an empty history" + "Cannot obtain the penultimate set of " + f"coordinates, only had {len(self._memory)}" ) - - return self[-1] + return self._memory[-2] @property - def minimum(self) -> OptCoordinates: + def _n_stored(self) -> int: + """Number of coordinates stored on disk""" + if self._filename is None: + return 0 + + with ZipFile(self._filename, "r") as file: + names = file.namelist() + n_coords = 0 + for name in names: + if name.startswith("coords_") and int(name[7:]) >= 0: + n_coords += 1 + return n_coords + + def __len__(self): + """How many coordinates have been put into this trajectory""" + return self._len + + def open(self, filename: str): """ - Minimum energy coordinates in the history + Initialise the trajectory file and write it on disk. + + Args: + filename (str): The name of the trajectory file, + should be .zip, and NOT a path + """ + if self._filename is not None: + raise RuntimeError("Already initialised, cannot initialise again!") + + # filename should not be a path + assert "\\" not in filename and "/" not in filename + if not filename.lower().endswith(".zip"): + filename += ".zip" + + if os.path.exists(filename): + logger.warning(f"File {filename} already exists, overwriting") + os.remove(filename) + + # get the full path so that it is robust to os.chdir + self._filename = os.path.abspath(filename) + + # write a header like file to help identify + with ZipFile(filename, "w") as file: + with file.open("ade_opt_trj", "w") as fh: + fh.write("Trajectory from autodE".encode("utf-8")) + return None + + @classmethod + def load(cls, filename: str): + """ + Reload the state of the trajectory from a file + + Args: + filename: The name of the trajectory .zip file, + could also be a relative path - ----------------------------------------------------------------------- Returns: - (OptCoordinates): + + """ + trj = cls() + if not filename.lower().endswith(".zip"): + filename += ".zip" + if not os.path.isfile(filename): + raise FileNotFoundError(f"The file {filename} does not exist!") + if not is_zipfile(filename): + raise ValueError( + f"The file {filename} is not a valid trajectory file" + ) + with ZipFile(filename, "r") as file: + names = file.namelist() + if "ade_opt_trj" not in names: + raise ValueError( + f"The file {filename} is not an autodE trajectory!" + ) + # handle paths with env vars or w.r.t. home dir + trj._filename = os.path.abspath( + os.path.expanduser(os.path.expandvars(filename)) + ) + trj._len = trj._n_stored + trj._is_closed = True + # load the last two into memory + if trj._len < 2: + load_idxs = [trj._len - 1] + else: + load_idxs = [trj._len - 2, trj._len - 1] + with ZipFile(trj._filename, "r") as file: + for idx in load_idxs: + with file.open(f"coords_{idx}") as fh: + trj._memory.append(pickle.load(fh)) + + return trj + + def clean_up(self): + """Remove the disk file associated with this history""" + os.remove(self._filename) + return None + + def save_opt_params(self, params: dict): """ - if len(self) == 0: - raise IndexError("No minimum with no history") + Save optimiser parameters given as a dict into the trajectory + savefile - return self[np.argmin([coords.e for coords in self])] + Args: + params (dict): + """ + assert isinstance(params, dict) + if self._filename is None: + raise RuntimeError("File not opened - cannot store data") + + # python's ZipFile does not allow overwriting files + with ZipFile(self._filename, "a") as file: + names = file.namelist() + if "opt_params" in names: + raise FileExistsError( + "Optimiser parameters are already stored -" + " cannot overwrite!" + ) + with file.open("opt_params", "w") as fh: + pickle.dump(params, fh, pickle.HIGHEST_PROTOCOL) - @property - def contains_energy_rise(self) -> bool: - r""" - Does this history contain a 'well' in the energy?:: + return None - | - E | ----- / <-- Does contain a rise in the energy - | \/ - |_________________ - Iteration + def get_opt_params(self) -> dict: + """ + Retrieve the stored optimiser parameters from the trajectory + file - ----------------------------------------------------------------------- Returns: - (bool): Presence of an explicit minima + (dict): Dictionary of optimiser parameters """ + if self._filename is None: + raise RuntimeError("File not opened - cannot get data") - for idx in range(1, len(self) - 1): - if self[idx].e < self[idx + 1].e: - return True + # python's ZipFile does not allow overwriting files + with ZipFile(self._filename, "r") as file: + names = file.namelist() + if "opt_params" not in names: + raise FileNotFoundError("Optimiser parameters are not found!") + with file.open("opt_params", "r") as fh: + data = pickle.load(fh) - return False + return data - def print_geometries(self, species: "Species", filename: str) -> None: + def add(self, coords: Optional["OptCoordinates"]) -> None: """ - Print geometries from this history of coordinates as a .xyz - trajectory file + Add a new set of coordinates to this trajectory Args: - species: The Species for which this coordinate history is generated - filename: Name of file + coords (OptCoordinates): The set of coordinates to be added """ - from autode.species import Species + if coords is None: + return None + elif not isinstance(coords, OptCoordinates): + raise ValueError("item added must be OptCoordinates") + + if self._is_closed: + raise RuntimeError("Cannot add to closed OptimiserHistory") - assert isinstance(species, Species) + self._len += 1 + # check if we need to push last coords to disk or can skip + if len(self._memory) < self._maxlen or self._filename is None: + self._memory.append(coords) + return None - if not filename.lower().endswith(".xyz"): - filename = filename + ".xyz" + n_stored = self._n_stored + with ZipFile(self._filename, "a") as file: + with file.open(f"coords_{n_stored}", "w") as fh: + pickle.dump(self._memory[0], fh, pickle.HIGHEST_PROTOCOL) + self._memory.append(coords) + return None - if os.path.isfile(filename): - logger.warning(f"{filename} already exists, overwriting...") - os.remove(filename) + def close(self): + """ + Close the Optimiser history by putting the coordinates still + in memory onto disk + """ + if self._filename is None: + raise RuntimeError("Cannot close - had no trajectory file!") - # take a copy so that original is not modified - tmp_spc = species.copy() - for coords in self: - tmp_spc.coordinates = coords.to("cart") - tmp_spc.energy = coords.e - tmp_spc.print_xyz_file(filename=filename, append=True) + idx = self._n_stored + with ZipFile(self._filename, "a") as file: + for coords in self._memory: + with file.open(f"coords_{idx}", "w") as fh: + pickle.dump(coords, fh, pickle.HIGHEST_PROTOCOL) + idx += 1 + self._is_closed = True return None + def __getitem__(self, item: int) -> Optional[OptCoordinates]: + """ + Access a coordinate from this trajectory, either from stored + data on disk, or from the memory. Only returns Cartesian + coordinates to ensure type consistency. + + Args: + item (int): Must be integer and not a slice + + Returns: + (CartesianCoordinates|None): The coordinates if found, None + if the file does not exist and coordinate is not + in memory + + Raises: + NotImplementedError: If slice is used + IndexError: If requested index does not exist + """ + if isinstance(item, slice): + raise NotImplementedError + elif isinstance(item, int): + pass + else: + raise ValueError("Index has to be type int") + + if item < 0: + item += self._len + if item < 0 or item >= self._len: + raise IndexError("Array index out of range") + + # read directly from memory if possible + if item >= (self._len - self._maxlen): + return self._memory[item - self._len] + + # have to read from disk now, return None if no file + if self._filename is None: + return None + + with ZipFile(self._filename, "r") as file: + with file.open(f"coords_{item}") as fh: + coords = pickle.load(fh) + + return coords + + def __iter__(self): + """ + Iterate through the coordinates of this trajectory + """ + for i in range(len(self)): + yield self[i] + + def __reversed__(self): + """ + Reversed iteration through the coordinates + """ + for i in reversed(range(len(self))): + yield self[i] + class ExternalOptimiser(BaseOptimiser, ABC): @property @@ -1024,3 +1102,37 @@ def __call__(self, coordinates: Optional[OptCoordinates]) -> Any: def _energy_method_string(species: "Species") -> str: return "" if species.energy is None else species.energy.method_str + + +def print_geometries_from( + coords_trj: Union[Iterator[OptCoordinates], OptimiserHistory], + species: "Species", + filename: str, +) -> None: + """ + Print geometries from an iterator over a series of coordinates + + Args: + coords_trj: The iterator for coordinates, can be OptimiserHistory + species: The Species for which the coordinate history is generated + filename: Name of file + """ + from autode.species import Species + + assert isinstance(species, Species) + + if not filename.lower().endswith(".xyz"): + filename = filename + ".xyz" + + if os.path.isfile(filename): + logger.warning(f"{filename} already exists, overwriting...") + os.remove(filename) + + # take a copy so that original is not modified + tmp_spc = species.copy() + for coords in coords_trj: + tmp_spc.coordinates = coords.to("cart") + tmp_spc.energy = coords.e + tmp_spc.print_xyz_file(filename=filename, append=True) + + return None diff --git a/autode/opt/optimisers/bfgs.py b/autode/opt/optimisers/bfgs.py deleted file mode 100644 index 7e7614831..000000000 --- a/autode/opt/optimisers/bfgs.py +++ /dev/null @@ -1,91 +0,0 @@ -import numpy as np -from abc import ABC -from typing import Type - -from autode.log import logger -from autode.values import GradientRMS, PotentialEnergy -from autode.opt.optimisers.base import NDOptimiser -from autode.opt.optimisers.hessian_update import BFGSUpdate, NullUpdate -from autode.opt.optimisers.line_search import ( - LineSearchOptimiser, - ArmijoLineSearch, -) - - -class BFGSOptimiser(NDOptimiser, ABC): - def __init__( - self, - maxiter: int, - gtol: GradientRMS, - etol: PotentialEnergy, - init_alpha: float = 1.0, - line_search_type: Type[LineSearchOptimiser] = ArmijoLineSearch, - **kwargs, - ): - """ - Broyden–Fletcher–Goldfarb–Shanno optimiser. Implementation taken - from: https://tinyurl.com/526yymsw - - ---------------------------------------------------------------------- - Arguments: - init_alpha (float): Length of the initial step to take in the line - search. Units of distance - - See Also: - - :py:meth:`NDOptimiser ` - """ - super().__init__(maxiter=maxiter, gtol=gtol, etol=etol, **kwargs) - - self._line_search_type = line_search_type - self._hessian_update_types = [BFGSUpdate, NullUpdate] - self._alpha = init_alpha - - def _step(self) -> None: - r""" - Perform a BFGS step. Requires an initial guess of the Hessian matrix - i.e. (self._coords.h must be defined). Steps follow: - - 1. Determine the inverse Hessian: :py:meth:`h_inv ` - - 2. Determine the search direction with: - - .. math:: - - \boldsymbol{p}_k = - H_k^{-1} \nabla E - - where H is the Hessian matrix, p is the search direction and - :math:`\nabla E` is the gradient of the energy with respect to the - coordinates. On the first iteration :math:`H_0` is either the true - or exact Hessian. - - 3. Performing a (in)exact line search to obtain a suitable step size - - .. math:: - - \alpha_k = \text{arg min} E(X_k + \alpha \boldsymbol{p}_k) - - and setting :math:`s_k = \alpha \boldsymbol{p}_k`, and updating the - positions accordingly (:math:`X_{k+1} = X_{k} + s_k`) - """ - assert ( - self._coords is not None - and self._coords.g is not None - and self._species is not None - and self._method is not None - ) - - self._coords.h_inv = self._updated_h_inv() - p = np.matmul(self._coords.h_inv, -self._coords.g) - - logger.info("Performing a line search") - ls = self._line_search_type( - direction=p, init_alpha=self._alpha, coords=self._coords.copy() - ) - - ls.run(self._species, self._method, n_cores=self._n_cores) - assert ls.minimum_e_coords is not None - - self._coords = ls.minimum_e_coords.copy() - return None diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py index 038f99666..a28aa8565 100644 --- a/autode/opt/optimisers/crfo.py +++ b/autode/opt/optimisers/crfo.py @@ -44,7 +44,11 @@ def _step(self) -> None: """Partitioned rational function step""" assert self._coords is not None, "Must have coords to take a step" assert self._coords.g is not None, "Must have a gradient" - self._coords.h = self._updated_h() + if self.iteration != 0: + self._coords.update_h_from_old_h( + self._history.penultimate, self._hessian_update_types + ) + assert self._coords.h is not None n, m = len(self._coords), self._coords.n_constraints logger.info(f"Optimising {n} coordinates and {m} lagrange multipliers") diff --git a/autode/opt/optimisers/line_search.py b/autode/opt/optimisers/line_search.py deleted file mode 100644 index 9b905d1fb..000000000 --- a/autode/opt/optimisers/line_search.py +++ /dev/null @@ -1,290 +0,0 @@ -""" -Line search optimisers used to solve the 1D optimisation problem by -taking steps :math:`X_{i+1} = X_i + \\alpha p` , where α is a step size and p is -a search direction. See e.g. -https://www.numerical.rl.ac.uk/people/nimg/oumsc/lectures/uepart2.2.pdf -""" -import numpy as np - -from abc import ABC, abstractmethod -from typing import Optional, Any, TYPE_CHECKING - -from autode.log import logger -from autode.opt.coordinates.cartesian import CartesianCoordinates -from autode.opt.optimisers.base import Optimiser - -if TYPE_CHECKING: - from autode.opt.coordinates import OptCoordinates - from autode.species.species import Species - from autode.wrappers.methods import Method - - -class LineSearchOptimiser(Optimiser, ABC): - """Optimiser for a 1D line search in a direction""" - - def __init__( - self, - maxiter: int = 10, - direction: Optional[np.ndarray] = None, - coords: Optional["OptCoordinates"] = None, - init_alpha: float = 1.0, - ): - """ - Line search optimiser - - ----------------------------------------------------------------------- - Keyword Arguments: - maxiter (int): Maximum number of iterations to perform - - direction (np.ndarray | None): Direction which to move the - coordinates. Shape must be broadcastable to the - coordinates. If None then will guess a sensible direction - - coords (autode.opt.coordinates.OptCoordinates | None): Initial - coordinates. If None then they will be initialised from the - species at runtime - - init_alpha (float): Initial step size - """ - super().__init__(maxiter=maxiter, coords=coords) - - self.p: Optional[np.ndarray] = direction # Direction - self.alpha = init_alpha # Step size - - @classmethod - def optimise( - cls, - species: "Species", - method: "Method", - n_cores: Optional[int] = None, - coords: Optional["OptCoordinates"] = None, - direction: Optional[np.ndarray] = None, - maxiter: int = 5, - **kwargs: Any, - ) -> None: - """ - Optimise a species along a single direction. If the direction is - unspecified then guess the direction as just the steepest decent - direction. - """ - - optimiser = cls(maxiter=maxiter, direction=direction, coords=coords) - optimiser.run(species, method, n_cores=n_cores) - - return None - - @abstractmethod - def _initialise_coordinates(self) -> None: - """Initialise the coordinates, if not already specified""" - - def _initialise_run(self) -> None: - """ - Initialise running the line search. Allows for both the coordinates - and search direction to be unspecified. - """ - if self._coords is None: - self._initialise_coordinates() - - if self.p is None: - logger.warning( - "Line search optimiser was initialised without a " - "search direction. Using steepest decent direction" - ) - assert self._coords is not None - if self._coords.g is None: - self._update_gradient_and_energy() - - assert self._coords.g is not None, "Gradient must be set" - self.p = -self._coords.g - - return None - - @property - def minimum_e_coords( - self, - ) -> Optional["OptCoordinates"]: - """Minimum energy coordinates""" - return None if len(self._history) == 0 else self._history.minimum - - @property - def _init_coords( - self, - ) -> Optional["OptCoordinates"]: - """ - Initial coordinates from which this line search was initialised from - - ----------------------------------------------------------------------- - Returns: - (autode.opt.coordinates.OptCoordinates | None): - """ - return None if len(self._history) == 0 else self._history[0] - - -class ArmijoLineSearch(LineSearchOptimiser): - def __init__( - self, - maxiter: int = 10, - direction: Optional[np.ndarray] = None, - beta: float = 0.1, - tau: float = 0.5, - init_alpha: float = 1.0, - coords: Optional["OptCoordinates"] = None, - ): - """ - Backtracking line search by Armijo. Reduces the step size iteratively - until the convergence condition is satisfied - - [1] L. Armijo. Pacific J. Math. 16, 1966, 1. DOI:10.2140/pjm.1966.16.1. - - ---------------------------------------------------------------------- - Arguments: - maxiter (int): Maximum number of iteration to perform. Should be - small O(1) - - Keyword Arguments: - direction (np.ndarray): Direction to search along - - beta (float): β parameter in the line search - - tau (float): τ parameter. Multiplicative factor when reducing the - step size. - - init_alpha (float): α_0 parameter. Initial value of the step size. - """ - super().__init__(maxiter, direction=direction, coords=coords) - - self.beta = float(beta) - self.tau = float(tau) - self.alpha = float(init_alpha) - - def _step(self) -> None: - """Take a step in the line search""" - assert self._coords is not None and self.p is not None - - self.alpha *= self.tau - self._coords = self._init_coords + self.alpha * self.p - - return None - - def _initialise_coordinates(self) -> None: - """Initialise the coordinates if they are not defined already. - Defaults to CartesianCoordinates""" - assert self._species is not None - self._coords = CartesianCoordinates(self._species.coordinates) - return None - - @property - def _satisfies_wolfe1(self) -> bool: - """First Wolfe condition:""" - assert ( - self._coords is not None - and self._init_coords is not None - and self.p is not None - ) - term_2 = self.alpha * self.beta * np.dot(self._init_coords.g, self.p) - return self._coords.e < self._init_coords.e + term_2 - - @property - def converged(self) -> bool: - r""" - Is the line search converged? Defined by the Wolfe condition - - .. math:: - - f(x + \alpha^{(l)} p_k) \le f(x) + \alpha^{(l)} \beta g\cdot p - - where α is the step size at the current iteration (denoted by l) and - β is a variable parameter. The search direction p, gradient are defined - for the initial point only. See: - https://en.wikipedia.org/wiki/Wolfe_conditions - - ----------------------------------------------------------------------- - Returns: - (bool): Line search converged? - """ - return self._has_coordinates_and_gradient and self._satisfies_wolfe1 - - -class SArmijoLineSearch(ArmijoLineSearch): - """Speculative Armijo line search""" - - @property - def converged(self) -> bool: - """Is the line search converged? For the speculative search to be - converged requires at least one iteration, there to be a well in the - search and that the Wolfe condition is satisfied. - - ----------------------------------------------------------------------- - Returns: - (bool): Line search converged? - """ - - if self.iteration == 0: - return False - - if not self._history.contains_energy_rise: - logger.warning( - f"Line search has not reached a well yet, " f"{self.alpha:.4f}" - ) - return False - - return self._has_coordinates_and_gradient and self._satisfies_wolfe1 - - def _step(self) -> None: - r""" - Take a step in the speculative line search. If the energy is monotonic - decreasing in the search then take steps that of the form - - .. math:: - - \alpha = \tau \alpha_\text{init} - - where :math:`\tau > 1`. But as soon as the energy rises then switch to - :math:`\tau_{k+1} = 1/\tau_{k}` - """ - assert self._init_coords is not None and self.p is not None - - func = min if self._history.contains_energy_rise else max - - self.tau = func(self.tau, 1 / self.tau) - self.alpha *= self.tau - - self._coords = self._init_coords + self.alpha * self.p - return None - - -class NullLineSearch(LineSearchOptimiser): - r"""A dummy line search - - .. math:: - - \alpha = \alpha_\text{init} - """ - - def _initialise_coordinates(self) -> None: - """No coordinates to initialise in a null line search""" - - def _initialise_run(self) -> None: - """Nothing required to initialise a null line search""" - - def _step(self) -> None: - """No step required in a null line search""" - - @property - def minimum_e_coords( - self, - ) -> Optional["OptCoordinates"]: - """ - Minimum energy coordinates are defined to be the true step - - ----------------------------------------------------------------------- - Returns: - (OptCoordinates): Coordinates - """ - assert self._init_coords is not None and self.p is not None - return self._init_coords + self.alpha * self.p - - @property - def converged(self) -> bool: - """A null line search is converged on the first iteration""" - return True diff --git a/autode/opt/optimisers/prfo.py b/autode/opt/optimisers/prfo.py index c2e0e4de1..4b8ae4100 100644 --- a/autode/opt/optimisers/prfo.py +++ b/autode/opt/optimisers/prfo.py @@ -47,8 +47,10 @@ def _step(self) -> None: if self.should_calculate_hessian: self._update_hessian() - else: - self._coords.h = self._updated_h() + elif self.iteration != 0: + self._coords.update_h_from_old_h( + self._history.penultimate, self._hessian_update_types + ) assert self._coords.h is not None # _update_hessian must set .h idxs = list(range(len(self._coords))) diff --git a/autode/opt/optimisers/rfo.py b/autode/opt/optimisers/rfo.py index f0304a479..dbb053309 100644 --- a/autode/opt/optimisers/rfo.py +++ b/autode/opt/optimisers/rfo.py @@ -44,8 +44,11 @@ def _step(self) -> None: assert self._coords is not None and self._coords.g is not None logger.info("Taking a RFO step") - self._coords.h = self._updated_h() - + if self.iteration != 0: + self._coords.update_h_from_old_h( + self._history.penultimate, self._hessian_update_types + ) + assert self._coords.h is not None h_n, _ = self._coords.h.shape # Form the augmented Hessian, structure from ref [1], eqn. (56) diff --git a/autode/opt/optimisers/trm.py b/autode/opt/optimisers/trm.py deleted file mode 100644 index a44a48f73..000000000 --- a/autode/opt/optimisers/trm.py +++ /dev/null @@ -1,487 +0,0 @@ -""" -A hybrid RFO (Rational Function Optimisation) and -Trust Radius Model (TRM) optimiser. Based upon the -optimisers available in multiple popular QM softwares. - -Only minimiser, this is not TS search/constrained optimiser -""" -import numpy as np -from scipy.optimize import root_scalar -from itertools import combinations -from typing import TYPE_CHECKING, Tuple, Any, Union - -from autode.opt.coordinates import CartesianCoordinates, DIC, OptCoordinates -from autode.opt.coordinates.primitives import PrimitiveDistance -from autode.opt.optimisers import CRFOptimiser -from autode.opt.optimisers.hessian_update import BFGSSR1Update -from autode.exceptions import OptimiserStepError -from autode.values import GradientRMS, PotentialEnergy, Distance -from autode.log import logger - -if TYPE_CHECKING: - from autode.species.species import Species - from autode.wrappers.methods import Method - - -class HybridTRMOptimiser(CRFOptimiser): - """ - A hybrid of RFO and TRM/QA optimiser, with dynamic trust - radius. If the RFO step is smaller in magnitude than trust radius, - then the step-length is controlled by using the Trust-Radius Model - (TRM), also sometimes known as Quadratic Approximation (QA). - It falls-back on simple scalar scaling of RFO step, if the - TRM/QA method fails to converge to the required step-size. The - trust radius is updated by a modification of Fletcher's method. - - See References: - - [1] P. Culot et. al, Theor. Chim. Acta, 82, 1992, 189-205 - - [2] T. Helgaker, Chem. Phys. Lett., 182(5), 1991, 503-510 - - [3] J. T. Golab et al., Chem. Phys., 78, 1983, 175-199 - - [4] R. Fletcher, Practical Methods of Optimization, Wiley, Chichester, 1981 - """ - - def __init__( - self, - init_trust: Union[Distance, float] = 0.1, - update_trust: bool = True, - min_trust: Union[Distance, float] = 0.01, - max_trust: Union[Distance, float] = 0.3, - damp: bool = True, - *args: Any, - **kwargs: Any, - ): - """ - Hybrid RFO and TRM/QA optimiser with dynamic trust radius (and - optional damping if oscillation in energy and gradient is detected) - - --------------------------------------------------------------------- - Args: - init_trust: Initial value of trust radius in Angstrom, this determines - the maximum norm of Cartesian step size - update_trust: Whether to update the trust radius or not - min_trust: Minimum bound for trust radius (only for trust update) - max_trust: Maximum bound for trust radius (only for trust update) - damp: Whether to apply damping if oscillation is detected - *args: Additional arguments for ``NDOptimiser`` - **kwargs: Additional keyword arguments for ``NDOptimiser`` - - Keyword Args: - maxiter (int): Maximum number of iterations - gtol (GradientRMS): Tolerance on RMS(|∇E|) in Cartesian coordinates - etol (PotentialEnergy): Tolerance on |E_i+1 - E_i| - - See Also: - :py:meth:`NDOptimiser ` - """ - _ = kwargs.pop("init_alpha", None) - super().__init__(init_alpha=init_trust, *args, **kwargs) # type: ignore - - assert 0.0 < min_trust < max_trust - self._max_alpha = Distance(max_trust, "ang") - self._min_alpha = Distance(min_trust, "ang") - self._upd_alpha = bool(update_trust) - - self._damping_on = bool(damp) - self._last_damped_iteration = 0 - - self._hessian_update_types = [BFGSSR1Update] - - @property - def _g_norm(self) -> GradientRMS: - """Calculate the RMS gradient norm in Cartesian coordinates""" - assert self._coords is not None, "Must have coordinates" - grad = self._coords.to("cart").g - return GradientRMS(np.sqrt(np.average(np.square(grad)))) - - @property - def converged(self) -> bool: - if self._species is not None and self._species.n_atoms == 1: - return True # Optimisation 0 DOF is always converged - - # gradient is better indicator of stationary point - if self._g_norm <= self.gtol / 10: - logger.warning( - f"Gradient norm criteria overachieved. " - f"{self._gtol.to('Ha/ang')/10:.3f} Ha/Å. " - f"Signalling convergence" - ) - return True - - # both gradient and energy must be converged!! - return self._abs_delta_e < self.etol and self._g_norm < self.gtol - - def _initialise_species_and_method( - self, - species: "Species", - method: "Method", - ) -> None: - """ - Initialise species and method while checking that - there are no constraints in the species - """ - super()._initialise_species_and_method(species, method) - assert self._species is not None, "Must have a species now" - assert ( - not self._species.constraints.any - ), "HybridTRMOptimiser cannot work with constraints!" - return None - - def _step(self) -> None: - """ - Hybrid RFO/TRM step; if the RFO step is larger than the trust - radius, it switches to TRM model to get the best step within - the trust radius, if that calculation fails, it simply scales - back the RFO step to lie within trust radius - """ - assert self._coords is not None, "Must have coordinates" - - self._coords.h = self._updated_h() - self._update_trust_radius() - - if self._damping_on and self._is_oscillating(): - self._damped_step() - logger.info("Skipping quasi-NR step after damping") - return None - - self._coords.allow_unconverged_back_transform = True - - rfo_h_eff = self._get_rfo_minimise_h_eff(self._coords) - rfo_step, rfo_step_size = self._get_step_and_cart_size_from_h_eff( - self._coords, rfo_h_eff - ) - - if rfo_step_size < self.alpha: - logger.info("Taking a pure RFO step") - step = rfo_step - else: - try: - qa_h_eff = self._get_trm_minimise_h_eff(self._coords) - qa_step, _ = self._get_step_and_cart_size_from_h_eff( - self._coords, qa_h_eff - ) - logger.info("Taking a TRM/QA step optimised to trust radius") - step = qa_step - - except OptimiserStepError as exc: - logger.debug(f"TRM/QA step failed: {str(exc)}") - # if TRM/QA failed, take scaled RFO step - scaled_step = rfo_step * (self.alpha / rfo_step_size) - - logger.info( - "Taking a scaled RFO step approximately within " - "trust radius" - ) - step = scaled_step - - self._coords.allow_unconverged_back_transform = False - self._coords = self._coords + step - - step_size = np.linalg.norm( - self._coords.to("cart") - self._history.penultimate.to("cart") # type: ignore - ) - logger.info(f"Size of step taken (in Cartesian) = {step_size:.3f} Å") - return None - - @staticmethod - def _get_step_and_cart_size_from_h_eff( - old_coords: OptCoordinates, hess_eff: np.ndarray - ) -> Tuple[np.ndarray, float]: - """ - Obtains the Cartesian step size given a step in the current - coordinate system (e.g., Delocalised Internal Coordinates) - - Args: - old_coords (OptCoordinates): previous coordinates - hess_eff (np.ndarray): Effective (shifted) hessian - - Returns: - (float): The step size in Cartesian - """ - step = np.matmul(-np.linalg.inv(hess_eff), old_coords.g) - new_coords = old_coords + step - cart_delta = new_coords.to("cart") - old_coords.to("cart") - return step, float(np.linalg.norm(cart_delta)) - - @staticmethod - def _get_rfo_minimise_h_eff(coords) -> np.ndarray: - """ - Using current Hessian and gradient, obtain the level-shifted - Hessian that would provide a minimising RFO step - - Returns: - (np.ndarray): The level-shifted effective Hessian - """ - h_n = coords.h.shape[0] - - # form the augmented Hessian - aug_h = np.zeros(shape=(h_n + 1, h_n + 1)) - aug_h[:h_n, :h_n] = coords.h - aug_h[-1, :h_n] = coords.g - aug_h[:h_n, -1] = coords.g - - aug_h_lmda, aug_h_v = np.linalg.eigh(aug_h) - - # RFO step uses the lowest non-zero eigenvalue - if abs(aug_h_lmda[0]) < 1.0e-14: - raise OptimiserStepError( - "First mode is 0, handling such cases are not implemented" - ) - lmda = aug_h_lmda[0] - - # effective hessian = H - lambda * I - return coords.h - lmda * np.eye(h_n) - - def _get_trm_minimise_h_eff(self, coords) -> np.ndarray: - """ - Using current Hessian and gradient, get the level-shifted Hessian - for a minimising step, whose magnitude (norm) is approximately - equal to the trust radius (TRM or QA step) in Cartesian coordinates. - - Described in J. Golab, D. L. Yeager, and P. Jorgensen, - Chem. Phys., 78, 1983, 175-199 - - Args: - coords (OptCoordinates): current coordinates - - Returns: - (np.ndarray): The level-shifted Hessian for TRM/QA step - """ - h_n = coords.h.shape[0] - h_eigvals = np.linalg.eigvalsh(coords.h) - if abs(h_eigvals[0]) < 1.0e-14: - raise OptimiserStepError( - "First mode is 0, handling such cases are not implemented" - ) - first_b = h_eigvals[0] # first non-zero eigenvalue of H - - def get_internal_step_size_and_deriv(lmda): - """Get the internal coordinate step, step size and - the derivative for the given lambda""" - shifted_h = coords.h - lmda * np.eye(h_n) - inv_shifted_h = np.linalg.inv(shifted_h) - step = -inv_shifted_h @ coords.g - size = np.linalg.norm(step) - deriv = -np.linalg.multi_dot((step, inv_shifted_h, step)) - deriv = float(deriv) / size - return size, deriv, step - - def optimise_lambda_for_int_step(int_size): - """ - Given a step length in internal coords, get the - value of lambda that will give that step. - Here f(λ) represents the internal coords step size - as a function of lambda. - """ - # use bisection to ensure lambda < first_b - d = 1.0 - found_upper_lim = False - for _ in range(20): - size, _, _ = get_internal_step_size_and_deriv(first_b - d) - # find f(x) > 0, so going downhill has no risk of jumping over - if size > int_size: - found_upper_lim = True - break - d = d / 2.0 - if not found_upper_lim: - raise OptimiserStepError("Failed to find f(λ) > 0") - - lmda_guess = first_b - d - found_step_size = False - for _ in range(50): - size, der, _ = get_internal_step_size_and_deriv(lmda_guess) - if abs(size - int_size) / int_size < 0.001: # 0.1% error - found_step_size = True - break - lmda_guess -= (1 - size / int_size) * (size / der) - if not found_step_size: - raise OptimiserStepError("Failed in optimising internal step") - return lmda_guess - - last_lmda = 0.0 # initialize non-local var - - def cart_step_length_error(int_size): - """ - Deviation from trust radius in Cartesian, - given step-size in internal coordinates - """ - nonlocal last_lmda - last_lmda = optimise_lambda_for_int_step(int_size) - _, _, step = get_internal_step_size_and_deriv(last_lmda) - step_size = np.linalg.norm( - (coords + step).to("cart") - coords.to("cart") - ) - return step_size - self.alpha - - # The value of shift parameter lambda must lie within (-infinity, first_b) - # Need to find the roots of the 1D function cart_step_length_error - int_step_size = float(self.alpha) - # starting guess size - fac = 1.0 - size_min_bound = None - size_max_bound = None - found_bounds = False - for _ in range(10): - int_step_size = int_step_size * fac - err = cart_step_length_error(int_step_size) - - if err < -1.0e-3: # found where error is < 0 - fac = 2.0 # increase the step size - size_min_bound = int_step_size - - elif err > 1.0e-3: # found where error is > 0 - fac = 0.5 # decrease the step size - size_max_bound = int_step_size - - else: # found err ~ 0 already!, no need for root finding - return coords.h - last_lmda * np.eye(h_n) - - if (size_max_bound is not None) and (size_min_bound is not None): - found_bounds = True - break - - if not found_bounds: - raise OptimiserStepError( - "Unable to find bracket range for root finding" - ) - # NOTE: The secant method does not need bracketing, however - # the bracketing is done to ensure better convergence, as - # the guesses for secant method should be closer to root - - logger.debug("Using secant method to find root") - # Use scipy's root finder with secant method - res = root_scalar( - cart_step_length_error, - method="secant", - x0=size_min_bound, - x1=size_max_bound, - maxiter=20, - rtol=0.01, # 1% margin of error - ) - - if not res.converged: - raise OptimiserStepError("Unable to converge step to trust radius") - - if last_lmda > min(0, first_b): - raise OptimiserStepError("Unknown error in finding optimal lambda") - - logger.debug(f"Optimised lambda for QA step: {last_lmda}") - # use the final lambda to construct level-shifted Hessian - return coords.h - last_lmda * np.eye(h_n) - - def _update_trust_radius(self) -> None: - """ - Updates the trust radius before a geometry step - """ - assert self._coords is not None, "Must have coordinates" - - # skip on first iteration - if self.iteration < 1: - return None - - if not self._upd_alpha: - return None - - # current coord must have energy - assert self._coords.e is not None - assert np.isfinite(self.last_energy_change) - - coords_l, coords_k = self._history.final, self._history.penultimate - step = coords_l.raw - coords_k.raw - cart_step = coords_l.to("cart") - coords_k.to("cart") - cart_step_size = np.linalg.norm(cart_step) - - pred_e_change = np.dot(coords_k.g, step) - pred_e_change += 0.5 * np.linalg.multi_dot((step, coords_k.h, step)) - - # ratio between actual deltaE and predicted deltaE - ratio = self.last_energy_change / float(pred_e_change) - - if ratio < 0.25: - set_trust = 0.5 * min(self.alpha, cart_step_size) - self.alpha = max(Distance(set_trust), self._min_alpha) - elif 0.25 < ratio < 0.75: - pass - elif 0.75 < ratio < 1.25: - # if step taken is within 5% of the trust radius - if abs(cart_step_size - self.alpha) / self.alpha < 0.05: - self.alpha = min(Distance(1.414 * self.alpha), self._max_alpha) - elif 1.25 < ratio < 1.5: - pass - else: # ratio > 1.5 - set_trust = 0.5 * min(self.alpha, Distance(cart_step_size)) - self.alpha = max(set_trust, self._min_alpha) - - if (ratio < -1.0) or (ratio > 2.0): - logger.warning( - "Energy increased/decreased by large amount," - "rejecting last geometry step" - ) - # remove the last geometry from history - self._history.pop() - - logger.info( - f"Ratio = {ratio:.3f}, Current trust radius = {self.alpha:.3f} Å" - ) - return None - - def _is_oscillating(self) -> bool: - """ - Check whether the optimiser is oscillating instead of converging. - If both the energy is oscillating (i.e. up-> down or down->up) in - the last two steps, and the energy has not gone below the lowest - energy in the last 4 iterations, then it is assumed that the - optimiser is oscillating - - Returns: - (bool): True if oscillation is detected, False otherwise - """ - # allow the optimiser 4 free iterations before checking oscillation - if self.iteration - self._last_damped_iteration < 4: - return False - - # energy change two steps before i.e. -3, -2 - e_change_before_last = self._history[-2].e - self._history[-3].e - - # sign of changes should be different if E oscillating - if not self.last_energy_change * e_change_before_last < 0: - return False - - # check if energy has gone down since the last 4 iters - min_index = np.argmin([coord.e for coord in self._history]) - if min_index < (self.iteration - 4): - logger.warning( - "Oscillation detected in optimiser, energy has " - "not decreased in 4 iterations" - ) - return True - - return False - - def _damped_step(self) -> None: - """ - Take a damped step by interpolating between the last two coordinates - """ - assert self._coords is not None, "Must have coordinates" - - logger.info("Taking a damped step...") - self._last_damped_iteration = self.iteration - - # is halfway interpolation good? - new_coords_raw = (self._coords.raw + self._history[-2].raw) / 2 - damped_step = new_coords_raw - self._coords.raw - self._coords = self._coords + damped_step - return None - - -class CartesianHybridTRMOptimiser(HybridTRMOptimiser): - def _initialise_run(self) -> None: - """Initialise the optimisation with Cartesian coordinates""" - assert self._species is not None, "Must have a species" - self._coords = CartesianCoordinates(self._species.coordinates) - self._coords.update_h_from_cart_h(self._low_level_cart_hessian) - self._update_gradient_and_energy() - return None diff --git a/autode/opt/optimisers/trust_region.py b/autode/opt/optimisers/trust_region.py deleted file mode 100644 index 38fbba776..000000000 --- a/autode/opt/optimisers/trust_region.py +++ /dev/null @@ -1,447 +0,0 @@ -r""" -Trust region methods for performing optimisations, in contrast to line -search methods the direction within the 1D problem is optimised rather than -the distance in a particular direction. The sub-problem is: - -.. math:: - - \min(m_k(p)) = E_k + g_k^T p + \frac{1}{2}p^T H_k p - \quad : \quad ||p|| \le \alpha_k - -where :math:`p` is the search direction, :math:`E_k` is the energy at an -iteration :math:`k`, :math:`g_k` is the gradient and :math:`H` is the Hessian -and :math:`alpha_k` is the trust radius at step k. Notation follows -https://optimization.mccormick.northwestern.edu/index.php/Trust-region_methods -with :math:`\Delta \equiv \alpha` -""" -import numpy as np -from typing import Optional, TYPE_CHECKING -from abc import ABC, abstractmethod - -from autode.log import logger -from autode.values import GradientRMS, PotentialEnergy -from autode.opt.optimisers.base import NDOptimiser -from autode.opt import CartesianCoordinates - -if TYPE_CHECKING: - from autode.opt.coordinates import OptCoordinates - from autode.species.species import Species - from autode.wrappers.methods import Method - - -class TrustRegionOptimiser(NDOptimiser, ABC): - def __init__( - self, - maxiter: int, - gtol: GradientRMS, - etol: PotentialEnergy, - trust_radius: float, - coords: Optional["OptCoordinates"] = None, - max_trust_radius: Optional[float] = None, - eta_1: float = 0.1, - eta_2: float = 0.25, - eta_3: float = 0.75, - t_1: float = 0.25, - t_2: float = 2.0, - **kwargs, - ): - """Trust radius optimiser""" - super().__init__( - maxiter=maxiter, etol=etol, gtol=gtol, coords=coords, **kwargs - ) - - self.alpha = trust_radius - self.alpha_max = ( - max_trust_radius - if max_trust_radius is not None - else 10 * trust_radius - ) - - # Parameters for the TR optimiser - self._eta = _Eta(eta_1, eta_2, eta_3) - self._t = _T(t_1, t_2) - - self.m: Optional[float] = None # Energy estimate - self.p: Optional[np.ndarray] = None # Direction - - @classmethod - def optimise( - cls, - species: "Species", - method: "Method", - n_cores: Optional[int] = None, - coords: Optional["OptCoordinates"] = None, - maxiter: int = 5, - gtol: GradientRMS = GradientRMS(1e-3, "Ha Å-1"), - etol: PotentialEnergy = PotentialEnergy(1e-4, "Ha"), - trust_radius: float = 0.2, - **kwargs, - ) -> None: - """ - Construct and optimiser using a trust region optimiser - """ - - optimiser = cls( - maxiter=maxiter, - gtol=gtol, - etol=etol, - trust_radius=trust_radius, - coords=coords, - ) - - optimiser.run(species, method, n_cores=n_cores) - - return None - - def _step(self) -> None: - """ - Perform a TR step based on a solution to the 'sub-problem' to find a - direction which to step in, the distance for which is fixed by the - current trust region (alpha) - """ - assert self._coords is not None - - rho = self.rho - self._solve_subproblem() - - e, g, h, p = self._coords.e, self._coords.g, self._coords.h, self.p - self.m = e + np.dot(g, p) + 0.5 * np.dot(p, np.matmul(h, p)) - - if self.iteration == 0: - # First iteration, so take a normal step - self._coords = self._coords + self.p - return - - if rho < self._eta[2]: - self.alpha *= self._t[1] - - else: - if rho > self._eta[3] and self._step_was_close_to_max: - self.alpha = min(self._t[2] * self.alpha, self.alpha_max) - - else: - pass # No updated required: α_k+1 = α_k - - if rho > self._eta[1]: - self._coords = self._coords + self.p - - else: - logger.warning( - "Trust radius step did not result in a satisfactory" - " reduction. Not taking a step" - ) - self._coords = self._coords.copy() - - return None - - @property - def _step_was_close_to_max(self) -> bool: - """ - Is the current step close to the maximum allowed? - - ----------------------------------------------------------------------- - Returns: - (bool): |p| ~ α_max - """ - return np.allclose(np.linalg.norm(self.p), self.alpha) - - @abstractmethod - def _solve_subproblem(self) -> None: - """Solve the TR 'subproblem' for the ideal step to use""" - - @abstractmethod - def _update_hessian(self) -> None: - """Solve the TR 'subproblem' for the ideal step to use""" - - def _update_gradient_and_energy(self) -> None: - """ - Update the gradient and energy, along with a couple of other - properties, derived from the energy and gradient - """ - super()._update_gradient_and_energy() - self._update_hessian() - return None - - @property - def rho(self) -> float: - """ - Calculate ρ, the ratio of the actual and predicted reductions for the - previous step - - ----------------------------------------------------------------------- - Returns: - (float): ρ - """ - - if self.iteration == 0: - logger.warning( - "ρ is unknown for the 0th iteration with only one " - "energy and gradient having been evaluated" - ) - return np.inf - - if not self._last_step_updated_coordinates: - logger.warning( - f"Step {self.iteration} did not update the " - "coordinates, using ρ = 0.5" - ) - return 0.5 - - if self.m is None: - raise RuntimeError("Predicted energy update (m) undefined") - - penultimate_energy = self._history.penultimate.e - final_energy = self._history.final.e - assert final_energy is not None and penultimate_energy is not None - - true_diff: PotentialEnergy = penultimate_energy - final_energy - predicted_diff: PotentialEnergy = penultimate_energy - self.m - - return true_diff / predicted_diff - - @property - def _last_step_updated_coordinates(self) -> bool: - """ - Did the last step in the optimiser update the coordinates - - Returns: - (bool): If the last step updated the coordinates - """ - dx = np.array(self._history.final, copy=True) - np.array( - self._history.penultimate, copy=True - ) - - return np.linalg.norm(dx) > 1e-10 - - -class CauchyTROptimiser(TrustRegionOptimiser): - """Most simple trust-radius optimiser, solving the subproblem with a - cauchy point calculation""" - - def _initialise_run(self) -> None: - """Initialise a TR optimiser, so it can take the first step""" - - if self._coords is None: - assert self._species is not None - self._coords = CartesianCoordinates(self._species.coordinates) - - self._update_gradient_and_energy() - self._solve_subproblem() - return None - - def _update_hessian(self) -> None: - """Hessian is always the identity matrix""" - assert self._coords is not None - self._coords.h = np.eye(len(self._coords)) - return None - - def _solve_subproblem(self) -> None: - r""" - Solve for the optimum direction by a Cauchy point calculation - - .. math:: - - \tau = - \begin{cases} - 1 \qquad \text{ if } g^T H g <= 0\\ - \min\left( \frac{|g|^3}{\alpha g^T H g}, 1\right) - \qquad \text{otherwise} - \end{cases} - - and - - .. math:: - - p = -\tau \frac{\alpha}{|g|} g - - """ - assert self._coords is not None - g = self._coords.g - - self.p = -self.tau * (self.alpha / np.linalg.norm(g)) * g - - return None - - @property - def tau(self) -> float: - """ - Calculate τ - - ---------------------------------------------------------------------- - Returns: - (float): τ - """ - assert self._coords is not None - - e, g, h = self._coords.e, self._coords.g, self._coords.h - g_h_g = np.dot(g, np.matmul(h, g)) - - if g_h_g <= 0: - return 1.0 - else: - return min((np.linalg.norm(g) ** 3 / (self.alpha * g_h_g), 1.0)) - - -class DoglegTROptimiser(CauchyTROptimiser): - """Dogleg Method for solving the subproblem""" - - def _solve_subproblem(self) -> None: - """ - Solve the subproblem to generate a dogleg step - """ - assert self._coords is not None - tau, g, h = self.tau, self._coords.g, self._coords.h - - p_u = -(np.dot(g, g) / np.dot(g, np.matmul(h, g))) * g - - if 0 < tau <= 1: - self.p = tau * p_u - - elif 1 < tau <= 2: - raise NotImplementedError - # self.p = tau * p_u + (tau - 1) * (p_b - p_u) - - else: - raise RuntimeError(f"τ = {tau} was outside the acceptable region") - - return None - - -class CGSteihaugTROptimiser(TrustRegionOptimiser): - """Conjugate Gradient Steihaung's Method""" - - coordinate_type = CartesianCoordinates - - def __init__(self, *args, epsilon: float = 0.001, **kwargs): - """ - CG trust region optimiser - - ----------------------------------------------------------------------- - Arguments: - *args: Arguments to pass to TrustRegionOptimiser - - epsilon: ε parameter - - **kwargs: Keyword arguments to pass to TrustRegionOptimiser - """ - super().__init__(*args, **kwargs) - - self.epsilon = epsilon - - def _initialise_run(self) -> None: - """Initialise a TR optimiser, so it can take the first step""" - - if self._coords is None: - assert self._species is not None - self._coords = self.coordinate_type(self._species.coordinates) - - self._update_gradient_and_energy() - self._solve_subproblem() - return None - - def _update_hessian(self) -> None: - """Hessian is always the identity matrix""" - assert self._coords is not None - self._coords.h = np.eye(len(self._coords)) - return None - - def _discrete_optimised_p(self, z, d) -> np.ndarray: - """ - Optimise the direction of the step to take by performing an exact line - search over τ as to reduce the value of 'm' as much as possible, where - m is the predicted reduction in the energy by performing this step - - ----------------------------------------------------------------------- - Arguments: - z: Current estimate of the step - d: Direction to move the step in - - Returns: - (np.ndarray): Step direction (p) - """ - assert self._coords is not None - - e, g, h = self._coords.e, self._coords.g, self._coords.h - tau_arr, m_arr = np.linspace(0, 10, num=1000), [] - - for tau in tau_arr: - p = z + tau * d - m = e + np.dot(g, p) + 0.5 * np.dot(p, np.matmul(h, p)) - - m_arr.append(m) - - min_m_tau = tau_arr[np.argmin(m_arr)] - logger.info(f"Optimised τ={min_m_tau:.6f}") - - p = z + min_m_tau * d - step_length = np.linalg.norm(p) - - if step_length > self.alpha: - logger.warning( - f"Step size {step_length} was too large, " f"scaling down" - ) - p *= self.alpha / step_length - - return p - - def _solve_subproblem(self) -> None: - """ - Solve the subproblem for a direction - """ - assert self._coords is not None - - h = self._coords.h - z, r = 0.0, np.array(self._coords.g, copy=True) - d = -r.copy() - - if np.linalg.norm(r) < self.epsilon: - self.p = np.zeros_like(r) - return - - for j in range(100): - if np.dot(d, np.matmul(h, d)) <= 0: - self.p = self._discrete_optimised_p(z, d) - return - - alpha = np.dot(r, r) / (np.dot(d, np.matmul(h, d))) - z += alpha * d - - if np.linalg.norm(z) >= self.alpha: - self.p = self._discrete_optimised_p(z, d) - return - - r_old = r.copy() - r += alpha * np.matmul(h, d) - - if np.linalg.norm(r) < self.epsilon: - self.p = z - return - - beta = np.dot(r, r) / np.dot(r_old, r_old) - d = -r + beta * d - - raise RuntimeError("Failed to converge CG trust region solve") - - -class _ParametersIndexedFromOne: - def __getitem__(self, item): - """Internal array is indexed from 1""" - return self._arr[item - 1] - - def __init__(self, *args): - """Scalar float arguments""" - self._arr = [float(arg) for arg in args] - - -class _Eta(_ParametersIndexedFromOne): - """η parameters in the TR optimisers""" - - def __init__(self, p1, p2, p3): - super().__init__(p1, p2, p3) - - -class _T(_ParametersIndexedFromOne): - """t parameters in the TR optimisers""" - - def __init__(self, p1, p2): - super().__init__(p1, p2) diff --git a/autode/opt/optimisers/utils.py b/autode/opt/optimisers/utils.py new file mode 100644 index 000000000..936d7a2c7 --- /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/autode/reactions/reaction.py b/autode/reactions/reaction.py index 013d867a1..e24899bce 100644 --- a/autode/reactions/reaction.py +++ b/autode/reactions/reaction.py @@ -431,10 +431,17 @@ def delta(self, delta_type: str) -> Optional[Energy]: """ def delta_type_matches(*args): - return any(s in delta_type.lower() for s in args) + return any( + s + in delta_type.lower() + .replace("ddagger", "") + .replace("double dagger", "") + for s in args + ) def is_ts_delta(): - return delta_type_matches("ddagger", "‡", "double dagger") + ts_synonyms = ["ddagger", "‡", "double dagger"] + return any(s in delta_type.lower() for s in ts_synonyms) # Determine the species on the left and right-hand sides of the equation lhs: List[Species] = self.reacs diff --git a/autode/species/species.py b/autode/species/species.py index 2bff81335..411009e97 100644 --- a/autode/species/species.py +++ b/autode/species/species.py @@ -401,10 +401,14 @@ def gradient(self, value: Union[val.Gradient, np.ndarray, None]): return if hasattr(value, "shape") and value.shape != (self.n_atoms, 3): - raise ValueError( - "Could not set the gradient. Incorrect shape: " - f"{value.shape} != {(self.n_atoms, 3)}" - ) + try: + value = value.reshape((self.n_atoms, 3)) + except (ValueError, AttributeError): + raise ValueError( + "Could not set the gradient. Incorrect shape: " + f"{value.shape}. Must be either {(self.n_atoms, 3)}, " + f"or {(self.n_atoms * 3,)}" + ) if isinstance(value, val.Gradient): self._grad = value diff --git a/autode/transition_states/base.py b/autode/transition_states/base.py index 31f2e47ff..7dfbcd57f 100644 --- a/autode/transition_states/base.py +++ b/autode/transition_states/base.py @@ -463,13 +463,7 @@ def displaced_species_along_mode( disp_coords -= (disp_factor / 20) * mode_disp_coords - # Create a new species from the initial - disp_species = Species( - name=f"{species.name}_disp", - atoms=species.atoms.copy(), - charge=species.charge, - mult=species.mult, - ) + disp_species = species.new_species(name=f"{species.name}_disp") disp_species.coordinates = disp_coords return disp_species diff --git a/doc/changelog.rst b/doc/changelog.rst index 924a896a0..d70b8d361 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -1,6 +1,27 @@ Changelog ========= +1.4.3 +------ +------- + +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 +********* +- Fixes no solvent being added in QRC calculations +- Fixes cases where .xyz files are printed with no space between coordinates when the coordinate value is large. +- Fixes DHS and DHS-GS methods ignoring number of cores + +Usability improvements/Changes +****************************** +- Added consistent aliases for double dagger across all energies in :code:`autode.reaction.delta` +- Optimiser trajectories are saved on disk instead of keeping completely in memory +- Hessian updates the refactored into :code:`OptCoordinates` class + 1.4.2 ------ ------- diff --git a/setup.py b/setup.py index 300adb952..3854df28e 100644 --- a/setup.py +++ b/setup.py @@ -32,7 +32,7 @@ setup( name="autode", - version="1.4.2", + version="1.4.3", python_requires=">3.7", packages=[ "autode", diff --git a/tests/test_bracket/data/geometries.zip b/tests/test_bracket/data/geometries.zip index b84f83fa4..78a4dfaba 100644 Binary files a/tests/test_bracket/data/geometries.zip and b/tests/test_bracket/data/geometries.zip differ diff --git a/tests/test_bracket/test_dhs.py b/tests/test_bracket/test_dhs.py index 4dd8de80f..e7e39b4d0 100644 --- a/tests/test_bracket/test_dhs.py +++ b/tests/test_bracket/test_dhs.py @@ -11,9 +11,9 @@ DHS, DHSGS, DistanceConstrainedOptimiser, - TruncatedTaylor, DHSImagePair, ImageSide, + OptimiserStepError, ) from autode.opt.coordinates import CartesianCoordinates from autode import Config @@ -23,33 +23,6 @@ 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(): @@ -97,6 +70,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() @@ -151,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, ) @@ -165,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) @@ -187,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, ) @@ -195,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() @@ -230,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, ) @@ -293,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, @@ -320,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_bracket/test_ieip.py b/tests/test_bracket/test_ieip.py index ae428090c..aa64d1337 100644 --- a/tests/test_bracket/test_ieip.py +++ b/tests/test_bracket/test_ieip.py @@ -4,7 +4,7 @@ from autode.species import Molecule from autode.methods import XTB -from autode.bracket.ieip import IEIP, ElasticImagePair, IEIPMicroImagePair +from autode.bracket.ieip import IEIP, ElasticImagePair, IEIPMicroIters from autode.bracket.ieip import _calculate_low_sp_energy_for_species from autode.geom import calc_rmsd from ..testutils import requires_working_xtb_install, work_in_zipped_dir @@ -55,16 +55,17 @@ def test_low_sp_func(): @work_in_zipped_dir(datazip) def test_ieip_microiters(): micro_step_size = 1e-4 - rct = Molecule("da_reactant.xyz") - prod = Molecule("da_product.xyz") + # use almost converged images for quick calc + rct = Molecule("da_rct_image.xyz") + prod = Molecule("da_prod_image.xyz") imgpair = ElasticImagePair(rct, prod) imgpair.set_method_and_n_cores(method=XTB(), n_cores=1) imgpair.update_both_img_engrad() - imgpair.redistribute_imagepair(ll_neb_interp=False) - imgpair.update_both_img_engrad() - imgpair.update_both_img_hessian_by_calc() - micro_imgpair = IEIPMicroImagePair( - imgpair._left_image, + # load hessian from txt to save time + imgpair.left_coords.h = np.loadtxt("da_rct_image_hess.txt") + imgpair.right_coords.h = np.loadtxt("da_prod_image_hess.txt") + # micro iterations + micro_imgpair = IEIPMicroIters( imgpair.left_coords, imgpair.right_coords, micro_step_size=micro_step_size, @@ -80,9 +81,7 @@ def test_ieip_microiters(): for _ in range(3): micro_imgpair.update_both_img_engrad() micro_imgpair.take_micro_step() - # the previous gradient matrices should be flushed after three iters - assert micro_imgpair._left_history[1].g is None - assert micro_imgpair._right_history[1].g is None + assert micro_imgpair.n_micro_iters == 4 @requires_working_xtb_install diff --git a/tests/test_bracket/test_imagepair.py b/tests/test_bracket/test_imagepair.py index dd12cffb5..1f9bcdb5d 100644 --- a/tests/test_bracket/test_imagepair.py +++ b/tests/test_bracket/test_imagepair.py @@ -300,37 +300,3 @@ def test_hessian_update(): imgpair.update_both_img_hessian_by_formula() assert imgpair.left_coords.h is not None assert imgpair.right_coords.h is not None - - # calling Hessian update again will raise exception - with pytest.raises(AssertionError): - imgpair.update_both_img_hessian_by_formula() - - -@work_in_zipped_dir(datazip) -def test_imgpair_hessian_flushing(): - import autode.bracket.imagepair - - mol1 = Molecule( - atoms=[ - Atom("N", 0.5588, 0.0000, 0.0000), - Atom("N", -0.5588, 0.0000, 0.0000), - ] - ) - h = np.loadtxt("n2_hess.txt") - - imgpair = NullImagePair(mol1, mol1.copy()) - # turn off flushing Hessians - autode.bracket.imagepair._flush_old_hessians = False - imgpair.left_coords.h = h.copy() - # generate dummy coordinates - imgpair.left_coords = imgpair.left_coords * 0.99 - imgpair.left_coords = imgpair.left_coords * 0.99 - assert imgpair._left_history[0].h is not None - - # turn on flushing Hessians again - autode.bracket.imagepair._flush_old_hessians = True - imgpair = NullImagePair(mol1, mol1.copy()) - imgpair.right_coords.h = h.copy() - imgpair.right_coords = imgpair.right_coords * 0.99 - imgpair.right_coords = imgpair.right_coords * 0.99 - assert imgpair._right_history[0].h is None diff --git a/tests/test_input_output.py b/tests/test_input_output.py index 649a9f54a..2066cc643 100644 --- a/tests/test_input_output.py +++ b/tests/test_input_output.py @@ -1,3 +1,5 @@ +import sys + import numpy as np from autode.input_output import ( @@ -76,6 +78,16 @@ def test_making_xyz_file(): xyz_lines = open("test.xyz", "r").readlines() assert len(xyz_lines) == 4 + # Simulate situations when the coordinates are super large (or small) + atoms[0].coord = (sys.float_info.min, sys.float_info.max, -1e6) + atoms_to_xyz_file(atoms, filename="test.xyz", append=False) + + # check if the first atom is written correctly: the third line (index 2) + xyz_lines = open("test.xyz", "r").readlines() + assert len(xyz_lines[2].split()) == 4 and np.isclose( + float(xyz_lines[2].split()[-1]), -1e6 + ) + # With append should add the next set of atoms to the same file atoms_to_xyz_file(atoms, filename="test.xyz", append=True) diff --git a/tests/test_opt/test_bfgs.py b/tests/test_opt/test_bfgs.py deleted file mode 100644 index 0be5bab79..000000000 --- a/tests/test_opt/test_bfgs.py +++ /dev/null @@ -1,219 +0,0 @@ -""" -https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm -""" -import numpy as np -from autode.opt.coordinates.cartesian import CartesianCoordinates -from autode.opt.optimisers.bfgs import BFGSOptimiser -from autode.opt.optimisers.line_search import ( - NullLineSearch, - ArmijoLineSearch, - SArmijoLineSearch, -) -from autode.species import Molecule -from .setup import Method - - -class TestBFGSOptimiser2D(BFGSOptimiser): - """Simple 2D optimiser using a BFGS update step""" - - __test__ = False - - def __init__( - self, - e_func, - g_func, - init_x, - init_y, - maxiter=30, - etol=1e-4, - gtol=1e-3, - line_search_type=NullLineSearch, - init_alpha=1.0, - ): - super().__init__( - maxiter=maxiter, - line_search_type=line_search_type, - init_alpha=init_alpha, - etol=etol, - gtol=gtol, - ) - - init_arr = np.array([init_x, init_y]) - self._coords = CartesianCoordinates(init_arr) - - self.e_func = e_func - self.g_func = g_func - - def _space_has_degrees_of_freedom(self) -> bool: - return True - - def _log_convergence(self) -> None: - x, y = self._coords - print( - f"{x:.4f}, {y:.4f}", - f"E = {round(self._coords.e, 5)}," - # f' {np.round(self._coords.g, 5)}' - ) - - def _update_gradient_and_energy(self) -> None: - x, y = self._coords - self._coords.e = self.e_func(x, y) - self._coords.g = self.g_func(x, y) - - def _initialise_run(self) -> None: - # Guess the Hessian as the identity matrix - self._update_gradient_and_energy() - self._coords.h = np.eye(len(self._coords)) - - @property - def converged(self) -> bool: - return np.linalg.norm(self._coords.g) < self._gtol - - -class TestBFGSOptimiser(TestBFGSOptimiser2D): - def __init__(self, **kwargs): - super().__init__( - e_func=lambda x, y: x**2 + y**2, - g_func=lambda x, y: np.array([2.0 * x, 2.0 * y]), - init_y=1.0, - init_x=1.0, - **kwargs, - ) - - -def test_simple_quadratic_opt(): - optimiser = TestBFGSOptimiser() - method = Method() - assert method.is_available - - optimiser.run(Molecule(name="blank"), method=Method()) - assert optimiser.converged - - -def test_inv_hessian_update(): - # Initial imperfect guess of the Hessian matrix for E = x^2 + y^2 - init_h = np.array([[1.0, 0.0], [0.0, 1.0]]) - - optimiser = TestBFGSOptimiser() - optimiser._species, optimiser._method = Molecule(name="blank"), Method() - assert optimiser.iteration == 0 - - optimiser._update_gradient_and_energy() - optimiser._coords.h = init_h.copy() - - optimiser._step() - - # Should take the inverse of the guess after the first step - assert np.allclose(optimiser._history[-1].h_inv, np.linalg.inv(init_h)) - - # Then for the new set of coordinates generate a better guess of the - # inverse Hessian once the gradient has been updated - optimiser._update_gradient_and_energy() - - h_inv_true = np.array([[0.5, 0.0], [0.0, 0.5]]) - - assert np.linalg.norm( - optimiser._updated_h_inv() - h_inv_true - ) < np.linalg.norm(np.linalg.inv(init_h) - h_inv_true) - - optimiser.run(Molecule(name="blank"), method=Method()) - assert optimiser.converged - assert np.isclose(optimiser._coords.e, 0.0, atol=1e-5) - - # By the end of the optimisation the Hessian should be pretty good - assert np.allclose(optimiser._coords.h_inv, h_inv_true, atol=1) - - -def test_quadratic_opt(): - """E = x^2 + y^2 + xy/10, ∇E = (2x + 0.1y, 2y + 0.1x)""" - - optimiser = TestBFGSOptimiser2D( - e_func=lambda x, y: x**2 + y**2 + x * y / 10.0, - g_func=lambda x, y: np.array([2.0 * x + 0.1 * y, 2.0 * y + 0.1 * x]), - init_x=1.0, - init_y=1.0, - ) - optimiser.run(Molecule(name="blank"), method=Method()) - assert optimiser.converged - assert optimiser._coords.e < 1e-3 - assert np.allclose( - optimiser._coords, np.zeros(2), atol=1e-3 # Minimum is at (0, 0) - ) - - -def check_gaussian_well_opt(init_x, init_y): - """E = -exp(-(x^2 + y^2)), ∇E = (2xE, 2yE)""" - - def energy(x, y): - return -np.exp(-(x**2 + y**2)) - - def grad(x, y): - e = np.exp(-(x**2 + y**2)) - return np.array([2.0 * x * e, 2.0 * y * e]) - - def hessian(x, y): - h_xy = -4.0 * x * y * np.exp(-(x**2) - y**2) - return np.array( - [ - [2 * (1 - 2.0 * x**2) * np.exp(-(x**2) - y**2), h_xy], - [h_xy, 2 * (1 - 2.0 * y**2) * np.exp(-(x**2) - y**2)], - ] - ) - - class LineSearch(SArmijoLineSearch): - def _update_gradient_and_energy(self) -> None: - self._coords.e = energy(*self._coords) - self._coords.g = grad(*self._coords) - - optimiser = TestBFGSOptimiser2D( - e_func=energy, - g_func=grad, - init_x=init_x, - init_y=init_y, - line_search_type=LineSearch, - ) - - optimiser.run(Molecule(name="blank"), method=Method()) - - assert optimiser.converged - assert np.isclose(optimiser._coords.e, -1) # Minimum is -1 - assert np.allclose( - optimiser._coords, np.zeros(2), atol=1e-3 # Minimum is at (0, 0) - ) - - # Hessian should be close-ish to the true - assert np.allclose( - optimiser._coords.h, hessian(*optimiser._coords), atol=1 - ) - - -def test_gaussian_well_opt(): - for x, y in [(-0.1, 0.0), (-1.0, 0.0), (2.0, 1.0), (2.0, -2.0)]: - check_gaussian_well_opt(init_x=x, init_y=y) - - -def test_complex_2d_opt(): - def energy(x, y): - return 10 * (y - x**2) ** 2 + (x - 1) ** 2 - - def grad(x, y): - return np.array( - [2 * (20 * x**3 - 20 * x * y + x - 1), 20 * (y - x**2)] - ) - - optimiser = TestBFGSOptimiser2D( - e_func=energy, - g_func=grad, - init_x=-1.0, - init_y=1.0, - init_alpha=0.1, - maxiter=100, - ) - optimiser.run(Molecule(name="blank"), method=Method()) - - assert optimiser.converged - assert np.allclose( - optimiser._coords, - np.array([1.0, 1.0]), - atol=0.1, # Minimum is at (1, 1) - ) diff --git a/tests/test_opt/test_dimer.py b/tests/test_opt/test_dimer.py index 4f24294bf..da0caed4e 100644 --- a/tests/test_opt/test_dimer.py +++ b/tests/test_opt/test_dimer.py @@ -6,6 +6,7 @@ from autode.values import MWDistance from autode.opt.coordinates.dimer import DimerCoordinates, DimerPoint from autode.opt.optimisers.dimer import Dimer +from autode.utils import work_in_tmp_dir from .molecules import methane_mol from ..testutils import requires_working_xtb_install @@ -164,6 +165,7 @@ def _initialise_run(self) -> None: return None +@work_in_tmp_dir() def test_dimer_2d(): arr = np.array( [ @@ -186,6 +188,7 @@ def test_dimer_2d(): ) # optimise the rotation, should be able to be very accurate + dimer._history.open("dimer_test.zip") dimer._initialise_run() dimer._optimise_rotation() @@ -212,6 +215,7 @@ def test_dimer_2d(): @requires_working_xtb_install +@work_in_tmp_dir() def test_dimer_sn2(): left_point = Molecule( name="sn2_left", diff --git a/tests/test_opt/test_hessian_update.py b/tests/test_opt/test_hessian_update.py index 9594c4c45..be6952a19 100644 --- a/tests/test_opt/test_hessian_update.py +++ b/tests/test_opt/test_hessian_update.py @@ -2,6 +2,7 @@ import pytest import numpy as np from ..testutils import work_in_zipped_dir +from autode.opt.coordinates import CartesianCoordinates from autode.opt.optimisers.hessian_update import ( BFGSUpdate, BFGSPDUpdate, @@ -184,6 +185,16 @@ def test_bfgs_pd_update(eigval, expected): assert updater.conditions_met == expected +def test_update_fails_if_no_suited_scheme(): + coords1 = CartesianCoordinates(np.array([0.1, 0.1])) + coords1.g = np.array([0.1, 0.1]) + coords1.h = np.array([[1.0, 0.0], [0.0, -1]]) + coords2 = CartesianCoordinates(np.array([0.2, 0.2])) + coords2.g = np.array([0.2, 0.2]) + with pytest.raises(RuntimeError): + coords2.update_h_from_old_h(coords1, [BFGSPDUpdateNoUpdate]) + + @work_in_zipped_dir(os.path.join(here, "data", "hessians.zip")) def test_bfgs_and_bfgssr1_update_water(): h = np.loadtxt("water_cart_hessian0.txt") diff --git a/tests/test_opt/test_line_search.py b/tests/test_opt/test_line_search.py deleted file mode 100644 index f80e2b932..000000000 --- a/tests/test_opt/test_line_search.py +++ /dev/null @@ -1,158 +0,0 @@ -import shutil -import numpy as np -from autode.atoms import Atom -from autode.species import Molecule -from autode.utils import work_in_tmp_dir -from autode.methods import XTB -from autode.opt.coordinates.cartesian import CartesianCoordinates -from autode.opt.optimisers.line_search import ( - ArmijoLineSearch, - NullLineSearch, - LineSearchOptimiser, -) -from ..testutils import requires_working_xtb_install -from .setup import Method - - -def quadratic(x, y): - return x**2 + y**2, np.array([2 * x, 2 * y]) - - -class TestSDLineSearch(LineSearchOptimiser): - """Line search for E = x^2 + y^2""" - - __test__ = False - - def __init__( - self, - init_alpha=1.0, - energy_grad_func=quadratic, - direction=None, - coords=None, - ): - super().__init__( - maxiter=10, - direction=direction, - init_alpha=init_alpha, - coords=coords, - ) - - self.energy_grad_func = energy_grad_func - - @property - def converged(self) -> bool: - """Simple convergence criteria""" - return ( - self.iteration > 0 - and self._coords.e is not None - and self._coords.e < 0.01 - ) - - def _log_convergence(self) -> None: - pass # print(self._coords.e) - - def _initialise_coordinates(self) -> None: - self._coords = CartesianCoordinates(np.array([1.1, 0.2])) - - def _step(self) -> None: - self._coords = self._coords + self.alpha * self.p - - def _update_gradient_and_energy(self) -> None: - x, y = self._coords - self._coords.e, self._coords.g = self.energy_grad_func(x, y) - - -class TestArmijoLineSearch(ArmijoLineSearch): - __test__ = False - - def __init__(self, init_step_size=1.0, energy_grad_func=quadratic): - super().__init__(maxiter=10, init_alpha=init_step_size) - - self.energy_grad_func = energy_grad_func - - def _initialise_coordinates(self) -> None: - self._coords = CartesianCoordinates(np.array([-0.8, 1.0])) - - def _update_gradient_and_energy(self) -> None: - return TestSDLineSearch._update_gradient_and_energy(self) - - def _log_convergence(self) -> None: - pass # print(self._e_prev, self._species.energy) - - -def test_null_line_search(): - ls = NullLineSearch(init_alpha=0.1) - assert ls.converged - - # Dummy (empty) coordinates should be initialisable without a species - assert ls._initialise_coordinates() is None - - ls.run(Molecule(name="blank"), Method()) - assert np.isclose(ls.alpha, 0.1) - - -def test_simple_line_search(): - blank_mol = Molecule(name="blank") - blank_method = Method() - - optimiser = TestSDLineSearch(init_alpha=0.1) - optimiser.run(blank_mol, method=blank_method) - assert optimiser.converged - - # Minimum is at (0, 0). Should be close to that - assert np.allclose(optimiser._coords, np.array([0.0, 0.0]), atol=1e-1) - - -def test_armijo_line_search_default(): - optimiser = TestArmijoLineSearch() - assert not optimiser.converged - - optimiser.run(Molecule(name="blank"), method=Method()) - assert optimiser.converged - - # Minimum is at (0, 0). Should be close to that point with 0 energy - assert np.allclose(optimiser._coords, np.array([0.0, 0.0])) - assert np.isclose(optimiser._coords.e, 0.0) - - -def test_armijo_line_search_diff_step_sizes(): - # using different step sizes should also converge - for init_step_size in (0.1, 0.5, 1.0, 2.0, 4.0, 10.0): - optimiser = TestArmijoLineSearch(init_step_size=init_step_size) - optimiser.run(Molecule(name="blank"), method=Method()) - - assert optimiser.converged - - -def test_armijo_line_search_complex_func(): - def energy_grad(x, y): - energy = 10 * (y - x**2) ** 2 + (x - 1) ** 2 - gradient = np.array( - [2 * (20 * x**3 - 20 * x * y + x - 1), 20 * (y - x**2)] - ) - - return energy, gradient - - optimiser = TestArmijoLineSearch( - energy_grad_func=energy_grad, init_step_size=0.1 - ) - optimiser.run(Molecule(name="blank"), method=Method()) - - assert optimiser.converged - assert optimiser._coords.e < optimiser._init_coords.e - - -@work_in_tmp_dir() -@requires_working_xtb_install -def test_xtb_h2_cart_opt(): - optimiser = ArmijoLineSearch() - assert not optimiser.converged - - h2 = Molecule(name="h2", atoms=[Atom("H"), Atom("H", x=1.5)]) - - # Should not converge in only two steps - ArmijoLineSearch.optimise(species=h2, method=XTB(), maxiter=2) - assert not optimiser.converged - - # Line search should step in the direction to reduce the distance - assert h2.distance(0, 1) < 1.4 diff --git a/tests/test_opt/test_opt.py b/tests/test_opt/test_opt.py index aa801e819..af216defd 100644 --- a/tests/test_opt/test_opt.py +++ b/tests/test_opt/test_opt.py @@ -1,11 +1,14 @@ import copy import os +import zipfile + import pytest import numpy as np from autode.methods import XTB from autode.calculations.types import CalculationType from autode.values import GradientRMS, PotentialEnergy +from autode.species.molecule import Molecule from autode.hessians import Hessian from autode.utils import work_in_tmp_dir from ..testutils import requires_working_xtb_install @@ -79,8 +82,8 @@ def test_coords_set(): def test_abs_diff_e(): # Define a intermediate optimiser state with two sets of coordinates optimiser = sample_cartesian_optimiser() - optimiser._history.append(CartesianCoordinates([0.0, 1.0])) - optimiser._history.append(CartesianCoordinates([0.0, 1.1])) + optimiser._history.add(CartesianCoordinates([0.0, 1.0])) + optimiser._history.add(CartesianCoordinates([0.0, 1.1])) # 2nd iteration for a history of two, indexed from 0 assert optimiser.iteration == 1 @@ -120,12 +123,12 @@ def test_optimiser_h_update(): c1 = CartesianCoordinates([1.0, 0.0, 0.0]) c1.h = np.eye(3) - optimiser._history.append(c1) + optimiser._history.add(c1) c2 = CartesianCoordinates([1.1, 0.0, 0.0]) c2.h = np.eye(3) - optimiser._history.append(c2) + optimiser._history.add(c2) # and try and update the (inverse) hessian, which is impossible without # an updater @@ -146,13 +149,6 @@ def test_history(): with pytest.raises(IndexError): _ = optimiser._history.penultimate - # or minimum in energy - with pytest.raises(IndexError): - _ = optimiser._history.minimum - - # and cannot contain a well in the energy - assert not optimiser._history.contains_energy_rise - @work_in_tmp_dir() @requires_working_xtb_install @@ -180,6 +176,11 @@ def test_xtb_h2_cart_opt(): # Should not converge in only two steps optimiser.run(method=XTB(), species=h2()) assert not optimiser.converged + # a trajectory file should be written + assert os.path.isfile("h2_opt_trj.zip") + # cleaning up optimiser will remove trajectory + optimiser.clean_up() + assert not os.path.isfile("h2_opt_trj.zip") @work_in_tmp_dir() @@ -229,6 +230,7 @@ def func(coords, m=None): optimiser.run(species=mol, method=Method()) +@work_in_tmp_dir() def test_last_energy_change_with_no_steps(): mol = h2() optimiser = HarmonicPotentialOptimiser( @@ -317,38 +319,208 @@ def test_multiple_optimiser_saves_overrides_not_append(): gtol=GradientRMS(0.01), etol=PotentialEnergy(1e-3), ) - optimiser.run(method=XTB(), species=h2()) - optimiser.save("tmp.traj") - - def n_lines_in_traj_file(): - return len(open("tmp.traj", "r").readlines()) + optimiser.run(method=XTB(), species=h2(), name="tmp.zip") - n_init_lines = n_lines_in_traj_file() - optimiser.save("tmp.traj") + assert os.path.isfile("tmp.zip") + with zipfile.ZipFile("tmp.zip") as file: + names = file.namelist() - assert n_lines_in_traj_file() == n_init_lines + old_n_coords = sum([1 for name in names if name.startswith("coords_")]) - -def test_optimiserhistory_operations_maintain_subclass(): optimiser = CartesianSDOptimiser( maxiter=2, - gtol=GradientRMS(0.2), - etol=PotentialEnergy(0.1), + gtol=GradientRMS(0.01), + etol=PotentialEnergy(1e-3), ) - coord = CartesianCoordinates(np.arange(6)) - optimiser._coords = coord - optimiser._coords = coord * 0.1 - hist = optimiser._history - assert len(hist) == 2 + optimiser.run(method=XTB(), species=h2(), name="tmp.zip") + # the file "tmp.zip" should be overwritten by new optimiser + with zipfile.ZipFile("tmp.zip") as file: + names = file.namelist() + + n_coords = sum([1 for name in names if name.startswith("coords_")]) + assert old_n_coords == n_coords + + +@work_in_tmp_dir() +def test_optimiser_print_geometries(caplog): + mol = Molecule(smiles="C=C", name="mymolecule") + coords1 = CartesianCoordinates(mol.coordinates) + opt = CartesianSDOptimiser(maxiter=20, gtol=1e-3, etol=1e-4) + opt._coords = coords1 + # cannot print geom without species + with pytest.raises(AssertionError): + opt.print_geometries() + + opt._species = mol + assert opt.iteration == 0 + with caplog.at_level("WARNING"): + opt.print_geometries() + assert "Optimiser did no steps, not saving .xyz" in caplog.text + assert not os.path.isfile("mymolecule_opt.trj.xyz") + opt._coords = coords1.copy() + opt.print_geometries() + assert os.path.isfile("mymolecule_opt.trj.xyz") + old_size = os.path.getsize("mymolecule_opt.trj.xyz") + # running should overwrite the geometries + opt.print_geometries() + new_size = os.path.getsize("mymolecule_opt.trj.xyz") + assert old_size == new_size + + +def _get_4_random_coordinates(): + coords_list = [] + for _ in range(4): + coords_list.append(CartesianCoordinates(np.random.rand(6))) + return coords_list - new_hist = copy.deepcopy(hist) - assert new_hist[0] is not hist[0] # must be deepcopy - add_hist = hist + new_hist - assert isinstance(add_hist, OptimiserHistory) +@work_in_tmp_dir() +def test_optimiser_history_storage(): + coords1, coords2, coords3, coords4 = _get_4_random_coordinates() - hist_slice = hist[:-1] - assert isinstance(hist_slice, OptimiserHistory) + hist = OptimiserHistory(maxlen=3) + # cannot close without opening a file + with pytest.raises(RuntimeError): + hist.close() + hist.open("test.zip") + assert os.path.isfile("test.zip") + # cannot reinitialise + with pytest.raises(RuntimeError, match="cannot initialise again"): + hist.open("test.zip") + # cannot add something that is not coordinates + with pytest.raises(ValueError, match="must be OptCoordinates"): + hist.add("x") + hist.add(coords1) + hist.add(coords2) + hist.add(coords3) + # nothing should be on disk yet + assert len(hist) == 3 and hist._n_stored == 0 + # now last coord is put on disk + hist.add(coords4) + assert len(hist) == 4 and hist._n_stored == 1 + assert len(hist._memory) == 3 + hist.close() + # now should be 3 more stored on disk + assert len(hist) == 4 and hist._n_stored == 4 + # adding new coords is forbidden + with pytest.raises(RuntimeError): + hist.add(coords1) + # iterate through the history in reverse + iterator = reversed(hist) + last = next(iterator) + before_last = next(iterator) + assert np.allclose(last, coords4) and np.allclose(before_last, coords3) + # clean up + hist.clean_up() + assert not os.path.isfile("test.zip") + + +@work_in_tmp_dir() +def test_optimiser_history_getitem(): + coords0, coords1, coords2, coords3 = _get_4_random_coordinates() + hist = OptimiserHistory(maxlen=2) + hist.open("test.zip") + hist.add(coords0) + hist.add(coords1) + hist.add(coords2) + hist.add(coords3) + assert np.allclose(hist[0], coords0) # from disk + hist[0].e = PotentialEnergy(0.001, "Ha") + assert hist[0].e is None # cannot modify disk + assert np.allclose(hist[2], coords2) # from memory + assert hist[2].e is None + hist[2].e = PotentialEnergy(0.01, "Ha") + assert np.isclose(hist[2].e, 0.01) + # slicing does not work + with pytest.raises(NotImplementedError): + _ = hist[0:1] + # can only have integer indices + with pytest.raises(ValueError): + _ = hist["x"] + with pytest.raises(IndexError): + _ = hist[4] + with pytest.raises(IndexError): + _ = hist[-5] + # if no disk backend, then old coordinates are lost + hist_nodisk = OptimiserHistory(maxlen=2) + hist_nodisk.add(coords0) + hist_nodisk.add(coords1) + hist_nodisk.add(coords2) + assert hist_nodisk[0] is None + assert hist_nodisk._n_stored == 0 + + +@work_in_tmp_dir() +def test_optimiser_history_reload(): + coords0, coords1, coords2, coords3 = _get_4_random_coordinates() + hist = OptimiserHistory(maxlen=2) + hist.open("savefile") + assert os.path.isfile("savefile.zip") # extension added + hist.add(coords0) + hist.add(coords1) + hist.add(coords2) + hist.add(coords3) + hist.close() + hist = None + with pytest.raises(FileNotFoundError, match="test.zip does not exist"): + _ = OptimiserHistory.load("test") + with open("test.zip", "w") as fh: + fh.write("abcd") + # error if file is not zip + with pytest.raises(ValueError, match="not a valid trajectory"): + _ = OptimiserHistory.load("test") + # error if file does not have the autodE opt header + with zipfile.ZipFile("new.zip", "w") as file: + fh = file.open("testfile", "w") + fh.write("abcd".encode()) + fh.close() + with pytest.raises(ValueError, match="not an autodE trajectory"): + _ = OptimiserHistory.load("new.zip") + hist = OptimiserHistory.load("savefile") + assert np.allclose(hist[-1], coords3) + assert np.allclose(hist[-2], coords2) + assert np.allclose(hist[-3], coords1) + + +@work_in_tmp_dir() +def test_optimiser_history_reload_works_with_one(): + coords0 = CartesianCoordinates(np.random.rand(6)) + hist = OptimiserHistory(maxlen=2) + + hist.open("savefile") + # adding None will not do anything + hist.add(None) + assert len(hist) == 0 + # just add one more coordinate + hist.add(coords0) + hist.close() + assert os.path.isfile("savefile.zip") + hist = OptimiserHistory.load("savefile") + assert len(hist) == 1 + assert np.allclose(hist[0], coords0) + assert hist[0] is hist[-1] + + +@work_in_tmp_dir() +def test_optimiser_history_save_load_params(): + hist = OptimiserHistory() + # cannot save or load without having file backing + with pytest.raises(RuntimeError, match="File not opened"): + hist.save_opt_params({"maxiter": 10}) + with pytest.raises(RuntimeError, match="File not opened"): + hist.get_opt_params() + hist.open("test.zip") + # cannot load as it is not available + with pytest.raises(FileNotFoundError, match="not found!"): + hist.get_opt_params() + # can now save and load + hist.save_opt_params({"maxiter": 10, "gtol": 1e-3}) + params = hist.get_opt_params() + assert len(params) == 2 + assert params["maxiter"] == 10 and params["gtol"] == 1e-3 + # cannot save again - overwrite not allowed + with pytest.raises(FileExistsError, match="already stored"): + hist.save_opt_params({"maxiter": 10}) def test_mocked_method(): @@ -360,7 +532,7 @@ def test_mocked_method(): def test_null_optimiser_methods(): optimiser = NullOptimiser() optimiser.run() - optimiser.save(filename="None") # run and saving does nothing + # run does nothing with pytest.raises(RuntimeError): _ = optimiser.final_coordinates diff --git a/tests/test_opt/test_opt_utils.py b/tests/test_opt/test_opt_utils.py new file mode 100644 index 000000000..c0507ef92 --- /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 diff --git a/tests/test_opt/test_rfo.py b/tests/test_opt/test_rfo.py index 8b31d1565..d7907a305 100644 --- a/tests/test_opt/test_rfo.py +++ b/tests/test_opt/test_rfo.py @@ -49,6 +49,7 @@ def _initialise_run(self) -> None: self._update_gradient_and_energy() +@work_in_tmp_dir() def test_simple_quadratic_opt(): optimiser = TestRFOOptimiser2D( e_func=lambda x, y: x**2 + y**2, @@ -62,6 +63,7 @@ def test_simple_quadratic_opt(): assert optimiser.iteration < 10 +@work_in_tmp_dir() def test_branin_opt(): def energy(x, y): return (y - 0.129 * x**2 + 1.6 * x - 6) ** 2 + 6.07 * np.cos(x) + 10 diff --git a/tests/test_opt/test_trm.py b/tests/test_opt/test_trm.py deleted file mode 100644 index c685ac396..000000000 --- a/tests/test_opt/test_trm.py +++ /dev/null @@ -1,258 +0,0 @@ -import os -import numpy as np -import pytest -from autode.species import Molecule -from autode.atoms import Atom -from autode.methods import XTB -from autode.values import Energy -from autode.utils import work_in_tmp_dir -from ..testutils import requires_working_xtb_install, work_in_zipped_dir -from autode.opt.coordinates import CartesianCoordinates, DIC -from autode.opt.optimisers.trm import ( - HybridTRMOptimiser, - CartesianHybridTRMOptimiser, -) - -here = os.path.dirname(os.path.abspath(__file__)) - - -@work_in_zipped_dir(os.path.join(here, "data", "opt.zip")) -def test_trm_step1(): - mol = Molecule("opt-test.xyz") - opt = HybridTRMOptimiser(maxiter=10, gtol=0.001, etol=0.001) - opt._species = mol - opt._build_internal_coordinates() - assert isinstance(opt._coords, DIC) - - grad = np.loadtxt("opt-test_grad.txt") - hess = np.loadtxt("opt-test_hess.txt") - opt._coords.update_g_from_cart_g(grad) - opt._coords.update_h_from_cart_h(hess) - - opt._step() - step = opt._history.final.to("cart") - opt._history.penultimate.to("cart") - step_size = np.linalg.norm(step) - assert np.isclose(step_size, opt.alpha, rtol=0.01) # 1% error margin - - opt = CartesianHybridTRMOptimiser(maxiter=10, gtol=0.001, etol=0.001) - opt._species = mol - opt._coords = CartesianCoordinates(mol.coordinates) - - opt._coords.update_g_from_cart_g(grad) - opt._coords.update_h_from_cart_h(hess) - - opt._step() - assert isinstance(opt._coords, CartesianCoordinates) - step = opt._history.final.to("cart") - opt._history.penultimate.to("cart") - step_size = np.linalg.norm(step) - assert np.isclose(step_size, opt.alpha, rtol=0.001) # 0.1% error margin - - -def test_damping_in_hybridtrm_optimiser(): - coord1 = CartesianCoordinates([1.0, -2.0, 1.0, 3.0, 1.1, 1.2]) - coord1.e = 0.10 - coord1.g = np.array([0.2, 0.3, 0.1, 0.2, 0.3, 0.4]) - coord1.h = np.eye(6) - coord2 = CartesianCoordinates([1.1, -1.9, 1.1, 3.1, 1.2, 1.3]) - coord2.e = 0.14 - coord2.g = np.array([0.1, 0.2, 0.01, 0.1, 0.2, 0.3]) - coord2.h = np.eye(6) - opt = CartesianHybridTRMOptimiser(maxiter=2, gtol=0.1, etol=0.1) - - # simulate an oscillation happening - opt._coords = coord1 - opt._coords = coord2 - opt._coords = coord1 - opt._coords = coord2 - opt._coords = coord1 - opt._coords = coord2 - - assert opt._is_oscillating() # should register as oscillation - opt._step() # should take damped step - # 50% mixing - avg_coord = (coord1 + coord2) / 2.0 - assert np.allclose(avg_coord, opt._coords) - - -def test_trim_convergence_gtol_overachieved(): - # if rms grad is <= gtol/10, assumed to be converged, even if - # delta E criteria is not met - opt = CartesianHybridTRMOptimiser(maxiter=3, gtol=0.1, etol=0.01) - coord1 = CartesianCoordinates([1.0, -2.0, 1.0, 3.0, 1.1, 1.2]) - coord1.e = Energy(0.10, "Ha") - coord1.g = np.array([0.01, 0.03, 0.01, 0.02, 0.03, 0.04]) - coord2 = coord1.copy() - coord2.e = Energy(0.05, "Ha") - coord2.g = coord1.g / 10 - opt._coords = coord1 - opt._coords = coord2 - - # energy criteria not achieved - assert opt._abs_delta_e > opt.etol - # grad criteria overachieved - assert opt._g_norm <= opt.gtol / 10 - assert opt.converged - - -@work_in_tmp_dir() -@requires_working_xtb_install -def test_trm_molecular_opt(): - mol = Molecule(smiles="O") - assert [atom.label for atom in mol.atoms] == ["O", "H", "H"] - - HybridTRMOptimiser.optimise(mol, method=XTB(), gtol=1.0e-4) - - # Check optimised distances are similar to running the optimiser in XTB - for oh_atom_idx_pair in [(0, 1), (0, 2)]: - assert np.isclose( - mol.distance(*oh_atom_idx_pair).to("Å"), 0.9595, atol=1e-3 - ) - - assert np.isclose(mol.distance(1, 2), 1.5438, atol=1e-3) - - mol = Molecule(smiles="O") - CartesianHybridTRMOptimiser.optimise(mol, method=XTB(), gtol=1e-4) - - # Same check for Cartesian - for oh_atom_idx_pair in [(0, 1), (0, 2)]: - assert np.isclose( - mol.distance(*oh_atom_idx_pair).to("Å"), 0.9595, atol=1e-3 - ) - - assert np.isclose(mol.distance(1, 2), 1.5438, atol=1e-3) - - -@work_in_tmp_dir() -@requires_working_xtb_install -def test_trust_update(caplog): - init_trust = 0.05 - water_atoms = [ - Atom("O", -0.0011, 0.3631, -0.0000), - Atom("H", -0.8250, -0.1819, -0.0000), - Atom("H", 0.8261, -0.1812, 0.0000), - ] - water = Molecule(atoms=water_atoms) - opt = HybridTRMOptimiser( - maxiter=10, gtol=1.0e-3, etol=1.0e-4, init_trust=init_trust - ) - opt._species = water.copy() # take copy so original is not modified - opt._method = XTB() - opt._n_cores = 1 - opt._initialise_run() - # store last grad - last_g = opt._coords.g.reshape(-1, 1).copy() - last_h = opt._coords.h.copy() - - opt._step() - opt._update_gradient_and_energy() - last_step = opt._history[-1].raw - opt._history[-2].raw - last_cart_step = opt._history[-1].to("cart") - opt._history[-2].to("cart") - cart_step_size = np.linalg.norm(last_cart_step) - pred_delta_e = float(last_g.T @ last_step) - pred_delta_e += 0.5 * (last_step.T @ last_h @ last_step) - # pred_dE should be around -0.0025446056 Ha (may change with different version of xTB) - - def simulate_energy_change_ratio_update_trust(ratio): - opt.alpha = init_trust - opt._history.final.e = ( - opt._history.penultimate.e + ratio * pred_delta_e - ) - opt._update_trust_radius() - - # if update is not allowed, there shouldn't be update - opt._upd_alpha = False - simulate_energy_change_ratio_update_trust(0.2) - assert np.isclose(opt.alpha, init_trust) - opt._upd_alpha = True - - simulate_energy_change_ratio_update_trust(0.2) - assert np.isclose(opt.alpha, 0.5 * min(init_trust, cart_step_size)) - - simulate_energy_change_ratio_update_trust(0.5) - assert np.isclose(opt.alpha, init_trust) - - simulate_energy_change_ratio_update_trust(1.0) - assert abs(cart_step_size - init_trust) / init_trust < 0.05 - assert np.isclose(opt.alpha, 1.414 * init_trust) - - simulate_energy_change_ratio_update_trust(1.3) - assert np.isclose(opt.alpha, init_trust) - - simulate_energy_change_ratio_update_trust(1.7) - assert np.isclose(opt.alpha, 0.5 * min(init_trust, cart_step_size)) - - # if energy change too high > 2.0 or too low < -1.0, trust radius - # is decreased, and last step is rejected - with caplog.at_level("WARNING"): - simulate_energy_change_ratio_update_trust(2.2) - assert "rejecting last geometry step" in caplog.text - assert np.isclose(opt.alpha, 0.5 * min(init_trust, cart_step_size)) - # should remove last coords, bringing back original state - assert opt.iteration == 0 - assert np.allclose( - opt._history.final.to("cart"), water.coordinates.flatten() - ) - - -@work_in_tmp_dir() -@requires_working_xtb_install -def test_optimiser_plotting(): - mol = Molecule(smiles="O") - - opt = HybridTRMOptimiser(maxiter=100, gtol=1e-3, etol=1e-3) - opt.run(mol, method=XTB()) - - opt.plot_optimisation() - assert os.path.isfile(f"{mol.name}_opt_plot.pdf") - - -@work_in_tmp_dir() -def test_optimiser_xyz_trajectory(): - mol = Molecule(smiles="N#N") - opt = CartesianHybridTRMOptimiser(maxiter=10, gtol=1e-3, etol=1e-3) - coord1 = CartesianCoordinates(mol.coordinates) - opt._coords = coord1 - - with pytest.raises(AssertionError): - # no species should raise error - opt.print_geometries() - - opt._species = mol - assert opt.iteration == 0 - # xyz trajectory only if more than 1 coordinates present - opt.print_geometries() - assert not os.path.isfile(f"{mol.name}_opt.trj.xyz") - - coord2 = coord1.copy() * 0.9 - opt._coords = coord2 - assert opt.iteration == 1 - opt.print_geometries() - assert os.path.isfile(f"{mol.name}_opt.trj.xyz") - - -@work_in_tmp_dir() -def test_optimiser_plotting_sanity_checks(caplog): - mol = Molecule(smiles="N#N") - - opt = CartesianHybridTRMOptimiser(maxiter=10, gtol=1e-3, etol=1e-3) - coord1 = CartesianCoordinates(mol.coordinates) - coord1.e = Energy(0.1, "Ha") - coord1.update_g_from_cart_g(np.array([0.01, 0.02, 0.05, 0.06, 0.03, 0.07])) - opt._coords = coord1 - opt._species = mol - assert opt.iteration == 0 - assert not opt.converged - # plotting does not work if less than 2 points - with caplog.at_level("WARNING"): - opt.plot_optimisation(filename="test_plot0.pdf") - assert not os.path.isfile("test_plot0.pdf") - assert "Less than 2 points, cannot draw optimisation" in caplog.text - - opt._coords = coord1.copy() - # either rms_grad or energy plot has to be requested - with caplog.at_level("WARNING"): - opt.plot_optimisation("test_plot1.pdf", False, False) - assert not os.path.isfile("test_plot1.pdf") - assert not opt.converged - assert "Must plot either energies or RMS gradients" in caplog.text - assert "Optimisation is not converged, drawing a plot" in caplog.text diff --git a/tests/test_opt/test_trust_region_methods.py b/tests/test_opt/test_trust_region_methods.py deleted file mode 100644 index 704e9dfcd..000000000 --- a/tests/test_opt/test_trust_region_methods.py +++ /dev/null @@ -1,225 +0,0 @@ -import numpy as np -import pytest - -from autode.species import Molecule -from autode.opt.coordinates import CartesianCoordinates -from autode.opt.optimisers.trust_region import ( - CauchyTROptimiser, - DoglegTROptimiser, - CGSteihaugTROptimiser, -) -from .setup import Method - - -def branin_energy(x, y): - return (y - 0.129 * x**2 + 1.6 * x - 6) ** 2 + 6.07 * np.cos(x) + 10 - - -class BraninCauchyTROptimiser(CauchyTROptimiser): - __test__ = False - - def _space_has_degrees_of_freedom(self) -> bool: - return True - - def _update_gradient_and_energy(self) -> None: - """Update the gradient and energy for the Branin function - - f(x, y) = (y - 0.129x^2 + 1.6x - 6)^2 + 6.07cos(x) + 10 - """ - x, y = self._coords - - self._coords.e = branin_energy(x, y) - - grad = [ - ( - 2 * (1.6 - 0.258 * x) * (y - 0.129 * x**2 + 1.6 * x - 6) - - 6.07 * np.sin(x) - ), - (2 * (y - 0.129 * x**2 + 1.6 * x - 6)), - ] - - self._coords.g = np.array(grad) - - h_xx = ( - 2 * (1.6 - 0.258 * x) ** 2 - - 0.516 * (-0.129 * x**2 + 1.6 * x + y - 6) - - 6.07 * np.cos(x) - ) - h_xy = 2 * (1.6 - 0.258 * x) - h_yy = 2 - self._coords.h = np.array([[h_xx, h_xy], [h_xy, h_yy]]) - - def _log_convergence(self) -> None: - attrs = [ - self.iteration, - self._coords.e, - *self._coords, - self.rho, - self.alpha, - np.linalg.norm(self.p), - self._g_norm, - ] - - if hasattr(self, "tau"): - attrs.append(self.tau) - - for thing in attrs: - print(f"{round(thing, 3):10.3f}" f"", end=" ") - print() - - @property - def converged(self) -> bool: - return np.linalg.norm(self._coords.g) < self._gtol - - -class BraninDoglegTROptimiser(DoglegTROptimiser): - def _space_has_degrees_of_freedom(self) -> bool: - return True - - def _update_gradient_and_energy(self) -> None: - return BraninCauchyTROptimiser._update_gradient_and_energy(self) - - def _log_convergence(self) -> None: - return BraninCauchyTROptimiser._log_convergence(self) - - @property - def converged(self) -> bool: - return np.linalg.norm(self._coords.g) < self._gtol - - -class BraninCGSteihaugTROptimiser(CGSteihaugTROptimiser): - def _space_has_degrees_of_freedom(self) -> bool: - return True - - def _update_gradient_and_energy(self) -> None: - return BraninCauchyTROptimiser._update_gradient_and_energy(self) - - def _log_convergence(self) -> None: - return BraninCauchyTROptimiser._log_convergence(self) - - @property - def converged(self) -> bool: - return ( - self._coords.g is not None - and np.linalg.norm(self._coords.g) < self._gtol - ) - - -def test_trm_base_properties(): - init_coords = CartesianCoordinates([6.0, 14.0]) - - optimiser = BraninCauchyTROptimiser( - maxiter=20, - etol=100, # Some large value - trust_radius=2.0, - coords=init_coords, - gtol=0.01, - ) - - # Updating the Hessian should be possible, and be the identity matrix - CauchyTROptimiser._update_hessian(optimiser) - assert np.allclose(optimiser._coords.h, np.eye(len(init_coords))) - - # rho requires a gradient to be evaluated - optimiser._history.append(CartesianCoordinates([5.5, 13.0])) - with pytest.raises(Exception): - _ = optimiser.rho - - -def test_branin_minimisation(): - """Uses the example from: - https://optimization.mccormick.northwestern.edu/index.php/Trust-region_methods - """ - - init_coords = CartesianCoordinates([6.0, 14.0]) - - optimiser = BraninCauchyTROptimiser( - maxiter=20, - etol=100, # Some large value - trust_radius=2.0, - coords=init_coords, - gtol=0.01, - max_trust_radius=5.0, - t_1=0.25, - t_2=2.0, - eta_1=0.2, - eta_2=0.25, - eta_3=0.75, - ) - - optimiser.run(Molecule(name="blank"), method=Method()) - assert optimiser.converged - assert np.allclose(optimiser._coords, np.array([3.138, 2.252]), atol=0.01) - - -def test_branin_dogleg_minimisation(): - # Dogleg optimiser doesn't seem to be so efficient - TODO: Check - optimiser = BraninDoglegTROptimiser( - maxiter=1000, - etol=100, # Some large value - trust_radius=2.0, - coords=CartesianCoordinates([6.0, 14.0]), - gtol=0.01, - max_trust_radius=5.0, - t_1=0.25, - t_2=2.0, - eta_1=0.2, - eta_2=0.25, - eta_3=0.75, - ) - - optimiser.run(Molecule(name="blank"), method=Method()) - - assert optimiser.converged - assert np.allclose(optimiser._coords, np.array([3.138, 2.252]), atol=0.02) - - -def test_branin_cg_minimisation(): - optimiser = BraninCGSteihaugTROptimiser( - maxiter=1000, - etol=100, # Some large value - trust_radius=2.0, - coords=CartesianCoordinates([6.0, 14.0]), - gtol=0.01, - max_trust_radius=5.0, - t_1=0.25, - t_2=2.0, - eta_1=0.2, - eta_2=0.25, - eta_3=0.75, - ) - - optimiser.run(Molecule(name="blank"), method=Method()) - - assert optimiser.converged - assert np.allclose(optimiser._coords, np.array([3.138, 2.252]), atol=0.02) - - # Should also be able to optimise directly - coords = CartesianCoordinates([6.0, 14.0]) - BraninCGSteihaugTROptimiser.optimise( - Molecule(name="blank"), - method=Method(), - gtol=0.01, - etol=10, - coords=coords, - ) - - assert optimiser.converged - - -def test_base_cg_properties(): - optimiser = CGSteihaugTROptimiser( - maxiter=10, - trust_radius=1.0, - etol=1, - gtol=0.1, - coords=CartesianCoordinates([0.1, 0.0]), - ) - - assert not optimiser.converged - - optimiser._coords.g = np.zeros(shape=(2,)) - optimiser._solve_subproblem() - - # For an almost zero gradient the step size should be zero - assert np.linalg.norm(optimiser.p) < 1e-10 diff --git a/tests/test_opt/test_trust_region_molecule.py b/tests/test_opt/test_trust_region_molecule.py deleted file mode 100644 index 3aa8fc973..000000000 --- a/tests/test_opt/test_trust_region_molecule.py +++ /dev/null @@ -1,22 +0,0 @@ -from autode.species.molecule import Molecule -from autode.opt.optimisers.trust_region import CGSteihaugTROptimiser -from autode.opt.coordinates.cartesian import CartesianCoordinates -from autode.methods import XTB -from ..testutils import requires_working_xtb_install - - -@requires_working_xtb_install -def test_h2o_opt(): - mol = Molecule(smiles="O") - - optimiser = CGSteihaugTROptimiser( - maxiter=20, - trust_radius=0.02, - etol=1e-6, - gtol=1e-3, - coords=CartesianCoordinates(mol.coordinates), - ) - optimiser.run(mol, method=XTB()) - - assert optimiser.converged - assert optimiser.iteration < 20 diff --git a/tests/test_plotting.py b/tests/test_plotting.py index 03a96c646..e5a0a999b 100644 --- a/tests/test_plotting.py +++ b/tests/test_plotting.py @@ -230,8 +230,8 @@ def test_optimiser_plot(): x.e = -1.45 x.g = np.array([1.0, 0.9, 1.2, 0.4, 3.4, 0.3]) x2 = x.copy() - hist.append(x) - hist.append(x2) + hist.add(x) + hist.add(x2) plotting.plot_optimiser_profile( history=hist, plot_energy=True, diff --git a/tests/test_reaction_class.py b/tests/test_reaction_class.py index b942e4490..a7b5e31ca 100644 --- a/tests/test_reaction_class.py +++ b/tests/test_reaction_class.py @@ -481,6 +481,36 @@ def test_barrierless_h_g(): assert np.isclose(rxn.delta("G‡"), 0.7) # -2+0.6 -> -1+0.3 --> ∆ = 0.7 +@pytest.mark.parametrize("energy_syn", ["E", "Energy"]) +@pytest.mark.parametrize("enthalpy_syn", ["H", "Enthalpy"]) +@pytest.mark.parametrize("free_syn", ["G", "free energy", "free_energy"]) +@pytest.mark.parametrize("ts_syn", ["ddagger", "‡", "double dagger"]) +def test_energy_synonyms(energy_syn, enthalpy_syn, free_syn, ts_syn): + a = Reactant(atoms=[Atom("H"), Atom("H", x=0.7, y=0.7), Atom("H", x=1.0)]) + a.energies.extend( + [PotentialEnergy(-2), EnthalpyCont(0.2), FreeEnergyCont(0.6)] + ) + + b = Product(atoms=[Atom("H"), Atom("H", x=-1.0), Atom("H", x=1.0)]) + b.energies.extend( + [PotentialEnergy(-1), EnthalpyCont(0.1), FreeEnergyCont(0.3)] + ) + + rxn = reaction.Reaction(a, b) + # Test potential energies + assert np.isclose(rxn.delta(energy_syn + ts_syn), 1.0) + + # Test enthalpies + assert np.isclose( + rxn.delta(enthalpy_syn + ts_syn), 0.9 + ) # -2+0.2 -> -1+0.1 --> ∆ = 0.9 + + # Test free energies + assert np.isclose( + rxn.delta(free_syn + ts_syn), 0.7 + ) # -2+0.6 -> -1+0.3 --> ∆ = 0.7 + + def test_same_composition(): r1 = reaction.Reaction( Reactant(atoms=[Atom("C"), Atom("H", x=1)]), diff --git a/tests/test_species.py b/tests/test_species.py index 51cc6b0b3..c1e6df247 100644 --- a/tests/test_species.py +++ b/tests/test_species.py @@ -253,6 +253,11 @@ def test_set_gradients(): test_mol.gradient = np.zeros(shape=(2, 3)) assert test_mol.gradient.units == ha_per_ang + # will reshape numpy array if possible + arr = np.random.rand(6) + test_mol.gradient = arr + assert np.allclose(test_mol.gradient, arr.reshape(-1, 3)) + def test_species_solvent(): assert mol.solvent is None diff --git a/tests/test_ts/test_mode_checking.py b/tests/test_ts/test_mode_checking.py index 6c0acf6df..d4ab23012 100644 --- a/tests/test_ts/test_mode_checking.py +++ b/tests/test_ts/test_mode_checking.py @@ -1,3 +1,7 @@ +import os +import numpy as np + +from autode.species.molecule import Molecule from autode.calculations import Calculation from autode.transition_states.base import TSbase from autode.transition_states.base import imag_mode_generates_other_bonds @@ -8,7 +12,6 @@ from autode.bond_rearrangement import BondRearrangement from autode.methods import ORCA from .. import testutils -import os here = os.path.dirname(os.path.abspath(__file__)) orca = ORCA() @@ -70,6 +73,18 @@ def test_graph_no_other_bonds(): ) +def test_disp_molecule_has_same_solvent(): + mol = Molecule(smiles="[H][H]", solvent_name="water") + mol.hessian = np.eye(3 * mol.n_atoms, 3 * mol.n_atoms) + + disp_mol = displaced_species_along_mode( + species=mol, + mode_number=1, + ) + assert disp_mol.solvent is not None + assert disp_mol.solvent == mol.solvent + + def has_correct_mode(name, fbonds, bbonds): calc = Calculation( name=name, diff --git a/tests/test_wrappers/data/nwchem.zip b/tests/test_wrappers/data/nwchem.zip index 4691c8ae2..a10a54a11 100644 Binary files a/tests/test_wrappers/data/nwchem.zip and b/tests/test_wrappers/data/nwchem.zip differ