Environments

Directory: libfly/env

The env folder contains routines working with local environments. Everything is contained within the namespace fly::env. The local environment (LE) of an atom is the set of atoms within some neighbourhood. It is assumed in OLKMC that the mechanisms accessible to an atom are completely contained-within and solely a-function-of its LE.

Geometry

File: libfly/env/geometry.hpp

Utilities for manipulating local environments of atoms.

Functions

template<typename ...M>
auto fly::env::centroid(system::VoS<M...> const &x) noexcept -> Vec

Compute the centroid of a vector of atoms, x.

For a set of atoms with positions \(x_i\) this is defined as:

\[\frac{1}{N} \sum_{i=1}^{N}{x_i}\]

Parameters:

x – Set of atoms with the Position property.

template<typename Pr = Position, typename ...M1, typename ...M2>
auto fly::env::rmsd(system::VoS<M1...> const &x, system::VoS<M2...> const &y) -> double

Compute the RMSD between the positions of the atoms in x and y.

The RMSD (root mean squared distance) is the sqrt of the sum of the squared distances between atoms, for two sets of atoms with positions \(x_i\) and \(y_i\) this is:

\[\Delta_\text{RMSD} = \sqrt{\sum_{i=1}^N{\left\lVert x_i - y_i \right\rVert^2}}\]

Template Parameters:

Pr – The distance-like vector property to use as atomic position.

Parameters:

• x – Set of atoms with the Pr property.

• y – Set of atoms with the Pr property.

template<typename Pr = Position, typename E, typename ...M1, typename ...M2>
auto fly::env::grmsd(E const &M, system::VoS<M1...> const &x, system::VoS<M2...> const &y) -> double

Compute the generalised RMSD between two sets of atoms.

The generalised RMSD (root mean squared distance) between two sets of atoms, x and y, is the RMSD between them after applying the matrix transformation M to each atom in x. For two sets of atoms with positions \(x_i\) and \(y_i\) this is:

\[\Delta_\text{GRMSD} = \sqrt{\sum_{i=1}^N{\left\lVert Mx_i - y_i \right\rVert^2}}\]

Template Parameters:

Pr – The distance-like vector property to use as atomic position.

Parameters:

• x – Set of atoms with the Pr property.

• y – Set of atoms with the Pr property.

• M – Transformation matrix.

template<typename ...M1, typename ...M2>
auto fly::env::ortho_onto(system::VoS<M1...> const &x, system::VoS<M2...> const &y) -> Mat

Compute the orthogonal transformation matrix that best transforms x onto y.

This does not permute any of the atoms it just minimizes the RMSD between the two sets.

Effectively solves this optimisation problem:

\[\min_{O}{\sum_{i=1}^N{\left\lVert O x_i - y_i \right\rVert^2}} \quad \text{with} \quad OO^\intercal = I\]

With \(x_i\) and \(y_i\) the position of the i th atom in x onto y respectively.

Uses the “Kabsch algorithm” see: https://en.wikipedia.org/wiki/Kabsch_algorithm

Parameters:

• x – Set of atoms with the Position property.

• y – Set of atoms with the Position property.

template<typename ...M1, typename ...M2, typename F>
auto fly::env::for_equiv_perms(system::VoS<M1...> &mut, system::VoS<M2...> const &ref, double delta, int n, F const &f) -> void

Explore equivalent permutations of atoms in two sets.

This function is for greedily exploring the permutations of atoms in mut such that, after an orthogonal transformation generated by ortho_onto(), the grmsd() between the two sets of atoms is less than delta.

This effectively finds a matrix \(O\) and permutation \(\pi\) that satisfies:

\[\sum_{i=1}^N{\left\lVert O x_{\pi(i)} - y_i \right\rVert^2} < \delta^2\]

Subject to the constraints:

\[OO^\intercal = I \quad \text{and} \quad \pi(i) = i \,\,\,\, \forall \,\,\,\, i < n\]

With \(x_i\) and \(y_i\) the position of the i th atom in x onto y respectively.

Example:

#include "libfly/env/geometry.hpp"

#include "libfly/system/VoS.hpp"
#include "libfly/system/property.hpp"

using Set = fly::system::VoS<fly::Position, fly::Colour>;

template <typename... Args>
auto work(Args const&...) -> void {
  // Artificial work.
}

auto count_sym(Set const& ref, double delta) -> int {
  //
  Set mut = ref;  // Make a copy of the reference set

  int count = 0;

  // Explore permutations fixing the first atom.
  fly::env::for_equiv_perms(mut, ref, delta, 1, [&](fly::Mat const& O, double rmsd) {
    // If this lambda is called then "mut" has been permuted into an equivalent
    // permutation. "O" is the matrix that maps "mut" onto "ref" such that that:
    //
    //     "rmsd" == grmsd(O, mut, ref) and "rmsd" < delta.

    // We could now do something with "rmsd" and "O".

    work(O, rmsd);  // Suppresses warnings

    // We MUST NOT mutate "ref" or "mut".

    // If we "return true" then exploration of permutations terminates.
    // If we "return false" then exploration of permutations continues.

    // Here we are counting all the permutations so we will return false.
    count++;

    return false;
  });

  return count;  // The number of approximate symmetries that ref has.
}

for_equiv_perms() will not permute the first n (== 1 in the above example) atoms in mut. This is useful if there exist a number of atoms who’s permutation is known.

Parameters:

• mut – Set of atoms with the Position & Colour properties that will be permuted to match ref.

• ref – The reference set of atoms with the Position & Colour properties.

• delta – The maximum grmsd() between the permuted mut and ref

• n – The first n atoms in mut will not be permuted.

• f – An invokable invoked with signature (fly::Mat const& O, double rmsd) -> bool called for each permutation of mut that matches ref, see the example above for details.

Geometry class

template<typename ...Pr>
class Geometry : public fly::system::VoS<Position, Colour, Pr...>

A Geometry models a local distribution of periodically resolved Atoms.

The distribution of atoms is centred on the first atom, i.e. the first atom cannot move and is always fixed during permutation operations and the centroid of the geometry is the origin.

A Geometry is a fly::system::VoS with the default properties Position and Colour.

Public Functions

template<typename ...Ts>
inline std::optional<GeoInfo> best_perm_onto(system::VoS<Ts...> const &other, double delta, Geometry *scratch_best = nullptr)

Find the best permutation of the atoms in this Geometry onto other.

Parameters:

• other – The reference set of atoms with the Position & Colour properties.

• delta – The maximum grmsd() between the permuted this and ref for this operation to succeed.

• scratch_best – An optional pointer to scratch space that the algorithm can use (otherwise this function will allocate).

Returns:

std::optional<GeoInfo> An engaged GeoInfo if a permutation was found.

inline void centre() noexcept

Shift the origin to the centroid of the atoms.

template<typename ...Ts>
inline std::optional<GeoInfo> permute_onto(system::VoS<Ts...> const &other, double delta)

Find the first permutation of the atoms in this Geometry onto other within tolerance delta.

Parameters:

• other – The reference set of atoms with the Position & Colour properties.

• delta – The maximum grmsd() between the permuted this and ref if this operation succeeds.

Returns:

std::optional<GeoInfo> An engaged GeoInfo if a permutation was found.

template<class Archive>
inline void serialize(Archive &archive)

Lib cereal serialization support.

struct GeoInfo

Returned by permutation methods, contains extra info about the permutation.

Public Members

Mat O

Orthogonal transformation required to map this onto x.

double rmsd

Equal to grmsd(O, *this, x).

Heuristics

File: libfly/env/heuristics.hpp

Utilities for quickly comparing local environments without resorting to a full Geometry::permute_onto().

Fingerprints

class Fingerprint

An ordered representation of the intra-atomic distances in a Geometry.

The fingerprint is a continuous representation that is cheap(ish)ly comparable.

Public Functions

auto equiv(Fingerprint const &other, double delta) const -> bool

A fast test to see if two local environments may be equivalent.

Explicitly, for the LE this fingerprint represents, this function must return true if the LE is able to be permute_onto() other with the same delta.

Parameters:

• other – The other fingerprint.

• delta – The tolerance for the equivalence.

Returns:

true If two local environments may be equivalent.

Returns:

false If two local environments are not equivalent.

inline auto r_min() const -> double

Get the smallest intra-atomic separation in this environment.

template<typename ...T>
inline auto rebuild(Geometry<T...> const &geo) -> void

Rebuild the fingerprint from a Geometry.

template<class Archive>
inline void serialize(Archive &archive)

Lib cereal serialization support.

Graph hash

std::size_t fly::env::canon_hash(Geometry<Index> &geo, double r_edge, std::size_t c_max, Geometry<Index> *scratch = nullptr)

Reorder a geometry into a canonical order and produce a hash encoding the colours and topology of the geometry.

This functions builds a graph based representation of the geometry, encoding atoms as coloured nodes joined with edges if the atoms are closer than r_edge. This graph is then canonised using the nauty library (see: https://pallini.di.uniroma1.it/) and a hash is derived from the canonical adjacency matrix.

Parameters:

• geo – The geometry to be reordered.

• r_edge – The edge length used to determine if atoms in the geometry are bonded in the graph representation.

• c_max – The number of colours in the simulation, equal to 2 * TypeMap::num_types().

• scratch – An optional parameter, if provided and scratch->size() >= geo.size() then no allocation will occur.

Mechanisms

File: libfly/env/mechanisms.hpp

Representations of mechanisms.

class Mechanism

A representation of a local mechanism containing the the displacements of the atoms in a LE.

A mechanism describes the displacement vectors from a LE to an adjacent minima via a SP.

Frozen atoms must have a displacement equal to Vec::Zero().

Public Functions

template<class Archive>
inline void serialize(Archive &archive)

Lib cereal serialization support.

Public Members

system::VoS<Axis> axis = {}

The axis of the dimer at the saddle-point, this is normalised.

double barrier = -1

Energy barrier (eSp - e0).

double delta = -1

Energy change (eF - e0).

system::VoS<Delta> delta_fwd = {}

Displacement vectors from initial to final.

system::VoS<Delta> delta_sp = {}

Displacement vectors from initial to saddle-point.

double err_fwd = std::numeric_limits<double>::max()

Forward reconstruction error.

double err_sp = std::numeric_limits<double>::max()

Saddle point reconstruction error.

double kinetic_pre = -1

Arrhenius pre factor.

bool poison_fwd = false

If true this mechanism’s final state cannot be reconstructed.

bool poison_sp = false

If true this mechanism’s SP cannot be reconstructed.

The catalogue

File: libfly/env/catalogue.hpp

Classes responsible for mapping local environments to mechanisms.

class Catalogue

The Catalogue is a mapping from Geometry to local environments and mechanisms.

The Catalogue is responsible for building geometries from a Supercell and associating them with equivalent, historically-seen geometries via a three step matching process detailed in Catalogue::rebuild().

Public Functions

inline explicit Catalogue(Options const &opt)

Construct a new Catalogue object.

inline Catalogue(Options const &opt, std::ifstream &fin)

Load a catalogue from a binary archive.

auto calc_self_syms(int i) const -> std::vector<SelfSymetry>

Compute the permutation and transformation of the ith geometry that maps it onto itself.

This is according to the delta_max() of the environment it is currently equivalent to.

inline auto dump(std::ofstream &fout) const -> void

Write this catalogue to a binary archive.

inline auto get_geo(int i) const -> Geometry<Index> const&

Get a view of the (indexed) geometry, in canonical order, around atom i.

inline auto get_ref(int i) const -> Env const&

Get the reference geometry stored in the catalogue that is equivalent to the geometry around atom i.

inline std::size_t num_keys() const noexcept

Get the number of discrete keys in the catalogue.

void optimize()

Re-order the catalogue into frequency order.

template<typename Map, typename ...T>
inline auto rebuild(system::Supercell<Map, T...> const &cell, int num_threads = 1) -> std::vector<int>

Build the geometries and associate them with equivalent, historically-seen geometries.

This occurs via a three step matching process:

1. The new geometries are canonised (based on their graph representation) and a discrete key generated.

2. Geometries with equal discrete keys are compared for equivalence of their fingerprints.

3. If they are equivalent under (2) they are compared for equivalence via geometry::permute_onto().

Any new geometries that have not been encountered before are inserted into the catalogue and the index of the atom atom which the unknown geometry is centred on is returned.

This process is deterministic, providing none of the reference geometries are modified (allowing new ones to be inserted) the same match will always be returned.

Parameters:

• cell – The supercell, must have the Position, TypeID and Frozen properties.

• num_threads – The number of openMP threads to use.

Returns:

A std::vector of indices corresponding to the new geometries.

template<typename Map, typename ...T>
inline auto reconstruct(system::SoA<Position&> out, Mechanism const &mech, int i, system::Supercell<Map, T...> const &in, bool in_ready_state, int num_threads) -> Mat

Reconstruct a mechanisms mech onto the i atom of in.

Reconstruct a mechanisms.

Parameters:

• out – Write the reconstructed position here.

• mech – The mechanism to reconstruct.

• i – The index of the atom to reconstruct the mechanism onto.

• in – The initial state of the system before the reconstruction.

• in_ready_state – If true this function will assume the currently loaded geo/ref of the ith atom match the input in otherwise, the geometry will be rebuilt.

• num_threads – Number of openMP threads to use.

Returns:

Mat The matrix that transforms mech before reconstruction onto out.

auto refine_tol(int i, double min_delta = 0) -> double

Tighten the tolerance of the ith environment.

This is done such that the current geometry no longer matches the reference or delta_max = min_delta. Additionally, this function resets the frequency and false positive counters of the reference environment.

Returns:

The new delta_max of the environment reference environment that atom i was equivalent to.

template<class Archive>
inline void serialize(Archive &archive)

Lib cereal serialization support.

auto set_mechs(int i, std::vector<Mechanism> const &m) -> void

Set the vector of Mechanisms.

Pre-condition, this may only be called once per LE.

inline int size() const noexcept

Get the number of local environments in the catalogue.

struct Env : private fly::env::Geometry<>

The representation of a local environment in the catalogue.

Public Functions

Env() = default

For lib cereal.

inline auto cat_index() const noexcept -> int

Get an index (unique to this environment) in the catalogue.

inline auto delta_max() const noexcept -> double

Get the internal maximum norm between this environment and a geometry for them to be considered equivalent.

inline auto get_mechs() const -> std::vector<Mechanism> const&

Get a vector of the Mechanisms centred on this environment.

Pre-condition, set_mechs() must have been called on this Env.

inline auto ref_geo() const -> Geometry<> const&

Fetch the reference geometry this environment represents.

template<class Archive>
inline void serialize(Archive &archive)

Lib cereal serialization support.

struct Options

Used to configure the catalogue.

Public Functions

template<class Archive>
inline void serialize(Archive &archive)

Lib cereal serialization support.

Public Members

bool debug = false

If true prints debugging info.

double delta_max = std::numeric_limits<double>::max()

Initial maximum difference in L2 norm between LEs for them to be considered the same.

double min_delta_max = 1e-7

Minimum value of delta_max during refinement (smaller is considered an error).

double overfuzz = 0.5

Smaller to decrease false positives during fingerprint equivalence but, values < 1.0 introduce false negatives.

double r_edge = 3.0

Maximum distance for atoms to be considered connected in the canonisation neighbour graph.

double r_env = 5.2

Radius of a local environment.

struct SelfSymetry

Describes a transformation + permutation that maps a geometry onto itself.

Explicitly this a the transformation such that:

Public Members

Mat O

Transformation matrix.

std::vector<Index::scalar_t> perm

Permutation.