Nudged elastic band¶
The Nudged Elastic Band method is a technique for finding transition paths (and corresponding energy barriers) between given initial and final states. The method involves constructing a “chain” of “replicas” or “images” of the system and relaxing them in a certain way.
Relevant literature References:
H. Jonsson, G. Mills, and K. W. Jacobsen, in ‘Classical and Quantum Dynamics in Condensed Phase Systems’, edited by B. J. Berne, G. Cicotti, and D. F. Coker, World Scientific, 1998 [standard formulation]
‘Improved Tangent Estimate in the NEB method for Finding Minimum Energy Paths and Saddle Points’, G. Henkelman and H. Jonsson, J. Chem. Phys. 113, 9978 (2000) [improved tangent estimates]
‘A Climbing-Image NEB Method for Finding Saddle Points and Minimum Energy Paths’, G. Henkelman, B. P. Uberuaga and H. Jonsson, J. Chem. Phys. 113, 9901 (2000)
‘Improved initial guess for minimum energy path calculations.’, S. Smidstrup, A. Pedersen, K. Stokbro and H. Jonsson, J. Chem. Phys. 140, 214106 (2014)
‘Scaled and Dynamic Optimizations of Nudged Elastic Bands’, P. Lindgren, G. Kastlunger and A. A. Peterson, J. Chem. Theory Comput. 15, 11, 5787-5793 (2019)
The NEB class¶
This module defines one class:
- class ase.mep.neb.NEB(images, k=0.1, climb=False, parallel=False, remove_rotation_and_translation=False, world=None, method='aseneb', allow_shared_calculator=False, precon=None, **kwargs)¶
Nudged elastic band.
G. Henkelman and H. Jonsson, Chem. Phys, 113, 9978 (2000). doi: 10.1063/1.1323224
G. Henkelman, B. P. Uberuaga, and H. Jonsson, Chem. Phys, 113, 9901 (2000). doi: 10.1063/1.1329672
E. L. Kolsbjerg, M. N. Groves, and B. Hammer, J. Chem. Phys, 145, 094107 (2016) doi: 10.1063/1.4961868
S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys. 150, 094109 (2019) https://dx.doi.org/10.1063/1.5064465
- images: list of Atoms objects
Images defining path from initial to final state.
- k: float or list of floats
Spring constant(s) in eV/Ang. One number or one for each spring.
- climb: bool
Use a climbing image (default is no climbing image).
- parallel: bool
Distribute images over processors.
- remove_rotation_and_translation: bool
TRUE actives NEB-TR for removing translation and rotation during NEB. By default applied non-periodic systems
- method: string of method
Choice betweeen five methods:
aseneb: standard ase NEB implementation
improvedtangent: Paper I NEB implementation
eb: Paper III full spring force implementation
spline: Paper IV spline interpolation (supports precon)
string: Paper IV string method (supports precon)
- allow_shared_calculator: bool
Allow images to share the same calculator between them. Incompatible with parallelisation over images.
- precon: string,
ase.optimize.precon.Preconinstance or list of
instances. If present, enable preconditioing as in Paper IV. This is possible using the ‘spline’ or ‘string’ methods only. Default is no preconditioning (precon=None), which is converted to a list of
Example of use, between initial and final state which have been previously saved in A.traj and B.traj:
from ase import io from ase.mep import NEB from ase.optimize import MDMin # Read initial and final states: initial = io.read('A.traj') final = io.read('B.traj') # Make a band consisting of 5 images: images = [initial] images += [initial.copy() for i in range(3)] images += [final] neb = NEB(images) # Interpolate linearly the potisions of the three middle images: neb.interpolate() # Set calculators: for image in images[1:4]: image.calc = MyCalculator(...) # Optimize: optimizer = MDMin(neb, trajectory='A2B.traj') optimizer.run(fmax=0.04)
Be sure to use the copy method (or similar) to create new instances of atoms within the list of images fed to the NEB. Do not use something like [initial for i in range(3)], as it will only create references to the original atoms object.
Notice the use of the
interpolate() method to obtain an
initial guess for the path from A to B.
Interpolate path linearly from initial to final state.
- ase.mep.neb.interpolate(images, mic=False, interpolate_cell=False, use_scaled_coord=False, apply_constraint=None)¶
Given a list of images, linearly interpolate the positions of the interior images.
- mic: bool
Map movement into the unit cell by using the minimum image convention.
- interpolate_cell: bool
Interpolate the three cell vectors linearly just like the atomic positions. Not implemented for NEB calculations!
- use_scaled_coord: bool
Use scaled/internal/fractional coordinates instead of real ones for the interpolation. Not implemented for NEB calculations!
- apply_constraint: bool
Controls if the constraints attached to the images are ignored or applied when setting the interpolated positions. Default value is None, in this case the resulting constrained positions (apply_constraint=True) are compared with unconstrained positions (apply_constraint=False), if the positions are not the same the user is required to specify the desired behaviour by setting up apply_constraint keyword argument to False or True.
Interpolate path linearly from initial to final state. This standalone function can be used independently of the NEB class, but is functionally identical.
Create an improved path from initial to final state using the IDPP approach . This will start from an initial guess of a linear interpolation.
- ase.mep.neb.idpp_interpolate(images, traj='idpp.traj', log='idpp.log', fmax=0.1, optimizer=<class 'ase.optimize.mdmin.MDMin'>, mic=False, steps=100)¶
Interpolate using the IDPP method. ‘images’ can either be a plain list of images or an NEB object (containing a list of images).
Generate an IDPP pathway from a set of images. This differs from above in that more IDPP-specific parameters can be specified, and an initial guess for the IDPP other than linear interpolation can be provided.
Only the internal images (not the endpoints) need have calculators attached.
Information about energy minimization (optimization). Note that you cannot use the default optimizer, BFGSLineSearch, with NEBs. (This is the optimizer imported when you import QuasiNewton.) If you would like a quasi-newton optimizer, use BFGS instead.
How to use calculators.
If there are \(M\) images and each image has \(N\) atoms, then the NEB
object behaves like one big Atoms object with \(MN\) atoms, so its
get_positions() method will return a \(MN \times 3\)
from ase.optimize import BFGS opt = BFGS(neb, trajectory='A2B.traj')
will write all images to one file. The Trajectory object knows about NEB calculations, so it will write \(M\) images with \(N\) atoms at every iteration and not one big configuration containing \(MN\) atoms.
The result of the latest iteration can now be analysed with this command: ase gui A2B.traj@-5:.
For the example above, you can write the images to individual trajectory files like this:
for i in range(1, 4): opt.attach(io.Trajectory('A2B-%d.traj' % i, 'w', images[i]))
The result of the latest iteration can be analysed like this:
$ ase gui A.traj A2B-?.traj B.traj -n -1
Restart the calculation like this:
images = io.read('A2B.traj@-5:')
The “climbing image” variation involves designating a specific image to behave differently to the rest of the chain: it feels no spring forces, and the component of the potential force parallel to the chain is reversed, such that it moves towards the saddle point. This depends on the adjacent images providing a reasonably good approximation of the correct tangent at the location of the climbing image; thus in general the climbing image is not turned on until some iterations have been run without it (generally 20% to 50% of the total number of iterations).
To use the climbing image NEB method, instantiate the NEB object like this:
neb = NEB(images, climb=True)
Quasi-Newton methods, such as BFGS, are not well suited for climbing image NEB calculations. FIRE have been known to give good results, although convergence is slow.
Scaled and dynamic optimizations¶
- class ase.mep.dyneb.DyNEB(images, k=0.1, fmax=0.05, climb=False, parallel=False, remove_rotation_and_translation=False, world=None, dynamic_relaxation=True, scale_fmax=0.0, method='aseneb', allow_shared_calculator=False, precon=None)¶
Subclass of NEB that allows for scaled and dynamic optimizations of images. This method, which only works in series, does not perform force calls on images that are below the convergence criterion. The convergence criteria can be scaled with a displacement metric to focus the optimization on the saddle point region.
‘Scaled and Dynamic Optimizations of Nudged Elastic Bands’, P. Lindgren, G. Kastlunger and A. A. Peterson, J. Chem. Theory Comput. 15, 11, 5787-5793 (2019).
- dynamic_relaxation: bool
True skips images with forces below the convergence criterion. This is updated after each force call; if a previously converged image goes out of tolerance (due to spring adjustments between the image and its neighbors), it will be optimized again. False reverts to the default NEB implementation.
- fmax: float
Must be identical to the fmax of the optimizer.
- scale_fmax: float
Scale convergence criteria along band based on the distance between an image and the image with the highest potential energy. This keyword determines how rapidly the convergence criteria are scaled.
The convergence of images is often non-uniform, and a large fraction of computational resources can be spent calculating images that are below the convergence criterion. This can be avoided with a dynamic optimization method that monitors the convergence of each image. Dynamic optimization is implemented as a subclass of the NEB method:
from ase.mep import DyNEB neb = DyNEB(images, fmax=0.05, dynamic_relaxation=True)
fmax must be identical to the
fmax of the optimizer.
Dynamic optimization only works efficiently in series, and will not result
in reduced computational time when resources are parallelized over images.
dynamic_relaxation=False reverts to the default NEB implementation.
The saddle point is the important result of an NEB calculation, and the other interior images are typically not used in subsequent analyses. The convergence criteria can be scaled to focus on the saddle point and increase the tolerance in other regions of the PES. Convergence scaling is implemented as:
neb = DyNEB(images, fmax=0.05, dynamic_relaxation=True, scale_fmax=1.)
where the convergence criterion of each image is scaled based on the position
of the image relative to the highest point in the band. The rate (slope) of
convergence scaling is controlled by the keyword
A low scaling factor (
scale_fmax=1-3) is often enough to significantly
reduce the number of force calls required for convergence.
Parallelization over images¶
Some calculators can parallelize over the images of a NEB calculation. The script will have to be run with an MPI-enabled Python interpreter like GPAW’s gpaw-python. All images exist on all processors, but only some of them have a calculator attached:
from ase.parallel import world from ase.calculators.emt import EMT # Number of internal images: n = len(images) - 2 j = world.rank * n // world.size for i, image in enumerate(images[1:-1]): if i == j: image.calc = EMT()
Analysis of output¶
A class exists to help in automating the analysis of NEB jobs. See the Diffusion Tutorial for some examples of its use.
- class ase.mep.neb.NEBTools(images)¶
Class to make many of the common tools for NEB analysis available to the user. Useful for scripting the output of many jobs. Initialize with list of images which make up one or more band of the NEB relaxation.
- get_barrier(fit=True, raw=False)¶
Returns the barrier estimate from the NEB, along with the Delta E of the elementary reaction. If fit=True, the barrier is estimated based on the interpolated fit to the images; if fit=False, the barrier is taken as the maximum-energy image without interpolation. Set raw=True to get the raw energy of the transition state instead of the forward barrier.
Deprecated since version 3.23.0: Please use
Returns fmax, as used by optimizers with NEB.
Plots the NEB band on matplotlib axes object ‘ax’. If ax=None returns a new figure object.
- plot_bands(constant_x=False, constant_y=False, nimages=None, label='nebplots')¶
Given a trajectory containing many steps of a NEB, makes plots of each band in the series in a single PDF.
- constant_x: bool
Use the same x limits on all plots.
- constant_y: bool
Use the same y limits on all plots.
- nimages: int
Number of images per band. Guessed if not supplied.
- label: str
Name for the output file. .pdf will be appended.
You can also make NEB plots that show the relaxation of your trajectory directly from the command line; this will output the plots into a single PDF:
$ ase nebplot neb.traj
You can find more help with:
$ ase nebplot -h
- class ase.mep.autoneb.AutoNEB(attach_calculators, prefix, n_simul, n_max, iter_folder='AutoNEB_iter', fmax=0.025, maxsteps=10000, k=0.1, climb=True, method='eb', optimizer=<class 'ase.optimize.fire.FIRE'>, remove_rotation_and_translation=False, space_energy_ratio=0.5, world=None, parallel=True, smooth_curve=False, interpolate_method='idpp')¶
The AutoNEB algorithm streamlines the execution of NEB and CI-NEB calculations following the algorithm described in:
E. L. Kolsbjerg, M. N. Groves, and B. Hammer, J. Chem. Phys, 145, 094107, 2016. (doi: 10.1063/1.4961868)
The user supplies at minimum the two end-points and possibly also some intermediate images.
- The stages are:
- Define a set of images and name them sequentially.
Must at least have a relaxed starting and ending image User can supply intermediate guesses which do not need to have previously determined energies (probably from another NEB calculation with a lower level of theory)
AutoNEB will first evaluate the user provided intermediate images
AutoNEB will then add additional images dynamically until n_max is reached
A climbing image will attempt to locate the saddle point
All the images between the highest point and the starting point are further relaxed to smooth the path
All the images between the highest point and the ending point are further relaxed to smooth the path
Step 4 and 5-6 are optional steps!
Function which adds valid calculators to the list of images supplied.
- prefix: string or path
All files that the AutoNEB method reads and writes are prefixed with prefix
- n_simul: int
The number of relaxations run in parallel.
- n_max: int
The number of images along the NEB path when done. This number includes the two end-points. Important: due to the dynamic adding of images around the peak n_max must be updated if the NEB is restarted.
- climb: boolean
Should a CI-NEB calculation be done at the top-point
- fmax: float or list of floats
The maximum force along the NEB path
- maxsteps: int
The maximum number of steps in each NEB relaxation. If a list is given the first number of steps is used in the build-up and final scan phase; the second number of steps is used in the CI step after all images have been inserted.
- k: float
The spring constant along the NEB path
- method: str (see neb.py)
Choice betweeen three method: ‘aseneb’, standard ase NEB implementation ‘improvedtangent’, published NEB implementation ‘eb’, full spring force implementation (default)
- optimizer: object
Optimizer object, defaults to FIRE
- space_energy_ratio: float
The preference for new images to be added in a big energy gab with a preference around the peak or in the biggest geometric gab. A space_energy_ratio set to 1 will only considder geometric gabs while one set to 0 will result in only images for energy resolution.
The AutoNEB method uses a fixed file-naming convention. The initial images should have the naming prefix000.traj, prefix001.traj, … up until the final image in prefix00N.traj Images are dynamically added in between the first and last image until n_max images have been reached. When doing the i’th NEB optimization a set of files prefixXXXiter00i.traj exists with XXX ranging from 000 to the N images currently in the NEB.
- The most recent NEB path can always be monitored by:
$ ase-gui -n -1 neb???.traj