tools_geodesic — geodesic interpolation

Geodesic interpolation: a smooth first guess at a reaction path.

The band reactiontools.tools_reaction.prepare_neb() starts from, and what reactiontools.tools_reaction.quick_guess_path() and reactiontools.tools_reaction.quick_guess_ts() hand back. Given two end states – or more, when an intermediate is already known – geodesic_interpolate() returns a path of images that is short in a metric of scaled inter-atomic distances rather than in Cartesian space, which is what keeps atoms from passing through each other on the way across.

It runs in two stages. redistribute() builds a raw path with the right number of images, bisecting the largest gap until there are enough of them and dropping the least costly image while there are too many. Geodesic then smooths that path by minimising its length in the internal metric, holding the two end points fixed. Up to SWEEP_ABOVE_N_ATOMS atoms the whole path is optimised at once; past that the images are smoothed one at a time, sweeping back and forth, as SciPy’s least-squares optimisers slow down badly with problem size.

The coordinates are redundant internals: every atom pair within a cut-off, plus every pair within three bonds of each other however far apart they are, each distance passed through a scaling function that decays with distance – morse_scaler() by default. The optimisation is carried out in Cartesians all the same, because an arbitrary set of redundant internal coordinate values need not correspond to any real geometry.

This module is a leaf of the package: it reaches ase, numpy and scipy, and nothing else here beyond reactiontools.tools_io.write_xyz_frame(). Its alignment helpers – align_geom(), align_path() and align_path_to() – exist for that reason. They are the bare Kabsch fit the interpolation needs between adjacent images; reactiontools.tools_geometry.align_atom_sets() is the fuller one, with weights, subset selection and validation, for superposing two structures you already have.

The interpolation itself is not periodic. The internal coordinates are plain inter-atomic distances with no minimum image convention, so a bond that crosses a cell boundary is not handled. A periodic path is nonetheless moved back onto the frame of reference it arrived in, since a cell only says where the atoms are in the frame it was given with.

Derived from geodesic-interpolate by Xiaolei Zhu, MIT licensed; the copyright notice travels with this package in LICENSE. It came by way of the fork at https://github.com/LouieSlocombe/geodesic_interpolate, commit 1e37b32aab77bcab0273e806b5a0a35df81160dd of 20 August 2026, which is the tree to diff this file against and the one to re-sync from. Cite zhu2019geodesic when you use it – see CITATIONS.bib.

SWEEP_ABOVE_N_ATOMS = 35

System size, in atoms, above which geodesic_interpolate() sweeps the images one at a time rather than optimising the whole path at once. SciPy’s least-squares optimisers slow down faster than linearly in the number of degrees of freedom, so past this the sweep is the cheaper of the two.

class Geodesic(atoms: list[str], path: ndarray | list[ndarray], scaler: float | Callable[[ndarray], tuple[ndarray, ndarray]] = 1.7, threshold: float = 3.0, min_neighbors: int = 4, friction: float = 0.001, rng: Generator | None = None)

Bases: object

Optimiser that finds a geodesic in redundant internal coordinates.

The heart of it is the path length in the internal metric: the length of each segment is measured through its midpoint, so the target function is a sum of squared displacements and can be handed straight to a least-squares optimiser.

smooth() optimises every image at once and sweep() takes them one at a time; both leave the result in path. The two end points are held fixed either way.

Variables:

• path (numpy.ndarray) – The current geometries, of shape (n_images, n_atoms, 3). Updated in place as the optimisation proceeds.

• n_atoms (n_images,) – Shape of the path, as images and atoms.

• rij_list (list of tuple of int) – The atom pairs making up the internal coordinates.

• n_rij (int) – How many of them there are.

• re (numpy.ndarray) – Equilibrium distance for each of those pairs.

• scaler (callable) – The scaling function setting the metric.

• friction (float) – Default weight of the friction term regularising the step size.

• length (float) – Path length in the internal metric, over the section last measured. Only set once the displacements have been computed.

• optimality (float) – Infinity norm of the gradient of the length, which is what convergence is judged on. Only set once the target function has been evaluated.

update_intc() → None

Fill in any internal coordinates currently marked unknown.

Missing entries are flagged with None in the internal storage; this finds them, works out the midpoint geometries where needed, and evaluates the coordinates and their gradients. Everything already known is left alone, so nothing is evaluated twice.

The derivatives are kept sparse, since each scaled distance only depends on the six Cartesian components of its own two atoms.

update_geometry(geom: ndarray, start: int, end: int) → bool

Move a segment of the path, invalidating what depended on it.

The internal coordinates, derivatives and midpoints of the affected images are reset to unknown so update_intc() recomputes them. Note that moving images start:end also invalidates the midpoint just before start, if there is one.

Parameters:

• geom (numpy.ndarray) – New Cartesian coordinates for the segment, flattened or otherwise.

• start (int) – The section of the path being replaced.

• end (int) – The section of the path being replaced.

Returns:

bool – True if the geometry actually changed, False if geom was already the current geometry and nothing needed invalidating.

compute_displacements(start: int = 1, end: int = -1, dx: ndarray | None = None, friction: float = 0.001) → None

Compute the displacements along a section of the path, and its length.

Each segment is split at its midpoint and measured in two halves, which is what makes the length a sum of squares. The friction term is appended to the same vector so the least-squares optimiser sees it as extra residuals.

Sets length and displacements.

Parameters:

• start (int, optional) – Section of the path to measure. A negative end counts back from the last image. start must be at least 1: the end points of the path are fixed, and each segment is measured against the image before it.

• end (int, optional) – Section of the path to measure. A negative end counts back from the last image. start must be at least 1: the end points of the path are fixed, and each segment is measured against the image before it.

• dx (numpy.ndarray, optional) – Displacement of the segment from its reference geometry. When given, it enters the target function scaled by friction.

• friction (float, optional) – Weight of the friction term.

Raises:

ValueError – If the section does not lie between the two fixed end points.

compute_disp_grad(start: int, end: int, friction: float = 0.001) → None

Differentiate the displacement vectors with respect to Cartesians.

Moving one image changes the two half-segments on either side of it directly, and the two midpoints it shares with its neighbours at half the rate, which is where the factors of a half come from. The friction residuals contribute a diagonal block at the bottom of the matrix.

An image only appears in the two half-segments it touches, so the matrix is block-bidiagonal, and each block is itself sparse because a distance only depends on its own two atoms. It is assembled sparse and left that way: for a few dozen atoms fewer than one entry in a hundred is non-zero, and handing a dense array to scipy.optimize.least_squares() makes it factorise the whole thing on every iteration.

Sets grad.

Parameters:

• start (int) – Section of the path being differentiated.

• end (int) – Section of the path being differentiated.

• friction (float, optional) – Weight of the friction term, matching compute_displacements().

compute_target_func(geom: ndarray | None = None, start: int = 1, end: int = -1, x0: ndarray | None = None, friction: float = 0.001) → None

Compute the vectorised target function the minimisation uses.

Sets displacements, grad and optimality.

Parameters:

• geom (numpy.ndarray, optional) – Geometry to evaluate at. If it matches the geometry already stored for this segment, the whole evaluation is skipped.

• start (int, optional) – Section of the path being optimised.

• end (int, optional) – Section of the path being optimised.

• x0 (numpy.ndarray, optional) – Reference geometry the friction term pulls back towards. Defaults to no pull at all.

• friction (float, optional) – Weight of the friction term.

target_func(geom: ndarray, **kwargs: Any) → ndarray

Residuals for the optimiser.

Wraps compute_target_func(), which skips the work if the geometry has not moved since the last call.

Parameters:

• geom (numpy.ndarray) – Geometry of the segment, flattened as SciPy hands it over.

• **kwargs – Passed straight through to compute_target_func().

Returns:

numpy.ndarray – The displacement vectors of the segment, with the friction residuals appended.

target_deriv(geom: ndarray, **kwargs: Any) → csr_matrix

Jacobian for the optimiser.

Wraps compute_target_func(), which skips the work if the geometry has not moved since the last call. Paired with target_func(), this means each geometry is only ever evaluated once even though SciPy asks for value and Jacobian separately.

Parameters:

• geom (numpy.ndarray) – Geometry of the segment, flattened as SciPy hands it over.

• **kwargs – Passed straight through to compute_target_func().

Returns:

scipy.sparse.csr_matrix – The derivatives of those displacement vectors, held sparse.

smooth(tol: float = 0.001, max_iter: int = 50, start: int = 1, end: int = -1, friction: float | None = None, xref: ndarray | None = None) → ndarray

Minimise the path length as one function of all the coordinates.

In principle this is very efficient, but it can get costly for large systems with many images, in which case sweep() is the better option.

Parameters:

• tol (float, optional) – Convergence tolerance on the optimality, i.e. the uniform gradient of the target function.

• max_iter (int, optional) – Ceiling on the number of function evaluations, passed to SciPy as max_nfev.

• start (int, optional) – Section of the path to optimise.

• end (int, optional) – Section of the path to optimise.

• friction (float, optional) – Weight of the friction term. Defaults to the value given to the constructor.

• xref (numpy.ndarray, optional) – Reference geometry for the friction term. Defaults to the starting geometry of the segment.

Returns:

numpy.ndarray – The optimised path. This is also stored in path.

sweep(tol: float = 0.001, max_iter: int = 50, micro_iter: int = 20, start: int = 1, end: int = -1) → ndarray

Minimise the path length one image at a time, sweeping back and forth.

Less efficient per iteration than smooth(), but it scales far more kindly with system size given how slow SciPy’s optimisers get on large problems. It also allows finer control, and makes it cheap to skip images that are already close to optimal.

Parameters:

• tol (float, optional) – Convergence tolerance on the optimality, i.e. the uniform gradient of the target function.

• max_iter (int, optional) – Maximum number of sweeps through the path.

• micro_iter (int, optional) – Number of micro-iterations spent optimising each image.

• start (int, optional) – Section of the path to optimise.

• end (int, optional) – Section of the path to optimise.

Returns:

numpy.ndarray – The optimised path. This is also stored in path.

align_geom(ref_geom: ndarray, geom: ndarray) → tuple[float, ndarray]

Find the rigid-body motion that maximally overlaps two geometries.

Implemented with the Kabsch algorithm: both geometries are centred, the SVD of their covariance matrix gives the optimal rotation, and the result is moved back onto the centre of the reference.

Every atom counts equally and all of them are used. See reactiontools.tools_geometry.align_atom_sets() for the weighted, subset-aware superposition of two structures.

Parameters:

• ref_geom (numpy.ndarray) – The reference geometry to rotate towards, of shape (n_atoms, 3).

• geom (numpy.ndarray) – The geometry to be rotated and shifted, of the same shape.

Returns:

• rmsd (float) – Root-mean-squared difference between the rotated geometry and the reference.

• aligned_geom (numpy.ndarray) – The rotated geometry that maximally overlaps the reference.

align_path(path: ndarray | list[ndarray]) → tuple[float, ndarray]

Rotate and translate the images of a path to minimise movement along it.

The first image is shifted so its geometric centre sits at the origin, and each following image is aligned against the one before it, which leaves every image centred on the origin as well.

Parameters:

path (array_like) – Sequence of geometries, of shape (n_images, n_atoms, 3).

Returns:

• max_rmsd (float) – Largest RMSD between any two adjacent images after alignment.

• path (numpy.ndarray) – The aligned path. This is a copy, so the input is left untouched.

align_path_to(ref_geom: ndarray, path: ndarray) → ndarray

Move a whole path rigidly so its first image sits on a reference.

align_path() leaves the path centred on the origin and rotated onto its own first image, which is a frame of the optimiser’s choosing rather than the one the caller handed in. That is immaterial for an isolated molecule, but a unit cell only describes where the atoms are in the frame it came with, so a periodic path has to be put back before its cell means anything.

One rotation and translation, taken from the first image, is applied to every image, so the shape of the path and its length are untouched.

Parameters:

• ref_geom (numpy.ndarray) – Geometry to put the first image back onto, of shape (n_atoms, 3).

• path (numpy.ndarray) – The path to move, of shape (n_images, n_atoms, 3).

Returns:

numpy.ndarray – The path in the reference geometry’s frame. This is a copy.

compute_rij(geom: ndarray, rij_list: list[tuple[int, int]]) → tuple[ndarray, ndarray]

Calculate inter-atomic distances and their Cartesian derivatives.

Parameters:

• geom (numpy.ndarray) – Cartesian geometry of all the atoms, of shape (n_atoms, 3).

• rij_list (list of tuple of int) – Indices of the atom pairs to evaluate.

Returns:

• rij (numpy.ndarray) – The distance for each pair.

• b_mat (numpy.ndarray) – Wilson B matrix, of shape (n_pairs, n_atoms, 3), holding the Cartesian gradient of every distance.

compute_wij(geom: ndarray, rij_list: list[tuple[int, int]], func: Callable[[ndarray], tuple[ndarray, ndarray]], sparse: bool = False) → tuple[ndarray, ndarray | csr_matrix]

Calculate scaled distances and their Cartesian derivatives.

Same as compute_rij(), except each distance is passed through a scaling function which sets the metric of the internal coordinates.

Parameters:

• geom (numpy.ndarray) – Cartesian geometry of all the atoms. Flattened input is accepted and reshaped to (n_atoms, 3).

• rij_list (list of tuple of int) – Indices of the atom pairs to evaluate.

• func (callable) – Scaling function returning both the scaled value and its derivative with respect to the raw distance. Must broadcast over arrays.

• sparse (bool, optional) – Return the B matrix as a sparse matrix rather than a dense array. Only six of its entries per row are ever non-zero, so this is what the optimisers are given: it saves building the dense array, and lets SciPy solve the least-squares steps iteratively instead of by dense factorisation.

Returns:

• wij (numpy.ndarray) – The scaled distance for each pair.

• b_mat (numpy.ndarray or scipy.sparse.csr_matrix) – Cartesian gradients of the scaled distances, of shape (n_pairs, 3 * n_atoms), with the atom and component axes flattened together so scipy.optimize can use it directly. A dense array by default, or a sparse matrix when sparse is set.

elu_scaler(re: float | ndarray = 2.0, alpha: float = 2.0, beta: float = 0.01) → Callable[[ndarray], tuple[ndarray, ndarray]]

Build a scaling function with an exponential tail and a linear core.

Shaped like an ELU: beyond re the scaled distance decays exponentially, as in morse_scaler(), while below re it continues linearly along the tangent at re instead of blowing up. A beta * re / r tail is added for the same reason as in the Morse case.

Parameters:

• re (float or numpy.ndarray, optional) – Distance at which the behaviour switches from linear to exponential.

• alpha (float, optional) – Decay constant of the exponential, which also sets the linear slope.

• beta (float, optional) – Weight of the long-range tail.

Returns:

callable – A function taking an array of distances to the scaled distances and their derivatives with respect to the raw distance, in the form compute_wij() expects.

from_ase_atoms(atoms: list[Atoms]) → tuple[list[str], list[ndarray]]

Split a list of ASE Atoms objects into symbols and coordinates.

Every frame must hold the same atoms in the same order, since the symbols are taken from the first one and the interpolation matches the frames up atom by atom. Two end states built separately are the usual way to get that wrong, so it is checked here rather than assumed: the coordinates alone carry no record of which atom is which, and a path interpolated between mismatched frames is silently a path to the wrong structure.

Parameters:

atoms (list of ase.Atoms) – Frames of the path.

Returns:

• atom_names (list of str) – Element symbols of all the atoms.

• coords (list of numpy.ndarray) – Cartesian coordinates for every frame.

Raises:

ValueError – If the frames do not all describe the same atoms in the same order.

geodesic_interpolate(atoms: list[Atoms] | str | PathLike, n_images: int = 17, output: str | PathLike = 'interpolated.xyz', tol: float = 0.002, max_iter: int = 50, micro_iter: int = 20, scaling: float = 1.7, friction: float = 0.01, dist_cutoff: float = 3.0, seed: int = 42) → list[Atoms] | None

Interpolate a reaction path between two or more geometries.

Runs the two stages in turn: redistribute() builds a raw path with the requested number of images, then Geodesic smooths it into a geodesic under the internal coordinate metric.

Which smoothing is used depends on the size of the system. Up to SWEEP_ABOVE_N_ATOMS atoms the whole path is optimised at once; beyond that the images are smoothed one at a time, sweeping back and forth along the path, because SciPy’s optimisers slow down badly as the problem grows.

Input and output mirror each other. Given ASE Atoms objects the interpolated path comes back as Atoms objects; given a filename it is written to output instead.

Given Atoms objects, everything the interpolation does not itself touch is taken from the first frame and carried onto every image: the unit cell, the boundary conditions, constraints, tags and so on. A periodic path is also moved back onto the frame of reference of the input, since the optimisation otherwise leaves it centred and rotated into a frame of its own, which would put the atoms in the wrong place relative to the cell. Note that the interpolation itself is not periodic: the internal coordinates are plain inter-atomic distances with no minimum image convention, so a bond that crosses a cell boundary is not handled.

Parameters:

• atoms (list of ase.Atoms, or str or os.PathLike) – Either a list of ASE Atoms objects, or the name of an XYZ file holding the end points. Only the first and last geometries need be meaningful, but intermediate ones are used if present.

• n_images (int, optional) – Number of images in the interpolated path.

• output (str or os.PathLike, optional) – XYZ file to write to. Only used when atoms is a filename.

• tol (float, optional) – Convergence tolerance for the smoothing, judged on the uniform gradient of the path length. The raw path is built at five times this, being only a starting guess for the smoothing.

• max_iter (int, optional) – Ceiling on the smoothing, counting function evaluations when the whole path is optimised at once and sweeps when sweeping. It is a ceiling rather than a target, as the optimisation stops as soon as it meets tol: the bundled test systems all converge somewhere between 10 and 35 at the default seed, though the count moves around with the seed. Setting it too low silently truncates the path part-way through the descent.

• micro_iter (int, optional) – Micro-iterations spent on each image. Only used when sweeping, so it has no effect on systems of SWEEP_ABOVE_N_ATOMS atoms or fewer.

• scaling (float, optional) – Alpha parameter of the Morse scaler setting the coordinate metric.

• friction (float, optional) – Weight of the friction term regularising the optimisation step size.

• dist_cutoff (float, optional) – Distance cut-off for building the internal coordinates. Atoms within three bonds of each other are included whatever their distance.

• seed (int, optional) – Seed for the random nudges and image sampling, so runs reproduce. The bisection is stochastic, and without a fixed seed larger systems will not give the same path twice. The seed goes to a numpy.random.Generator used only here, so the caller’s own random state is left alone.

Returns:

list of ase.Atoms or None – The interpolated path when atoms was a list of Atoms objects, otherwise None, with the path written to output.

Raises:

• TypeError – If atoms is neither a list of ASE Atoms nor a filename.

• ValueError – If fewer than two geometries are supplied.

get_bond_list(geom: ndarray | list[ndarray], atoms: list[str] | None = None, threshold: float = 4.0, min_neighbors: int = 4, snapshots: int = 30, bond_threshold: float = 1.8, enforce: Iterable[tuple[int, int]] = (), rng: Generator | None = None) → tuple[list[tuple[int, int]], ndarray]

Get the list of atom pairs that define the internal coordinate system.

Samples images from the path and collects every pair of atoms that comes within threshold of each other in any of them. Pairs linked by three or fewer bonds are always included, as are the pairs given in enforce, and further pairs are added for any atom left with fewer than min_neighbors neighbours.

Parameters:

• geom (array_like) – A single geometry or a path of them. Anything that is not already (n_images, n_atoms, 3) is promoted to that shape.

• atoms (list of str, optional) – Element symbols, used to look up covalent radii. When omitted, every pair is given a nominal equilibrium distance of 2.0 angstrom instead.

• threshold (float, optional) – Distance cut-off for including a pair in the coordinate system.

• min_neighbors (int, optional) – Minimum number of pairs each atom must take part in. Atoms below this get their nearest neighbours added.

• snapshots (int, optional) – Maximum number of images to sample. Keeps the cost down when the path is long and the atoms numerous.

• bond_threshold (float, optional) – Distance below which two atoms count as bonded, used to work out which pairs are within three bonds of each other.

• enforce (iterable of tuple of int, optional) – Pairs to include regardless of how far apart the atoms are.

• rng (numpy.random.Generator, optional) – Random source for choosing which images to sample. Defaults to a fresh unseeded generator. Pass a seeded one to keep the choice reproducible without disturbing anyone else’s random state.

Returns:

• rij_list (list of tuple of int) – The (i, j) atom pairs making up the coordinates.

• re (numpy.ndarray) – Equilibrium distance for each pair, taken as the sum of the two covalent radii.

morse_scaler(re: float | ndarray = 1.5, alpha: float = 1.7, beta: float = 0.01) → Callable[[ndarray], tuple[ndarray, ndarray]]

Build a scaling function based on a Morse potential.

This is the metric geodesic interpolation uses by default. A small beta / r tail is added so the coordinate keeps responding to atoms that are far apart, where the exponential has already decayed away.

Parameters:

• re (float or numpy.ndarray, optional) – Equilibrium distance. Usually the per-pair array returned by get_bond_list() rather than a single number.

• alpha (float, optional) – Decay constant of the exponential. Larger values are more localised, which tracks a sharp energy landscape better, while smaller values have longer range and give smoother paths from few images.

• beta (float, optional) – Weight of the long-range tail.

Returns:

callable – A function taking an array of distances to the scaled distances and their derivatives with respect to the raw distance, in the form compute_wij() expects.

read_xyz(filename: str | PathLike) → tuple[list[str], list[ndarray]]

Read an XYZ file and return the atom names and coordinates.

Parameters:

filename (str or os.PathLike) – Name of the XYZ data file. It may hold any number of frames, and blank lines between them are ignored.

Returns:

• atom_names (list of str) – Element symbols of all the atoms, read from the last frame.

• coords (list of numpy.ndarray) – Cartesian coordinates for every frame.

Raises:

ValueError – If the file is empty or does not parse as XYZ.

redistribute(atoms: list[str], geoms: list[ndarray], n_images: int, tol: float = 0.01, rng: Generator | None = None) → list[ndarray]

Add or remove images so the path has the requested number of them.

If there are too few, new points are added by bisecting the largest RMSD gap. If there are too many, images are dropped one at a time, each time choosing the one whose removal leaves the shortest merged segment.

Parameters:

• atoms (list of str) – Atom symbols, used to look up covalent radii.

• geoms (list of numpy.ndarray) – Geometries of the original path.

• n_images (int) – The desired number of images.

• tol (float, optional) – Convergence tolerance for the bisection.

• rng (numpy.random.Generator, optional) – Random source for the bisection, which is stochastic. Defaults to a fresh unseeded generator.

Returns:

list of numpy.ndarray – An aligned path with the correct number of images.

to_ase_atoms(atoms: list[str], coords: ndarray | list[ndarray], template: Atoms | None = None) → list[Atoms]

Rebuild a list of ASE Atoms objects from symbols and coordinates.

Parameters:

• atoms (list of str) – Element symbols of all the atoms.

• coords (array_like) – Cartesian coordinates, of shape (n_images, n_atoms, 3). A single frame of shape (n_atoms, 3) is accepted too.

• template (ase.Atoms, optional) – Atoms object to take everything other than the positions from, so the unit cell, boundary conditions, constraints, tags and the rest survive a round trip through this module. Each frame is a copy of the template with its positions replaced; as with ase.Atoms.copy(), the calculator is not carried over, since its results belong to the geometry it was attached to. Any constraints come along for later use but are not applied to the coordinates, which are always exactly the ones given. Without a template the frames are built from the symbols alone and carry no cell.

Returns:

list of ase.Atoms – One Atoms object per frame.

Raises:

ValueError – If template describes a different set of atoms from atoms.

write_xyz(filename: str | PathLike, atoms: list[str], coords: ndarray | list[ndarray]) → None

Write atom names and coordinate data to a multi-frame XYZ file.

Parameters:

• filename (str or os.PathLike) – Name of the XYZ data file to write.

• atoms (list of str) – Element symbol for each atom.

• coords (array_like) – Coordinates, of shape (n_images, n_atoms, 3). A single frame of shape (n_atoms, 3) is accepted too.