tools_geometry — building product end states¶

Geometry manipulation for building reaction end states.

A NEB needs a product as well as a reactant, and for anything bigger than a single bond rearrangement the product is the awkward one to draw by hand. The functions here build one instead: get_dimer_bonded_cluster_indices() works out which atoms form each half of a stacked dimer, and get_best_flip_and_face_bases() swaps those halves over to give the flipped structure, ready to pass to prepare_neb().

swap_bonding_configuration() does the same job for proton transfers, moving one or more hydrogens across their hydrogen bonds. Where those two build a product out of what the reaction is known to do, seed_product_from_ts() builds one out of a transition state instead, stepping past the saddle along the geodesic that reaches it – no calculator, and no need to have guessed the mechanism first. seed_minima_from_ts() needs no end state at all: it rattles the saddle at random and relaxes what comes out, from enough directions to fall into both of the basins it connects, and reports them as a SeedSummary. For structures that already describe the same atoms, align_atom_sets() superposes one on the other with the optimal rigid Kabsch transform, and atom_set_rmsd() measures what remains.

kabsch_transform(mobile_positions: Sequence[Sequence[float]] | ndarray, reference_positions: Sequence[Sequence[float]] | ndarray, weights: Sequence[float] | ndarray | None = None) → tuple[ndarray, ndarray][source]¶

Find the best proper rigid transform from one point set to another.

The points correspond by row. The returned transform uses NumPy’s row vector convention, so it is applied as mobile_positions @ rotation + translation. Reflections are excluded: the rotation is always right-handed, even when a reflected fit would have a lower residual.

Parameters:

mobile_positions (array_like) – Corresponding Cartesian coordinates, each with shape (n, 3).
reference_positions (array_like) – Corresponding Cartesian coordinates, each with shape (n, 3).
weights (array_like, optional) – One non-negative weight per correspondence. At least one must be positive. By default every point has equal weight.

Returns:

tuple of numpy.ndarray – (rotation, translation) with shapes (3, 3) and (3,).

Raises:

TypeError – If positions or weights cannot be converted to numeric arrays.
ValueError – If the coordinate shapes differ, no positions are supplied, a value is non-finite, or the weights are invalid.

Rigidly superpose one atom set on a corresponding reference set.

The selected atoms determine the least-squares Kabsch fit, while the resulting rotation and translation are applied to every atom in mobile. This keeps the mobile structure internally rigid and makes it possible to fit on a stable substructure while carrying spectators or a flexible region along with it. Neither input is modified.

Correspondence is positional: the first mobile selection index is matched to the first reference selection index, and so on. Atom identities are not automatically matched or reordered.

Parameters:

mobile (ase.Atoms) – Structure to move and structure to fit it onto.
reference (ase.Atoms) – Structure to move and structure to fit it onto.
mobile_indices (int or iterable of int, optional) – Corresponding selections used for the fit. Each defaults to every atom in its structure and the selections must have equal length.
reference_indices (int or iterable of int, optional) – Corresponding selections used for the fit. Each defaults to every atom in its structure and the selections must have equal length.
weights ({None, "masses"} or array_like, optional) – Fit weights. None gives every correspondence equal influence; "masses" uses the selected mobile atoms’ masses; a numeric sequence supplies one non-negative weight per correspondence.

Returns:

ase.Atoms – A copy of mobile with every position transformed into the reference frame.

Raises:

TypeError – If either structure is not an ase.Atoms, or selections or weights have the wrong type.
ValueError – If the selections differ in length, are empty or repeated, or weights are invalid.
IndexError – If a selection contains an atom index outside its structure.

Calculate the RMSD between corresponding atoms, optionally after fitting.

Parameters are the same as align_atom_sets(). When align is True, the optimal rigid transform is applied to the selected mobile coordinates before the RMSD is measured; the input structures still are not modified.

Parameters:

mobile (ase.Atoms) – Structures whose corresponding positions are compared.
reference (ase.Atoms) – Structures whose corresponding positions are compared.
mobile_indices (int or iterable of int, optional) – Corresponding selections to compare. Each defaults to every atom.
reference_indices (int or iterable of int, optional) – Corresponding selections to compare. Each defaults to every atom.
weights ({None, "masses"} or array_like, optional) – RMSD and fit weights, interpreted as in align_atom_sets().
align (bool, optional) – Remove the best rigid rotation and translation before measuring. Defaults to False so the function can also report displacement in the current coordinate frame.

Returns:

float – Root-mean-square Cartesian distance in Å.

bonded_cluster_indices_no_anchor_hub(atoms: Atoms, anchor: int, mult: float = 1.0, multi_h: float = 1.3) → list[int][source]¶

Collect the atoms bonded to an anchor without routing through it.

A plain flood fill over the bonded graph would leak through the anchor into whatever else it touches, which for a stacked dimer means swallowing both halves at once. Here the anchor is marked visited before the walk starts, so the search reaches its immediate neighbours and then spreads outwards from them only — the anchor is a starting point, never a thoroughfare.

Parameters:

atoms (ase.Atoms) – Structure whose bonded graph is walked.
anchor (int) – Index of the anchor atom, included in the returned cluster.
mult (float, optional) – Scale applied to ASE’s natural covalent cutoffs. Raise it to make bonding more permissive.
multi_h (float, optional) – Separate cutoff scale for hydrogen, applied to its covalent radius. Hydrogen’s natural cutoff is tight enough to miss hydrogen bonds, which are exactly the contacts that matter across a base pair.

Returns:

list of int – Sorted indices of the anchor and everything bonded to it, in ascending order.

Raises:

IndexError – If anchor is out of range.

get_dimer_bonded_cluster_indices(atoms: Atoms, anchors: list[int], mults: Sequence[float] | None = None, multi_h: float = 1.3) → list[int][source]¶

Collect the atoms of both halves of a dimer, as one index list.

Runs bonded_cluster_indices_no_anchor_hub() from each anchor and merges the results, so an atom bridging the two halves is counted once.

Parameters:

atoms (ase.Atoms) – Structure whose bonded graph is walked.
anchors (list of int) – Exactly two anchor indices, one per half.
mults (sequence of float, optional) – Per-anchor cutoff scales. Defaults to [1.0, 1.0].
multi_h (float, optional) – Hydrogen cutoff scale, shared by both walks.

Returns:

list of int – Sorted union of the two clusters.

Raises:

ValueError – If anchors or mults does not hold exactly two values.

flip_and_face_bases(atoms: Atoms, baseA_idxs: list, baseB_idxs: list, anchors: list, rot_matrix: list | None = None) → Atoms[source]¶

Swap two fragments over, each landing on the other’s anchor.

Fits a local frame to each fragment, maps one frame onto the other through a reflection, and moves each fragment onto the opposite anchor. The reflection is what makes the two end up facing each other rather than merely translated, so the result is a plausible flipped end state rather than a structure with the fragments back to back.

Parameters:

atoms (ase.Atoms) – Structure to rearrange. Not modified; a copy is returned.
baseA_idxs (list of int) – Indices of the first fragment.
baseB_idxs (list of int) – Indices of the second fragment.
anchors (list of int) – Two anchor indices, one per fragment, in the same order as the fragments. Fragment A lands on anchor B and vice versa.
rot_matrix (list of float, optional) – Diagonal of the reflection applied between the two frames, as three signs. Defaults to [-1.0, 1.0, -1.0]. Which one is right depends on the geometry; get_best_flip_and_face_bases() searches them.

Returns:

ase.Atoms – Copy of atoms with the two fragments swapped. Atoms outside both fragments keep their positions.

optimize_with_fixed_anchors(atoms: ~ase.atoms.Atoms, baseA_idxs: list, baseB_idxs: list, anchor_indices: list, calc: ~typing.Any, fmax: float = 0.05, steps: int = 1000, raise_on_unconverged: bool = False, optimiser: ~collections.abc.Callable[[...], ~typing.Any] = <class 'ase.optimize.bfgs.BFGS'>, logfile: str | ~typing.TextIO | None = '-') → Atoms[source]¶

Relax the two fragments while holding their anchors still.

A structure straight out of flip_and_face_bases() is rigid-body exact but strained, because the fragments were moved as blocks. This relieves the strain without letting the fragments drift back, which is what pinning the anchors buys.

Only the fragment atoms are optimised; anything outside both fragments is left exactly where it was.

Parameters:

atoms (ase.Atoms) – Structure to relax. Not modified; a copy is returned.
baseA_idxs (list of int) – Indices of the first fragment.
baseB_idxs (list of int) – Indices of the second fragment.
anchor_indices (list of int) – Indices held fixed during the relaxation.
calc (ase.calculators.Calculator) – Calculator used for the relaxation.
fmax (float, optional) – Maximum force criterion in eV/Å.
steps (int, optional) – Maximum number of optimiser steps. A flipped structure needing more than the default has usually landed somewhere unphysical, which the ConvergenceWarning then says out loud instead of grinding on against ASE’s effectively unlimited default.
raise_on_unconverged (bool, optional) – Raise ConvergenceError instead of warning when the relaxation hits steps without reaching fmax.
optimiser (callable, optional) – ASE optimiser class to relax with, or anything callable as optimiser(atoms, logfile=...). Defaults to BFGS. See optimise_geom().
logfile (str, file object or None, optional) – Where the optimiser writes its per-step table. '-', the default, is stdout; a filename writes there instead; None silences it.

Returns:

ase.Atoms – Copy of atoms with the relaxed fragment positions written back, and info["converged"] recording whether the relaxation reached fmax.

Raises:

ConvergenceError – If the relaxation did not converge and raise_on_unconverged is True.

get_best_flip_and_face_bases(atoms: ~ase.atoms.Atoms, baseA_idxs: list, baseB_idxs: list, anchors: list, optimise_after: bool = True, calc: ~typing.Any = None, raise_on_unconverged: bool = False, optimiser: ~collections.abc.Callable[[...], ~typing.Any] = <class 'ase.optimize.bfgs.BFGS'>, logfile: str | ~typing.TextIO | None = '-') → Atoms[source]¶

Search the reflection signs for the tightest flipped structure.

flip_and_face_bases() takes a reflection whose correct signs depend on how the fragments happen to sit, and a wrong choice throws them apart. This tries every sign combination that is a genuine reflection and keeps whichever leaves the two fragment centres of mass closest together.

Parameters:

atoms (ase.Atoms) – Structure to rearrange. Not modified.
baseA_idxs (list of int) – Indices of the first fragment.
baseB_idxs (list of int) – Indices of the second fragment.
anchors (list of int) – Two anchor indices, one per fragment.
optimise_after (bool, optional) – Relax the result with optimize_with_fixed_anchors().
calc (ase.calculators.Calculator, optional) – Calculator for that relaxation. Required when optimise_after is True.
raise_on_unconverged (bool, optional) – Passed to optimize_with_fixed_anchors(); ignored when optimise_after is False, since nothing is then relaxed.
optimiser – Passed to optimize_with_fixed_anchors() for the single relaxation at the end; the reflection search itself moves the fragments rigidly and never touches the calculator. Ignored when optimise_after is False.
logfile – Passed to optimize_with_fixed_anchors() for the single relaxation at the end; the reflection search itself moves the fragments rigidly and never touches the calculator. Ignored when optimise_after is False.

Returns:

ase.Atoms – Flipped structure, relaxed when optimise_after is True, in which case info["converged"] records whether that relaxation reached fmax.

Raises:

ValueError – If optimise_after is True and no calculator is given.
ConvergenceError – If the relaxation did not converge and raise_on_unconverged is True.

swap_bonding_configuration(atoms: Atoms, donor_index: int | Iterable[int], hydrogen_index: int | Iterable[int], acceptor_index: int | Iterable[int]) → Atoms[source]¶

Swap one or more donor-H…acceptor bonds to donor…H-acceptor.

Builds the product end state of one or more proton transfers. Each hydrogen is moved to the acceptor side of its hydrogen bond, keeping the bond length it had to its donor, so the result is a sensible starting geometry for optimise_geom() and then a band.

A scalar donor or acceptor index is shared by all the hydrogens. Otherwise, donor and acceptor iterables must contain one index per hydrogen.

Parameters:

atoms (ase.Atoms) – Structure holding the hydrogen bond. Not modified; a copy is returned.
donor_index (int or iterable of int) – Index of a shared donor atom, or one donor index per hydrogen.
hydrogen_index (int or iterable of int) – Index or indices of the hydrogens being moved.
acceptor_index (int or iterable of int) – Index of a shared acceptor atom, or one acceptor index per hydrogen.

Returns:

ase.Atoms – Copy of atoms with each hydrogen on its acceptor side.

Raises:

TypeError – If an index argument is neither an integer nor an iterable of integers.
ValueError – If no hydrogen is supplied, donor or acceptor counts do not match the hydrogen count, a hydrogen is repeated, a transfer reuses an atom in more than one role, or a donor and acceptor occupy the same position.
IndexError – If an atom index is out of range.

exception SeedWarning[source]¶

Bases: UserWarning

A seeded end state did not come out looking like a new structure.

Warned rather than raised by seed_product_from_ts(), because a seed that stopped short of where it was aimed is still a starting point, and because whether it landed in the right basin is settled by relaxing it rather than by measuring it. Turn it into an error with warnings.simplefilter("error", SeedWarning) in a batch script, where a seed that is really a second copy of the reactant would otherwise be carried into everything downstream.

seed_minima_from_ts() warns it for the same reason: a rattle that found one side of a saddle rather than both has still found a minimum, and a search that has not answered the question it was run to ask should say so without discarding what it did find.

seed_product_from_ts(reactant: Atoms, ts: Atoms, n_images: int = 25, push: float = 1.0, n_steps: int = 10, tangent_images: int = 2, weights: str | Sequence[float] | ndarray | None = 'masses', clash_scale: float | None = 0.7, return_path: bool = False) → Atoms | tuple[Atoms, list[Atoms]][source]¶

Seed the far end state by stepping past a transition state.

Geodesically interpolates from reactant to ts, reads the direction the path is travelling in as it arrives at the saddle, and keeps stepping that way past it. What comes out is the structure on the other side: a product guess built from a reactant and a transition state alone.

This is the cheap counterpart to optimise_irc(), which answers the same question properly by following the true reaction coordinate downhill, at the cost of a converged saddle, a calculator and hundreds of gradients. Nothing here is evaluated: the seed is unrelaxed, and the way to turn it into an end state worth using is to relax it with optimise_geom(). Running prepare_neb() between the reactant and the relaxed seed, and checking the band comes back over a barrier near the transition state it started from, is what confirms the seed landed in the basin it was aimed at.

Nothing about this is specific to products. The reactant is only the end state you already have, so passing a product as reactant seeds the reactant instead, by stepping off the other side of the same saddle.

Parameters:

reactant (ase.Atoms) – End state to start from. Not modified.
ts (ase.Atoms) – Transition state, or a guess at one – it is used for its geometry and nothing else, so it need not be converged. Not modified. Must describe the same atoms, in the same order, as reactant.
n_images (int, optional) – Number of images in the geodesic interpolation from reactant to ts. More images resolve the arriving direction more finely.
push (float, optional) – How far past the transition state to step, as a multiple of the geodesic length from reactant to ts. The default of 1.0 steps as far beyond the saddle as the reactant sits before it, which for a near-symmetric reaction – a proton transfer, say – puts the seed roughly where the reactant would be reflected through the saddle. Raise it for a reaction whose product lies further out, lower it for a saddle that sits late along the path.
n_steps (int, optional) – Number of equal increments the push is taken in. Nothing but the clash check looks between them, so this sets both how finely a push that has to stop early can stop and how well that check sees: increments long enough to carry an atom clean through another one step over the clash without noticing it.
tangent_images (int, optional) – Number of trailing images the arrival direction is measured over. The default of 2 is the difference between the last two images; a larger chord is less sensitive to the geodesic wobbling near its endpoint, at the cost of averaging in curvature from further back.
weights ({"masses", None} or array_like, optional) – How the two images are superposed before their difference is taken, interpreted as in align_atom_sets(). The default weights by mass, which is the frame a reaction coordinate is conventionally defined in and which keeps the heavy atoms still while a hydrogen does the travelling. Weighting every atom equally instead spreads a proton’s motion back over the atoms it left, and extrapolating that pulls them apart. Pass None for equal weights, or an array to fix the frame on part of the structure.
clash_scale (float or None, optional) – Stop stepping before any two atoms come closer than this fraction of the sum of their covalent radii, and warn SeedWarning saying how far it got. The default of 0.7 sits well inside a normal bond, so a bond forming across the saddle does not trip it, while a step driving atoms through each other does. None steps the whole way regardless.
return_path (bool, optional) – Also return every structure the seed was built from.

Returns:

ase.Atoms or tuple – The seeded end state, carrying info["seeded"], whether both the checks below passed; info["seed_push"], info["seed_rmsd_reactant"] and info["seed_rmsd_ts"], how far it actually travelled, in Å; and info["seed_alignment"], how much the direction it was stepped along had to do with the direction the path travelled overall, as a cosine. With return_path it is instead (seed, path), where path is the geodesic from the reactant to the transition state followed by every extrapolated image, ending on the seed – a band to plot, or to hand to restart_neb(). It is n_images + n_steps long unless the push stopped early. The whole band is put in the frame of ts, so its image n_images - 1 is exactly the transition state as given, while its first image is the reactant up to the rigid drift the interpolation accumulated along the way.

Raises:

ValueError – If reactant and ts do not describe the same atoms, if any of the arguments are out of range, or if the two structures are so alike that the path between them gives no direction to extrapolate along.

Warns:

SeedWarning – If the push had to stop early to avoid a clash; if the seed did not end up further from the reactant than the transition state already was, which means the step went nowhere useful; or if the direction the path arrived in bore little relation to the direction it travelled overall, which means there was no reaction coordinate there to read.

Notes

Geodesic interpolation measures plain Cartesian distances, with no minimum-image convention and no knowledge of the cell, and hands its path back rigidly rotated. For a periodic system that rotation is meaningless, since the cell does not follow it – which is why the band is put back in the frame of ts before anything is measured off it, and why the seed comes out where the cell expects it. What that does not rescue is a reaction whose atoms cross a cell boundary: the interpolation never sees the periodic image, so there is nothing sensible to extrapolate and this is the wrong tool for it.

Examples

>>> seed = seed_product_from_ts(reactant, ts)
>>> product = optimise_geom(seed, calc)
>>> summarise_neb(optimise_neb(prepare_neb(reactant, product, calc)))
...

class SeedSummary(ts_energy: float, ts_fmax: float, minima: list[Atoms], connecting: tuple[int, int] | None, n_directions: int, n_relaxations: int, n_stalled: int, n_unconverged: int, n_failed: int)[source]¶

Bases: object

What rattling a transition state found on either side of it.

Returned by seed_minima_from_ts(). The numbers describing an individual minimum ride on that structure’s info, as everywhere else in this package; what is collected here is the search as a whole.

Variables:

ts_energy (float) – Potential energy of the transition state as handed over, in eV. Every barrier is measured down from it.
ts_fmax (float) – Maximum force on the transition state as handed over, in eV/Å. Free, since the energy above had to be evaluated anyway, and the one number that separates “the rattle failed” from “that was never a saddle”: a structure sitting at 0.3 eV/Å is somewhere on a slope, and the structures found below it are wherever the slope happened to run out.
minima (list of ase.Atoms) – The distinct minima the rattles relaxed into, lowest energy first. Empty when every relaxation stalled or failed. Each carries info["converged"] from its relaxation, plus info["seed_count"], how many relaxations landed in it; info["seed_barrier"], how far it sits below the transition state in eV; and info["seed_rmsd_ts"], how far it moved to get there in Å.
connecting (tuple of int or None) – Indices into minima of the pair the search most often landed either side of, or None when no mirrored pair ever bracketed the saddle – which is the headline failure, and what mirror=False guarantees.
n_directions (int) – Number of random directions drawn.
n_relaxations (int) – Number of relaxations run, which is twice n_directions when the directions were mirrored.
n_stalled (int) – Relaxations that came back onto the transition state instead of rolling off it, because the displacement was too small to escape it.
n_unconverged (int) – Relaxations that hit their step limit.
n_failed (int) – Relaxations the calculator raised on.

property barriers: ndarray¶

how far each minimum sits below the saddle, in eV.

Parallel to minima, and so also in ascending energy order, which makes this descending. A negative entry is a minimum that came back above the transition state, which seed_minima_from_ts() warns about.

Type:: numpy.ndarray

property connecting_minima: tuple[Atoms, Atoms] | None¶

the bracketing pair, ready for a band.

The two structures connecting names, or None when nothing bracketed the saddle. Handing these to prepare_neb() and checking the band comes back over a barrier near the transition state they were found from is what settles whether the rattle landed where it looked like it did.

Type:: tuple of ase.Atoms or None

seed_minima_from_ts(ts: ~ase.atoms.Atoms, calc: ~ase.calculators.calculator.Calculator, n_directions: int = 5, stdev: float = 0.1, mirror: bool = True, indices: int | ~collections.abc.Iterable[int] | None = None, fmax: float = 0.05, steps: int = 200, rmsd_tol: float = 0.3, energy_tol: float = 0.01, align: bool | None = None, seed: int | ~numpy.random._generator.Generator | None = 0, optimiser: ~collections.abc.Callable[[...], ~typing.Any] = <class 'ase.optimize.bfgs.BFGS'>, logfile: str | ~typing.TextIO | None = None, seed_traj: str = 'seed_minima.traj', keep_traj: bool = False, raise_on_unconverged: bool = False) → SeedSummary[source]¶

Rattle a transition state downhill to find the minima it connects.

Displaces the saddle at random, relaxes what comes out, and repeats, collecting the distinct basins the relaxations fall into. A saddle is downhill in one direction and uphill in every other, so a structure nudged off it rolls into one of the two states it connects; doing that from enough directions finds both.

Each direction is stepped both ways. From a saddle a displacement and its negative roll to opposite sides, so a mirrored pair brackets the reaction rather than sampling one basin twice, and the pair the search most often brackets is the answer this exists to give – reported in connecting, and available as two structures ready for a band through SeedSummary.connecting_minima.

This is the crude counterpart to optimise_irc(), which answers the same question properly by following the true reaction coordinate downhill, and needs Sella, a tightly converged saddle and hundreds of gradients to do it. Where seed_product_from_ts() builds one end state out of a saddle and the other end state, with no calculator at all, this needs no end state and finds both – for the price of 2 * n_directions full geometry relaxations, ten by default.

Nothing here checks that ts is a first-order saddle; that is what get_vibrations() is for, and it costs more than this whole search. ts_fmax on the summary is the cheap stand-in: a structure handed over at 0.3 eV/Å was never a stationary point, and whatever comes back is wherever the slope ran out.

Parameters:

ts (ase.Atoms) – Transition state, or a guess at one. Not modified; every relaxation runs on its own copy.
calc (ase.calculators.Calculator) – Calculator attached to the transition state and to every relaxation.
n_directions (int, optional) – Number of random directions to draw. Each costs two relaxations unless mirror says otherwise, so this is half the calculator budget, not all of it.
stdev (float, optional) – Standard deviation of the displacement drawn for each atom and each Cartesian component, in Å. The one knob that matters. Too small and every relaxation converges straight back onto the saddle, because the force there is already below any usable fmax; too large and the structure lands somewhere the saddle never connected to. The default of 0.1 Å clears a converged saddle reliably. Note that ASE’s own Atoms.rattle defaults to 0.001 Å, which would stall every time – the two are not meant to match.
mirror (bool, optional) – Relax each direction from both sides of the saddle. On by default, because a mirrored pair is the only thing here that brackets anything: with this off nothing can be shown to lie either side of the saddle, connecting is always None, and what is left is a way of enumerating the minima near a structure for half the cost.
indices (int or iterable of int, optional) – Atoms to displace, and to compare structures over. Every atom by default. Worth setting for a large system on both counts: an isotropic rattle spends most of its amplitude on atoms with nothing to do with the reaction, and an RMSD taken over everything dilutes as one over the square root of the atom count, until a real hop between basins measures smaller than the scatter a loose fmax leaves behind.
fmax (float, optional) – Maximum force criterion for the relaxations, in eV/Å. Loose by the standards of the rest of the package, because what is wanted from these is which basin, not the bottom of it – but see rmsd_tol.
steps (int, optional) – Maximum number of optimiser steps per relaxation. A structure starting a hair off a saddle that has not settled within a couple of hundred steps has gone somewhere unintended, and the warning saying so is more use than grinding on.
rmsd_tol (float, optional) – Largest RMSD in Å at which two relaxed structures still count as the same minimum. Coupled to fmax: a loose force criterion leaves structures scattered around the bottom of a basin, and this has to be wider than that scatter or one basin is reported as several. Tighten fmax before tightening this.
energy_tol (float, optional) – Largest energy difference in eV at which two relaxed structures still count as the same minimum, and the width of the band around the transition state’s own energy that counts as not having left it. Roughly kBT/2 at room temperature.
align (bool or None, optional) – Remove the best rigid rotation and translation before measuring any RMSD. None, the default, decides per structure: off when the cell is periodic or an atom is held by FixAtoms, on otherwise. A free molecule needs it, because a rattle imparts net translation and rotation that the relaxation will never remove – the force along a rigid mode is zero – so without it two structures in the same basin look different and the whole search fragments. A periodic cell or a pinned substrate must not have it: the frame is already fixed, and rotating a structure the cell does not follow is meaningless.
seed (int, numpy.random.Generator or None, optional) – Seed for the displacements, passed to numpy.random.default_rng(), so an integer, a generator or None all work. Defaults to 0: what this function finds is entirely a function of the draw, and a result that cannot be reproduced is hard to argue with.
optimiser (callable, optional) – ASE optimiser class used for the relaxations, as in optimise_geom().
logfile (str, file object or None, optional) – Where the optimisers write their per-step tables. None, the default, silences them – deliberately unlike the rest of the package, because ten interleaved tables with nothing marking where one ends and the next begins are worse than no tables at all. Pass '-' for stdout.
seed_traj (str, optional) – Base trajectory filename. Each relaxation gets its own name derived from it, seed_minima_0_plus.traj and so on, so that keeping them keeps all of them rather than only the last.
keep_traj (bool, optional) – Keep those trajectories instead of deleting them. Off by default. A relaxation the calculator raised on leaves its trajectory behind either way, since that is the evidence for what went wrong.
raise_on_unconverged (bool, optional) – Raise ConvergenceError instead of warning when any relaxation hit steps without reaching fmax. Checked once, after every relaxation has run, so that turning it on does not throw away the work already done – the same order optimise_irc() runs its two halves in.

Returns:

SeedSummary – The distinct minima found, which pair of them the search bracketed the saddle with, and what was discarded getting there.

Raises:

TypeError – If indices is neither an integer nor an iterable of integers.
ValueError – If any argument is out of range, or indices is empty or repeats.
IndexError – If indices names an atom outside ts. Negative indices are not accepted; use len(ts) - 1 for the last atom.
RuntimeError – If the calculator raised on every single relaxation, chained to the last exception it raised.
ConvergenceError – If any relaxation hit steps without reaching fmax and raise_on_unconverged is True.

Warns:

SeedWarning – If any relaxation came back onto the transition state instead of rolling off it; if no mirrored pair bracketed the saddle, so there is no answer to the question this was run to ask; if the calculator raised on some of the relaxations; or if a minimum came back above the transition state, which means either that the input was not a maximum along that direction or that the relaxation stopped short.

Notes

Everything measured here goes through atom_set_rmsd(), which matches atoms by index, applies no minimum-image convention and does not consider permutations. So two minima that differ only by a lattice translation are reported as distinct, as are two related by swapping identical atoms; and a mechanism in which two atoms exchange places reads as an enormous displacement. A reaction that crosses a cell boundary is the wrong shape for this function.

The displacement is not mass-weighted, unlike the direction seed_product_from_ts() extrapolates along, so a hydrogen and a tungsten are moved by the same amount and a proton transfer gets less of the amplitude than it wants. indices is the way to put it back.

Growing the system does not call for a larger stdev. The component of an isotropic displacement along any one direction is distributed as N(0, stdev) however many atoms there are; what grows is the total, stdev * sqrt(3 * n_atoms), which is strain dumped into every mode that is not the reaction coordinate. That costs optimiser steps, so steps is the argument a large system runs out of, not stdev.

Examples

>>> summary = seed_minima_from_ts(ts, calc)
>>> print(summary)
>>> reactant, product = summary.connecting_minima
>>> summarise_neb(optimise_neb(prepare_neb(reactant, product, calc)))
...