7. Minimization

Searching for minimum energy geometries of molecules is preferably done using second order methods that employ hessian information. Using the plain hessian \(H\) for step prediction through \(p=-H^{-1}g\), with \(g\) being the gradient vector, may yield erroneous uphill steps when \(H\) has negative eigenvalues.

This can be understood by calculating the step in the basis of the hessian eigenvectors \(V\).

\[\begin{split}\begin{align} \tilde{H} &= V^T H V \\ \tilde{g} &= V^T g \\ \tilde{p_i} &= -\frac{\tilde{g}_i}{\tilde{H}_{ii}} \\ \end{align}\end{split}\]

\(\tilde{H}, \tilde{g}\) and \(\tilde{p}\) are transformed into the eigenvector basis and the subscript \(i\) indicates the component belongs to the \(i\)-th eigenvalue (-vector) of \(H\). As the gradient always points into the direction of greater function values dividing it by a negative eigenvalues \(\tilde{H}_{ii}\) will lead to a step in uphill direction along the \(i\)-th eigenvector. The step in the original basis is obtained by a simple back-transformation:

\[p = V \tilde{p}\]

A step in downhill direction can be ensured by introducing an appropriate shift parameter \(\lambda\) that must be smaller than the smallest eigenvalue of \(H\):

\[\begin{split}\tilde{p_i} = -\frac{\tilde{g}_i}{\tilde{H}_{ii} - \lambda} \\\end{split}\]

In pysisyphus the shift parameter \(\lambda\) is obtained by the Rational Function Optimization approach.

When the root-mean-square of the predicted step \(p\) drops below a certain threshold (default is 0.0025 au or rad) controlled GDIIS is tried. If GDIIS fails or is not yet possible a polynomial line-search is conducted. Using the latest two energies and the projection of the latest two gradients onto the last step the coefficients of a cubic polynomial are defined unambiguously. By requiring \(\mathrm{d}^2 f(x)/\mathrm{d}x^2 >= 0\) and the equality holding at exactly one point also a quartic polynomial can be determined.

For now line-searches using a cubic polynomial are disabled as they seem to degrade the optimizer performance, so only the constrained quartic polynomial is tried. By also using projected hessian information the coefficients of a quintic polynomial can be determined, but this also seems to be inferior compared to the constrained quartic polynomial.

Convergence is indicated when the root-mean-square and the maximum value of the gradient and step drop below a certain threshold. It is not uncommon that the gradient convergence is achieved before step convergence. By using overachieve_factor: [n, float > 1] in the YAML input convergence will be signalled, when gradient convergence is overachieved by factor n.

7.1. YAML Example

Below you can find an example YAML-input including the most important options that the user may want to modify when using the RFOptimizer.

opt:
 type: rfo                      # Optimization algorithm
 max_cycles: 50                 # Maximum number of optimization cycles

 overachieve_factor: 2          # Indicate convergence, regardless of the
                                # proposed step when max(grad) and rms(grad)
                                # are overachieved by factor [n]

 do_hess: True                  # Calculate the hessian at the final geometry
                                # after the optimization.

 #hessian_recalc: None          # Recalculate exact hessian every n-th cylce

 #hessian_recalc_adapt: None    # Expects a float. Recalculate exact hessian
                                # whenever the gradient norms drops below
                                # 1/[n] of the gradient norm at the last hessian
                                # recalculation.

 #hessian_init: fischer         # Type of model hessian. Other options are: 'calc,
                                # simple, xtb, lindh, swart, unit'

 #hessian_update: bfgs          # Hessian-update. Other options are: 'flowchart,
                                # damped_bfgs, bofill'. bofill is not recommended
                                # for minimum searches.
 #small_eigval_thresh: 1e-8     # Neglect eigenvalues and corresponding eigenvectors
                                # below this threshold.

 #max_micro_cycles: 50          # No. of micro cycles for the RS-variants. Does not apply
                                # to TRIM.

 #trust_radius: 0.3             # Initial trust radius.
 #trust_max: 1.0                # Max. trust radius
 #trust_min: 0.1                # Min. trust radius

 #line_search: True             # Do line search

 #gdiis_thresh: 0.0025          # May do GDIIS if rms(step) falls below this threshold
 #gediis_thresh: 0.01           # May do GEDIIS if rms(grad) falls below this threshold
 #gdiis: True                   # Do controlled GDIIS after 'gdiis_thresh' is reached
 #gediis: False                 # Do GEDIIS after 'gediis_thresh' is reached

calc:
 type: turbomole
 control_path: control_path_pbe0_def2svp_s1     # Path to the prepared calculation
 track: True                                    # Activate excited state tracking
 ovlp_type: tden                                # Track with transition density matrix overlaps
 charge: 0
 mult: 1
 pal: 4
 mem: 2000

geom:
 type: redund
 fn: cytosin.xyz

Further examples for optimizations from .yaml input can be found here.

7.2. Convergence Thresholds

Optimization convergence is indicated when four critera are met. Each cycle the root-mean-square and the absolute maximum of force and proposed step vectors are checked. The default criteria in pysisyphus correspond to the opt=loose in Gaussian. Different thresholds can easily be requested by supplying the respective keyword in the YAML input.

Convergence is also indicated, when max(force) and rms(force) are overachieved, by a certain factor (overachieve_factor, default=3), similar to ORCAs optimizer. Further keywords, controlling the convergence check are found below.

opt:
 ...
 thresh: gau_loose           # See table below for different thresholds.
 overachieve_factor: 3
 rms_force: null             # Derive four critera from this value
 rms_force_only: False       # Only check rms(force)
 max_force_only: False       # Only check max(force)
 # Chain-of-States specific
 coord_diff_thresh: 0.01     # Terminate on insignificant coordinate change
 reparam_thresh: 0.001       # Terminate on insignificant coordinate change upon reparametrization
 ...

Threshold	max(force)	rms(forces)	max(step)	rms(step)	Comment
gau_loose	2.5e-3	1.7e-3	1.0e-2	6.7e-3	default
gau	4.5e-4	3.0e-4	1.8e-3	1.2e-3
gau_tight	1.5e-5	1.0e-5	6.0e-5	4.0e-5
gau_vtight	2.0e-6	1.0e-6	6.0e-6	4.0e-6
baker	3.0e-4	2.0e-4	3.0e-4	2.0e-4	energy difference and step are also checked
never	2.0e-6	1.0e-6	6.0e-6	4.0e-6	dummy thresholds; never converges

7.3. Available Optimizers

Base class for optimizers that utilize Hessian information.

class pysisyphus.optimizers.Optimizer.ConvInfo(cur_cycle, energy_converged, max_force_converged, rms_force_converged, max_step_converged, rms_step_converged, desired_eigval_structure)

Bases: object

cur_cycle: int

desired_eigval_structure: bool

energy_converged: bool

get_convergence()

is_converged()

max_force_converged: bool

max_step_converged: bool

rms_force_converged: bool

rms_step_converged: bool

class pysisyphus.optimizers.Optimizer.Optimizer(geometry, thresh='gau_loose', max_step=0.04, max_cycles=150, min_step_norm=1e-08, assert_min_step=True, rms_force=None, rms_force_only=False, max_force_only=False, force_only=False, converge_to_geom_rms_thresh=0.05, align=False, align_factor=1.0, dump=False, dump_restart=False, print_every=1, prefix='', reparam_thresh=0.001, reparam_check_rms=True, reparam_when='after', overachieve_factor=0.0, check_eigval_structure=False, restart_info=None, check_coord_diffs=True, coord_diff_thresh=0.01, fragments=None, monitor_frag_dists=0, out_dir='.', h5_fn='optimization.h5', h5_group_name='opt')

Bases: object

__init__(geometry, thresh='gau_loose', max_step=0.04, max_cycles=150, min_step_norm=1e-08, assert_min_step=True, rms_force=None, rms_force_only=False, max_force_only=False, force_only=False, converge_to_geom_rms_thresh=0.05, align=False, align_factor=1.0, dump=False, dump_restart=False, print_every=1, prefix='', reparam_thresh=0.001, reparam_check_rms=True, reparam_when='after', overachieve_factor=0.0, check_eigval_structure=False, restart_info=None, check_coord_diffs=True, coord_diff_thresh=0.01, fragments=None, monitor_frag_dists=0, out_dir='.', h5_fn='optimization.h5', h5_group_name='opt')

Optimizer baseclass. Meant to be subclassed.

Parameters:

• geometry (Geometry) -- Geometry to be optimized.

• thresh (Literal['gau_loose', 'gau', 'gau_tight', 'gau_vtight', 'baker', 'never'], default: 'gau_loose') -- Convergence threshold.

• max_step (float, default: 0.04) -- Maximum absolute component of the allowed step vector. Utilized in optimizers that don't support a trust region or line search.

• max_cycles (int, default: 150) -- Maximum number of allowed optimization cycles.

• min_step_norm (float, default: 1e-08) -- Minimum norm of an allowed step. If the step norm drops below this value a ZeroStepLength-exception is raised. The unit depends on the coordinate system of the supplied geometry.

• assert_min_step (bool, default: True) -- Flag that controls whether the norm of the proposed step is check for being too small.

• rms_force (Optional[float], default: None) -- Root-mean-square of the force from which user-defined thresholds are derived. When 'rms_force' is given 'thresh' is ignored.

• rms_force_only (bool, default: False) -- When set, convergence is signalled only based on rms(forces).

• max_force_only (bool, default: False) -- When set, convergence is signalled only based on max(|forces|).

• force_only (bool, default: False) -- When set, convergence is signalled only based on max(|forces|) and rms(forces).

• converge_to_geom_rms_thresh (float, default: 0.05) -- Threshold for the RMSD with another geometry. When the RMSD drops below this threshold convergence is signalled. Only used with Growing Newton trajectories.

• align (bool, default: False) -- Flag that controls whether the geometry is aligned in every step onto the coordinates of the previous step. Must not be used with internal coordinates.

• align_factor (float, default: 1.0) -- Factor that controls the strength of the alignment. 1.0 means full alignment, 0.0 means no alignment. The factor mixes the rotation matrix of the alignment with the identity matrix.

• dump (bool, default: False) -- Flag to control dumping/writing of optimization progress to the filesystem

• dump_restart (bool, default: False) -- Flag to control whether restart information is dumped to the filesystem.

• print_every (int, default: 1) -- Report optimization progress every nth cycle.

• prefix (str, default: '') -- Short string that is prepended to several files created by the optimizer. Allows distinguishing several optimizations carried out in the same directory.

• reparam_thresh (float, default: 0.001) -- Controls the minimal allowed similarity between coordinates after two successive reparametrizations. Convergence is signalled if the coordinates did not change significantly.

• reparam_check_rms (bool, default: True) -- Whether to check for (too) similar coordinates after reparametrization.

• reparam_when (Optional[Literal['before', 'after']], default: 'after') -- Reparametrize before or after calculating the step. Can also be turned off by setting it to None.

• overachieve_factor (float, default: 0.0) -- Signal convergence when max(forces) and rms(forces) fall below the chosen threshold, divided by this factor. Convergence of max(step) and rms(step) is ignored.

• check_eigval_structure (bool, default: False) -- Check the eigenvalues of the modes we maximize along. Convergence requires them to be negative. Useful if TS searches are started from geometries close to a minimum.

• restart_info (default: None) -- Restart information. Undocumented.

• check_coord_diffs (bool, default: True) -- Whether coordinates of chain-of-sates images are checked for being too similar.

• coord_diff_thresh (float, default: 0.01) -- Unitless threshold for similary checking of COS image coordinates. The first image is assigned 0, the last image is assigned to 1.

• fragments (Optional[Tuple], default: None) -- Tuple of lists containing atom indices, defining two fragments.

• monitor_frag_dists (int, default: 0) -- Monitor fragment distances for N cycles. The optimization is terminated when the interfragment distances falls below the initial value after N cycles.

• out_dir (str, default: '.') -- String poiting to a directory where optimization progress is dumped.

• h5_fn (str, default: 'optimization.h5') -- Basename of the HDF5 file used for dumping.

• h5_group_name (str, default: 'opt') -- Groupname used for dumping of this optimization.

check_convergence(step=None, multiple=1.0, overachieve_factor=None)

Check if the current convergence of the optimization is equal to or below the required thresholds, or a multiple thereof. The latter may be used in initiating the climbing image.

dump_restart_info()

final_summary()

fit_rigid(*, vectors=None, vector_lists=None, hessian=None)

get_path_for_fn(fn, with_prefix=True)

get_restart_info()

log(message, level=50)

make_conv_dict(key, rms_force=None, rms_force_only=False, max_force_only=False, force_only=False)

abstract optimize()

postprocess_opt()

prepare_opt()

print_opt_progress(conv_info)

procrustes()

Wrapper for procrustes that passes additional arguments along.

report_conv_thresholds()

run()

scale_by_max_step(steps)

set_restart_info(restart_info)

write_cycle_to_file()

write_image_trjs()

write_results()

write_to_out_dir(out_fn, content, mode='w')

pysisyphus.optimizers.Optimizer.get_data_model(geometry, is_cos, max_cycles)

class pysisyphus.optimizers.HessianOptimizer.HessianOptimizer(geometry, trust_radius=0.5, trust_update=True, trust_min=0.1, trust_max=1, max_energy_incr=None, hessian_update='bfgs', hessian_init='fischer', hessian_recalc=None, hessian_recalc_adapt=None, hessian_xtb=False, hessian_recalc_reset=False, small_eigval_thresh=1e-08, line_search=False, alpha0=1.0, max_micro_cycles=25, rfo_overlaps=False, **kwargs)

Bases: Optimizer

__init__(geometry, trust_radius=0.5, trust_update=True, trust_min=0.1, trust_max=1, max_energy_incr=None, hessian_update='bfgs', hessian_init='fischer', hessian_recalc=None, hessian_recalc_adapt=None, hessian_xtb=False, hessian_recalc_reset=False, small_eigval_thresh=1e-08, line_search=False, alpha0=1.0, max_micro_cycles=25, rfo_overlaps=False, **kwargs)

Baseclass for optimizers utilizing Hessian information.

Parameters:

• geometry (Geometry) -- Geometry to be optimized.

• trust_radius (float, default: 0.5) -- Initial trust radius in whatever unit the optimization is carried out.

• trust_update (bool, default: True) -- Whether to update the trust radius throughout the optimization.

• trust_min (float, default: 0.1) -- Minimum trust radius.

• trust_max (float, default: 1) -- Maximum trust radius.

• max_energy_incr (Optional[float], default: None) -- Maximum allowed energy increased after a faulty step. Optimization is aborted when the threshold is exceeded.

• hessian_update (Literal['none', None, False, 'bfgs', 'damped_bfgs', 'flowchart', 'bofill', 'ts_bfgs', 'ts_bfgs_org', 'ts_bfgs_rev'], default: 'bfgs') -- Type of Hessian update. Defaults to BFGS for minimizations and Bofill for saddle point searches.

• hessian_init (Literal['calc', 'unit', 'fischer', 'lindh', 'simple', 'swart', 'xtb', 'xtb1', 'xtbff'], default: 'fischer') -- Type of initial model Hessian.

• hessian_recalc (Optional[int], default: None) -- Recalculate exact Hessian every n-th cycle instead of updating it.

• hessian_recalc_adapt (Optional[float], default: None) -- Use a more flexible scheme to determine Hessian recalculation. Undocumented.

• hessian_xtb (bool, default: False) -- Recalculate the Hessian at the GFN2-XTB level of theory.

• hessian_recalc_reset (bool, default: False) -- Whether to skip Hessian recalculation after reset. Undocumented.

• small_eigval_thresh (float, default: 1e-08) -- Threshold for small eigenvalues. Eigenvectors belonging to eigenvalues below this threshold are discardewd.

• line_search (bool, default: False) -- Whether to carry out a line search. Not implemented by a subclassing optimizers.

• alpha0 (float, default: 1.0) -- Initial alpha for restricted-step (RS) procedure.

• max_micro_cycles (int, default: 25) -- Maximum number of RS iterations.

• rfo_overlaps (bool, default: False) -- Enable mode-following in RS procedure.

• **kwargs -- Keyword arguments passed to the Optimizer baseclass.

filter_small_eigvals(eigvals, eigvecs, mask=False)

get_alpha_step(cur_alpha, rfo_eigval, step_norm, eigvals, gradient)

get_augmented_hessian(eigvals, gradient, alpha=1.0)

static get_newton_step(eigvals, eigvecs, gradient)

get_newton_step_on_trust(eigvals, eigvecs, gradient, transform=True)

Step on trust-radius.

See Nocedal 4.3 Iterative solutions of the subproblem

get_rs_step(eigvals, eigvecs, gradient, name='RS')

static get_shifted_step_trans(eigvals, gradient_trans, shift)

get_step_func(eigvals, gradient, grad_rms_thresh=0.01)

housekeeping()

Calculate gradient and energy. Update trust radius and hessian if needed. Return energy, gradient and hessian for the current cycle.

log_negative_eigenvalues(eigvals, pre_str='')

prepare_opt(hessian_init=None)

property prev_eigvec_max

property prev_eigvec_min

static quadratic_model(gradient, hessian, step)

reset()

rfo_dict = {'max': (-1, 'max'), 'min': (0, 'min')}

static rfo_model(gradient, hessian, step)

save_hessian()

set_new_trust_radius(coeff, last_step_norm)

solve_rfo(rfo_mat, kind='min', prev_eigvec=None)

update_hessian()

update_trust_radius()

pysisyphus.optimizers.HessianOptimizer.dummy_hessian_update(H, dx, dg)

class pysisyphus.optimizers.RFOptimizer.RFOptimizer(geometry, line_search=True, gediis=False, gdiis=True, gdiis_thresh=0.0025, gediis_thresh=0.01, gdiis_test_direction=True, max_micro_cycles=25, adapt_step_func=False, **kwargs)

Bases: HessianOptimizer

__init__(geometry, line_search=True, gediis=False, gdiis=True, gdiis_thresh=0.0025, gediis_thresh=0.01, gdiis_test_direction=True, max_micro_cycles=25, adapt_step_func=False, **kwargs)

Rational function Optimizer.

Parameters:

• geometry (Geometry) -- Geometry to be optimized.

• line_search (bool, default: True) -- Whether to carry out implicit line searches.

• gediis (bool, default: False) -- Whether to enable GEDIIS.

• gdiis (bool, default: True) -- Whether to enable GDIIS.

• gdiis_thresh (float, default: 0.0025) -- Threshold for rms(forces) to enable GDIIS.

• gediis_thresh (float, default: 0.01) -- Threshold for rms(step) to enable GEDIIS.

• gdiis_test_direction (bool, default: True) -- Whether to the overlap of the RFO step and the GDIIS step.

• max_micro_cycles (int, default: 25) -- Number of restricted-step microcycles. Disabled by default.

• adapt_step_func (bool, default: False) -- Whether to switch between shifted Newton and RFO-steps.

• **kwargs -- Keyword arguments passed to the Optimizer/HessianOptimizer baseclass.

optimize()

postprocess_opt()