Saddle-points#

Directory: libfly/saddle

The saddle folder contains a collection of utilities for saddle-point finding. Everything is contained within the namespace fly::saddle.

A (strictly first-order) saddle-point is a stationary point where the Hessian matrix has exactly one negative eigenvalue.

Master saddle-point finder#

File: libfly/saddle/find.hpp

Class to coordinate a group of saddle-point searches.

Saddle point (SP) searching is an involved process below is an outline of the stages:

For each central atom i in the input the following steps are carried out.
Saddle point searches are carried out (in batches) until a stopping criterion (one of: reconstruction failure, consecutive failed search limit or hitting a total search limit) is reached. Each SP search requires the following series of steps.
A random perturbation of atomic positions and the dimer axis (centred on i) is made.
The dimer method is used to move the perturbed state to a SP of the potential. It uses a cache of previously discovered saddle points (SPs)to avoid previous SPs by cancelling the search if the angle between the current search and the historical SP falls below theta. The angle theta is a decaying function of the number of SP searches.
It is ensured that the maximally displaced atom is i and the minimally displaced atom is frozen (if there are translational degrees of freedom).
The LBFGS minimiser is used to find the adjacent minima, due to freezing the atom in the previous step we can guarantee no translation in the min->sp->min pathway.
It is ensured than the pathway starts at the initial basin and that all the stationary points are distinct points in space.
A local mechanism is constructed from the displacements of the atoms around the central atom.
It is ensured that the local mechanism can be reconstructed in the orientation and location it was discovered. If not it is marked as a poisoned mechanism.
It is ensured that the mechanism has not been encountered before.
The set of mechanisms related by the symmetries of the local environment are constructed and it is ensured that when these mechanisms are reconstructed the energy barriers are accurate.

class Master#

Coordinate saddle-point searches on a (compute) node.

The Master coordinates saddle-point finding for a set of openMP threads, typically one Master should be instantiated per compute node. Each Master stores all the threads reusable objects: minimiser, prng, etc.

Public Types

using SoA = system::SoA<Position const&, Frozen const&, TypeID const&>#: Convenience alias.

Public Functions

Master(Options const &opt, system::Box const &box, system::TypeMap<Mass> map, potential::Generic const &pot, minimise::LBFGS const &min, saddle::Dimer const &dimer)#

Construct a new Master Finder object to manage opt.num_threads openMP threads.

This will create opt.num_threads copies of each of the input parameters, one for each thread.

Parameters:

opt – The configuration options.
box – Description of the simulation space.
map – The mass map for mass-weighting the Hessians.
pot – The potential to use for minimisations and SP searches.
min – The minimiser to relax either side of a SP.
dimer – Used to find SPs.

auto find_mechs(std::vector<LocalisedGeo> const &geos, SoA in, std::optional<Hint> const &hint = {}) -> std::vector<Found>#

Find all the mechanisms centred on the unknown geometries.

This will recursively schedule SP searches on the slave threads.

Parameters:

geos – A list of localised geometries encoding the atoms to centre the SP searches on and their local environments.
in – Description of system to search in.
hint – A hint to aid in

inline auto get_options() const noexcept -> Options const&#: Get the options object.

Public Static Functions

static auto package(std::vector<int> const &ix, env::Catalogue const &cat, int num_threads = 1) -> std::vector<LocalisedGeo>#

Transform a list of atom indices, ix, into a list of LocalisedGeo for use in find_mechs().

Parameters:

ix – List of indexes of environments which the mechanisms must be centred on.
cat – Catalogue in ready state.
num_threads – Number of openMP threads to use for this operation.

class Found#

A collection of mechanisms centred on a central atom.

Public Functions

inline std::vector<env::Mechanism> &mechs()#

Fetch the mechanisms on this->centre.

Throws if SPS failed.

inline explicit operator bool() const noexcept#

Truthy if search was successful.

If during the SPS, a reconstruction was attempted onto a supposedly symmetric state and the dimer-method failed to converge that reconstruction to a SP then the searches in that environment will be cancelled and this will return false.

This probably occurred because the symmetry tolerance was not tight enough e.g. this geometry was associated to the wrong reference LE. The best course of action is to tighten the corresponding reference LE’s delta_max, reclassify system and re-run the SP searches.

struct Hint#

A hint for ensuring ergodic pathways.

Public Members

Index::scalar_t centre#: Central atom of mech that brought us to current state.

system::VoS<Delta> delta_sp#: To SP of mechanism that brought us to current state.

env::Geometry<Index> geo#: Geometry around central atom in previous state.

system::SoA<Position> prev_state#: A reference to the previous state of the system.

struct LocalisedGeo#

A structure which packages the data require to perform saddle point searches.

Public Members

Index::scalar_t centre#: Central atom.

env::Geometry<Index> geo#: Geometry around central atom.

std::vector<env::Catalogue::SelfSymetry> syms#: Symmetries of geo onto itself.

struct Options#

Configure the saddle-point finding algorithm.

Public Members

double basin_tol = 0.05#: Tolerance for basins to be considered distinct.

int batch_size = std::min(max_failed_searches, num_threads)#: Number of simultaneous SP searches per new environment.

double capture_E_tol = 0.01#

Maximum energy error between reconstructed+relaxed minima/saddle-points and their original discovery.

After a min->sp->min pathway is found it is: converted into a mechanism; reconstructed onto the initial minima and then relaxed. These relaxed SP/minima must have energies within capture_E_tol of the pathway’s sp and final min.

double capture_r_tol = 0.05#

Maximum distance error between reconstructed+relaxed minima/saddle-points and their original discovery.

After a min->sp->min pathway is found it is: converted into a mechanism; reconstructed onto the initial minima and then relaxed. These relaxed SP/minima must be within capture_r_tol of the pathway’s sp and final min.

bool debug = false#: Print out debugging info.

io::BinaryFile *fout = nullptr#

If provided write debugging data here.

It is the users responsibility to ensure the lifetime of fout is at least as long as the lifetime of the this object and ensure only a single thread writes to fout at any one time. This really only exist for debugging purposes…

double freeze_tol = 0.01#

Tolerance for minimally displaced atom to be frozen (to remove translational error accumulation).

If the minimally displaced atom displaces further than this then an error will be raised.

double hessian_eigen_zero_tol = 1e-6#: Absolute eigen values of the mass-weighted hessian matrix, smaller than this, are considered zero.

int max_failed_searches = 100#: Maximum number of consecutive failed searches per environment.

int max_searches = 250#: Maximum number of searches per environment.

double mech_tol = 0.1#: Tolerance for mechanisms to be considered distinct.

double nudge_frac = 0.01#: Fraction of distance between initial position and SP that dimer is displaced along its axis before minimisation.

int num_threads = omp_get_max_threads()#: Number of openMP threads to dispatch saddle-point finding to.

double r_pert = 4#

The cut-off radius of the random perturbation (centred on an atom).

The initial system randomly perturbed to produce a starting-point for saddle-point finding. The perturbation is Gaussian in each coordinate axis with standard deviation stddev. An envelope function will linearly decrease the size of each atoms perturbation based off its distance to a central atom. Only atoms closer than r_pert will be perturbed.

double recon_e_tol_abs = 0.01#: The absolute energy tolerance for a mechanisms energy barrier & delta when reconstructed onto symmetries.

double recon_e_tol_frac = 0.01#: The fractional energy tolerance for a mechanisms energy barrier & delta when reconstructed onto symmetries.

double recon_norm_abs_tol = 0.25#: The absolute difference between the expected relaxation and the measured relaxation of a reconstructed mechanism onto its symmetries.

double recon_norm_frac_tol = 0.5#: The fractional difference between the expected relaxation and the measured relaxation of a reconstructed mechanism onto its symmetries.

double stationary_tol = basin_tol / 2#: Tolerance for stationary points to be considered distinct.

double stddev = 0.6#: The standard deviation of the random perturbations, see r_pert.

Dimer#

File: libfly/saddle/dimer.hpp

Saddle-point finding with the dimer method.

class Dimer#

Saddle point locator.

A dimer is two images of a system. These images have a centre point and a unit axis. The images are (conceptually) located delta_r in the plus/minus directions along the axis.

This implementation uses the superlinear dimer method to find saddle-points. This alternates optimizing the orientation of the dimer along the minimum eigen-mode then translating the dimer along an effective potential to find the SP.

Public Types

enum Exit#

Return codes for step()`.

Values:

enumerator success#: Found saddle-point.

enumerator convex#: Failed: stuck in +curvature region.

enumerator iter_max#: Exceeded the maximum number of iterations but in -curvature region.

enumerator collision#: Approached a known SP during search.

enumerator uninit#: Stand in for uninitialised exit code.

Public Functions

Dimer(Dimer&&) = default#: Move construct a new Dimer object.

Dimer(Dimer const&) = default#: Copy construct a new Dimer object.

inline Dimer(Options const &opt, Rotor::Options const &r_opt, system::Box const &box)#

Construct a new Dimer object.

Parameters:

opt – Options for translations pass.
r_opt – Options for rotation pass.
box – Description of simulation space.

auto find_sp(system::SoA<Position&, Axis&> out, system::SoA<Position const&, Axis const&, TypeID const&, Frozen const&> in, system::SoA<Position const&> in_min, potential::Generic &pot, std::vector<system::SoA<Position>> const &hist_sp, double theta_tol, int count_frozen, int num_threads = 1) -> Exit #

Advance the dimer towards a saddle-point.

Performs at most max_find_sp translation steps. The translation uses the LBFGS algorithm with a trust-radius limited step size and an early-exit condition if the curvature is positive for too long.

Parameters:

out – Final position and axis orientation written here.
in – Initial position, axis and per-particle data forwarded to potential.
pot – Potential energy function.
in_min – initial minima that dimer is climbing FROM.
hist_sp – Previously discovered saddle points.
theta_tol – If the angle (Radians) between this SPS and a previous SP is less than theta_tol then abort.
count_frozen – The number of frozen atoms in in.
num_threads – Number of openMP threads to use.

Returns:

Exit-code.

struct Options#

Configure the dimers internal translation minimisation pass.

Public Members

double convex_max = 3#: If convex_max steps with +Ve curvature then exit early.

bool debug = false#: Print out debug info.

double f2norm = 1e-5#: Force convergence criterion (eV/Angstroms).

io::BinaryFile *fout = nullptr#

If provided at each frame of the translation out will be written to this file.

It is the users responsibility to ensure the lifetime of fout is at least as long as the lifetime of the LBFGS object and ensure only a single LBFGS object writes to fout at any one time. This really only exist for debugging purposes…

double grow_trust = 1.5#: Trust radius expansion rate.

int hist_check_freq = 5#: Every hist_check_freq number of steps will check not converging to known SP.

int max_steps = 1000#: Maximum number of steps before failing.

double max_trust = 0.5#: Maximum trust radius e.g max steps size (Angstroms).

double min_trust = 0.05#: Minimum trust radius e.g initial step size (Angstroms).

int n = 10#: L-BFGS core translation-history size.

double proj_tol = 0#: Trust tolerance, set larger to reduce trust radius change.

double shrink_trust = 0.5#: Trust radius contraction rate.

double skin_frac = 1.1#

Used to determine the skin size.

The skin thickness is determined such that the expected number of neighbours is skin_frac * (num_neighbours if no skin used). Hence skin_frac must be larger than one. Neighbour lists are built less often when skin_frac is higher but there will be more more non-neighbours in neighbour lists.

bool use_history = true#: If true the search will abort if it approaches too close to to a previous SP.

Rotor#

File: libfly/saddle/rotor.hpp

Class to orient a dimer, exposes part of the Dimer classes implementation for re-use.

class Rotor#

Orients a dimer to align it with the minimum eigen-mode.

Public Functions

inline explicit Rotor(Options const &opt)#

Construct a new Rotor object.

Parameters:: opt – Options for rotor minimisation pass.

Rotor(Rotor&&) = default#: Move construct a new Rotor object.

Rotor(Rotor const&) = default#: Copy construct a new Rotor object.

auto eff_gradient(system::SoA<PotentialGradient&> out, system::SoA<Axis&> inout, system::SoA<TypeID const&, Frozen const&> in, potential::Generic &pot, neigh::List &nl, int count_frozen, int threads = 1) -> double#

Compute the effective potential’s gradient.

This is achieved by inverting the component of the gradient parallel to the minimum eigen-mode of the potential. The dimer is rotated using the LBFGS algorithm to make it align with the minimum eigen-mode.

Assumes the neighbour list is ready, force on frozen atoms will be zero.

Note

This function does actually modify neighbour but returns it to an identical state after it is done. It guarantees not to call neigh::List::rebuild().

Parameters:

in – Per-atom input properties
out – Write effective gradient here.
inout – The axis input and output.
nl – Neighbour list (in ready state i.e. neigh::List::update() or neigh::List::rebuild() called) configured with a cut-off at least r_cut().
pot – The potential energy function.
count_frozen – The number of frozen atoms in in.
threads – Number of openMP threads to use.

Returns:

The approximate curvature of the wrapped potential along the output Axis.

inline auto r_cut(potential::Generic const &pot) const noexcept#

Compute the cut-off radius.

This is the cut-off that nl in call to eff_grad() with potential pot must be configured with.

Parameters:: pot – Potential that eff_grad() will be called with.
Returns:: auto A slightly larger cut-off that enables eff_grad to displace atoms a little without rebuilding the nl.

struct Options#

Configure the dimers rotation minimisation pass.

Public Members

bool debug = false#: Print out debug info.

double delta_r = 0.001#: Half dimer length.

int iter_max_rot = 20#: Maximum number of iterations during rotation minimization.

int n = 6#: L-BFGS core rotation-history size.

bool relax_in_convex = true#: If false when in convex region we return only the component parallel to the min mode.

double theta_tol = 1 * M_PI / 360.#: (Rad) rotation convergence criterion.