Source code for pysisyphus.optimizers.RFOptimizer

# [1] Optimization review
# [2] Trust region method
# [3] 10.1007/978-0-387-40065-5 Numerical optimization
# [4] 10.1007/s00214-016-1847-3 Explorations of some refinements

import numpy as np

from pysisyphus.Geometry import Geometry
from pysisyphus.helpers_pure import rms
from pysisyphus.optimizers.HessianOptimizer import HessianOptimizer
from pysisyphus.optimizers.poly_fit import poly_line_search
from pysisyphus.optimizers.gdiis import gdiis, gediis

[docs] class RFOptimizer(HessianOptimizer):
[docs] def __init__( self, geometry: Geometry, line_search: bool = True, gediis: bool = False, gdiis: bool = True, gdiis_thresh: float = 2.5e-3, gediis_thresh: float = 1e-2, gdiis_test_direction: bool = True, max_micro_cycles: int = 25, adapt_step_func: bool = False, **kwargs, ) -> None: """ Rational function Optimizer. Parameters ---------- geometry Geometry to be optimized. line_search Whether to carry out implicit line searches. gediis Whether to enable GEDIIS. gdiis Whether to enable GDIIS. gdiis_thresh Threshold for rms(forces) to enable GDIIS. gediis_thresh Threshold for rms(step) to enable GEDIIS. gdiis_test_direction Whether to the overlap of the RFO step and the GDIIS step. max_micro_cycles Number of restricted-step microcycles. Disabled by default. adapt_step_func Whether to switch between shifted Newton and RFO-steps. Other Parameters ---------------- **kwargs Keyword arguments passed to the Optimizer/HessianOptimizer baseclass. """ super().__init__(geometry, max_micro_cycles=max_micro_cycles, **kwargs) self.line_search = line_search self.gediis = gediis self.gdiis = gdiis self.gdiis_thresh = gdiis_thresh # Will be compared to rms(step) self.gediis_thresh = gediis_thresh # Will be compared to rms(forces) self.gdiis_test_direction = gdiis_test_direction self.adapt_step_func = adapt_step_func self.successful_gediis = 0 self.successful_gdiis = 0 self.successful_line_search = 0
[docs] def optimize(self): energy, gradient, H, big_eigvals, big_eigvecs, resetted = self.housekeeping() step_func, pred_func = self.get_step_func(big_eigvals, gradient) ref_gradient = gradient.copy() # Reference step, used for judging the proposed GDIIS step ref_step = step_func(big_eigvals, big_eigvecs, gradient) # Right everything is in place to check for convergence. If all values are below # the thresholds, there is no need to do additional inter/extrapolations. if self.check_convergence(ref_step)[0]: # Drop conv_info self.log("Convergence achieved! Skipping inter/extrapolation.") return ref_step # Try to interpolate an intermediate geometry, either from GDIIS or line search. # # Set some defaults ip_gradient = None ip_step = None diis_result = None # Check if we can do GDIIS or GEDIIS. If we (can) do a line search is decided # after trying GDIIS. rms_forces = rms(gradient) rms_step = rms(ref_step) can_diis = (rms_step <= self.gdiis_thresh) and (not resetted) can_gediis = (rms_forces <= self.gediis_thresh) and (not resetted) # GDIIS / GEDIIS, prefer GDIIS over GEDIIS if self.gdiis and can_diis: # Gradients as error vectors err_vecs = -np.array(self.forces) diis_result = gdiis( err_vecs, self.coords, self.forces, ref_step, test_direction=self.gdiis_test_direction, ) self.successful_gdiis += 1 if diis_result else 0 # Don't try GEDIIS if GDIIS failed. If GEDIIS should be tried after GDIIS failed # comment the line below and uncomment the line following it. elif self.gediis and can_gediis: # if self.gediis and can_gediis and (diis_result == None): diis_result = gediis(self.coords, self.energies, self.forces, hessian=H) self.successful_gediis += 1 if diis_result else 0 try: ip_coords = diis_result.coords ip_step = ip_coords - self.geometry.coords ip_gradient = -diis_result.forces # When diis_result is None except AttributeError: self.log("GDIIS didn't succeed.") # Try line search if GDIIS failed or not requested if self.line_search and (diis_result is None) and (not resetted): ip_energy, ip_gradient, ip_step = poly_line_search( energy, self.energies[-2], gradient, -self.forces[-2], self.steps[-1], cubic_max_x=-1, quartic_max_x=2, logger=self.logger, ) self.successful_line_search += 1 if ip_gradient is not None else 0 # Use the interpolated gradient for the RFO step if interpolation succeeded if (ip_gradient is not None) and (ip_step is not None): gradient = ip_gradient # Keep the original gradient when the interpolation failed, but recreate # ip_step, as it will be returned as None from poly_line_search(). else: ip_step = np.zeros_like(gradient) step = step_func(big_eigvals, big_eigvecs, gradient) # Form full step. If we did not interpolate or interpolation failed, # ip_step will be zero. step = step + ip_step # Use the original, actually calculated, gradient prediction = pred_func(ref_gradient, H, step) self.predicted_energy_changes.append(prediction) return step
[docs] def postprocess_opt(self): msg = ( f"Successful invocations:\n" f"\t GEDIIS: {self.successful_gediis}\n" f"\t GDIIS: {self.successful_gdiis}\n" f"\tLine Search: {self.successful_line_search}\n" ) self.log(msg)