Introduction to atomman: Minimum energy paths

Lucas M. Hale, lucas.hale@nist.gov, Materials Science and Engineering Division, NIST.

1. Introduction

New version 1.3.7

This Notebook outlines the mep subpackage of atomman that provides tools for identifying minimum energy pathways. These tools were designed primarily for finding relaxed energy paths along 2D energy maps, but may also work for more complex problems.

NOTE #1 This documentation uses the term “energy” for the output of the evaluation function and the term “force” for the derivative of the evaluation function with respect to the coordinates. This does not mean that the evaluation function needs to be in units of energy, just that it is a function analogous to energy.

Note #2 The atomman.mep module is relatively new and only contains a few options at the moment. The code has been designed with modularity in mind, so additional methods can hopefully be added in the future without requiring too many changes to the existing design.

2. Theory and Methodology

Minimum energy paths are transition paths from one (meta)stable state to another in which each point along the path is located at a minimum for all directions tangent to the path’s line. Another way of stating this is that every point along the path has a gradient of zero for all directions except along the path’s line. The minimum energy path for a transition is important as it marks the most likely states that the system will occupy during a transition and it provides an accurate means of identifying the transition barrier.

The atomman.mep module contains Path classes. Each Path class provides - A representation of a path along an energy surface as a series of coordinates, - Methods that allow for the integration of the path towards the minimum energy path using a particular relaxation algorithm, and - Tools for viewing the path, as well as the energies and forces along the path.

3. BasePath class and common setup

3.1. Path initialization

All Path classes are children of BasePath, which defines the commonly shared methods and attributes.

A Path can be initialized with the following parameters

coord (array-like object) The list of coordinates associated with the points along the path.
energyfxn (function) The function that evaluates the energy associated with the different point coordinates.
gradientfxn (str or function, optional) The function to use to estimate the gradient of the energy. A str value of ‘central_difference’ or ‘cdiff’ (default) will use atomman.mep.gradient.central_difference.
gradientkwargs (dict, optional) The keyword arguments (i.e. settings) to use wit the gradientfxn. Default is an empty dictionary, i.e. default settings of gradientfxn.
integratorfxn (str or function, optional) The function to use to integrate relaxation steps. A str value of euler will use atomman.mep.integrator.euler, while a str value of ‘rungekutta’ or ‘rk’ (default) will use atomman.mep.integrator.rungekutta.

3.2. gradientfxn, gradientkwargs and integratorfxn

The gradientfxn and gradientkwargs define the function to use for evaluating the gradient of coord. The atomman.mep.gradient module collects the following gradient methods - atomman.mep.gradient.central_difference(energyfxn, coord, shift=1e-5) computes the gradient of coord using central difference. For each dimension, the energy is evaluated at coord\(\pm\)shift and the gradient computed as the difference of the shifted energies divided by 2*shift.

The integratorfxn defines the integration scheme to use to iterate coord forward by timestep in the direction given by the ratefxn to produce a new set of coord. The atomman.mep.integrator module collects the following integration methods - atomman.mep.integrator.euler(ratefxn, coord, timestep, **kwargs) uses the Euler integration method - newcoord = coord + timestep * ratefxn(kwargs) - **atomman.mep.integrator.rungekutta(ratefxn, coord, timestep, **kwargs) uses the Runge-Kutta integration method, which uses a higher order estimate based on four shifts of coord. - k1 = timestep * ratefxn(coord, kwargs) - k2 = timestep * ratefxn(coord - 0.5 * k1, kwargs) - k3 = timestep * ratefxn(coord - 0.5 * k2, **kwargs) - k4 = timestep * ratefxn(coord - k3, **kwargs) - newcoord = coord + k1 / 6 + k2 / 3 + k3 / 3 + k4 / 6

3.3. Common Path attributes

coord (numpy.NDArray) is the array of coordinates along the path.
arccoord (numpy.NDArray) computes the arc length coordinates for coord. This is the cumulative sum of the radial distances between each set of neighboring points along the path.
energyfxn (function) is the function to evaluate the energy associated with the coord points.
gradientfxn (function) is the function for computing the energy gradient.
gradientkwargs (dict) is a dict containing additional terms to pass to the gradientfxn.
integratorfxn (function) is the function for iterating coord forward by a timestep.
unittangent (numpy.NDArray)computes the unit tangent vector along the path at each point. Note that this is not computed the same way for all subclasses.
force (numpy.NDArray) computes the force associated with moving along the path at each coord.

3.4. Common Path methods

energy(coord=None) computes the energy using energyfxn for the set/given coord points.
grad_energy(coord=None) computes the gradient of the energy using energyfxn and gradientfxn for the set/given coord points.
step(*args, **kwargs) returns a new Path object with coord evolved forward by one timestep.
relax(*args, **kwargs) performs multiple steps until a tolerance or a max number of steps is reached. Returns a new Path object with the final coord.
plot_energy(energy_unit=None, length_unit=None, ax=None, **kwargs) plots the energy vs. arccoord for all points along the path. ax is an optional existing matplotlib.pyplot.axis object, allowing for subplots to be constructed. All additional kwargs are passed on to matplotlib.pyplot.figure(). Note that energy_unit should correspond to the units returned by energyfxn, which may not be units of energy per se.
plot_force(force_unit=None, length_unit=None, ax=None, **kwargs) plots the force vs. arccoord for all points along the path. ax is an optional existing matplotlib.pyplot.axis object, allowing for subplots to be constructed. All additional kwargs are passed on to matplotlib.pyplot.figure(). Note that force_unit should correspond to the units returned by energyfxn divided by length, which may not be units of force per se.

4. ISMPath

The ISMPath class uses the improved string method (ISM) [link?] to compute the minimum energy pathways.

4.1. Improved string method theory

In the improved string method, each relaxation step involves

An initial path guess is given where all coordinates along the path are equally spaced according to the arc length coordinates.
Each point is integrated by one time step along the negative gradient direction.
The coordinates of the integrated points are then shifted such that they are once again equally spaced according to the arc length coordinates. This is done by fitting a cubic spline interpolation relating arc length coordinates to the full coordinates of the relaxed points, then interpolating full coordinates associated with equally spaced arc length coordinates.

When well behaved, step 2 relaxes the coordinates according to the gradient, and step 3 pulls the points back up to remain along the path.

After relaxing, a “climbing” algorithm can be used to better identify the critical barrier point along the path. This follows the same basic concepts of the relaxation steps except that 1. One or more climbing points are identified as having local maximum values after the relaxation steps.

2. The rate function for integrating the climbing points is modified from the negative gradient by reversing the gradient only along the unit tangent direction of the line path. This means that the climbing points will climb up the unit tangent direction rather than relax down. 3. When applying the arc length adjustments, the full line path is broken into segments according to the climbing points. The points within each segment are then adjusted to be equally arc length spaced with the end/climbing points not adjusted.

4.2. ISMPath specific methods and attributes

default_timestep (float) is the default relaxation timestep calculated as 0.05 * min(0.2, N^-1), where N is the number of coord points.
default_tolerance (float) is the default relaxation tolerance calculated as max(N^-4, 1e-10), where N is the number of coord points.
unittangent (numpy.NDArray) are the unit tangent vectors along the path line at each coord point. For the ISM method, the tangent vector for the starting and ending points are calculated as the difference vector between those points and their neighbors, while the tangent vectors for all intermediate points are calculated as the difference vectors between their two neighboring points.
interpolate_point(arccoord) constructs a cubic spline interpolation based on the current coord and corresponding arccoord values. The given arccoord parameter values are then used to interpolate new coord values. The method returns a new ISMPath object with the interpolated coord values.

4.2.1. ISMPath.step()

Performs a single string relaxation step.

Parameters

timestep (float, optional) The size of the timestep to use. Will use the path’s default timestep if not given.
climbindex (int, list or None, optional) Indicates the indices of the path points to apply the climb algorithm to. If None (default), no climbing will be performed.

Returns

newpath (Path) A Path with coordinates evolved forward by one timestep.

4.2.2. ISMPath.relax()

Perform multiple relaxation and/or climb steps until either the maximum coordinate displacement per step drops below a tolerance or the maximum number of steps is reached.