OptimizeNudgedElasticBand

OptimizeNudgedElasticBand(neb, max_forces=None, max_stress=None, max_steps=None, max_step_length=None, constraints=None, trajectory_filename=None, trajectory_object_id=None, log_interval=None, spring_constant=None, spring_skew=None, climbing_image=None, climbing_image_tolerance=None, preoptimization=None, optimizer_method=None, log_filename_prefix=<class 'NL.ComputerScienceUtilities.NLFlag._NLFlag.Automatic'>, pre_step_hook=None, post_step_hook=None, optimize_cell=None, disable_stress=None, target_stress=None, trajectory_interval=None, restart_strategy=None, write_raw_stress=None, write_raw_forces=None, surrogate_model=None)

Function for performing a NEB optimization.

Parameters:

• neb (NudgedElasticBand) – The nudged elastic band configuration to optimize.

• max_forces (PhysicalQuantity of type eV/Ang) – The maximum forces when the optimization should stop.
  Default: 0.05*eV/Ang.

• max_stress (PhysicalQuantity of type pressure) – The convergence criterion for the maximum difference between the internal stress and the target stress.
  Default: 0.1*GPa.

• max_steps (int) – The maximum number of steps to take before the NEB relaxation stops. When restarting, the max step are added to the current step.
  Default: 200.

• max_step_length (PhysicalQuantity of type length) – The maximum step length the optimizer may take.
  Default: 0.2*Ang.

• constraints (list of ints) – List of atom indices that are kept fixed during optimization.
  Default: [].

• trajectory_filename (str | None) – The filename used to store the trajectory. If the value is None then no trajectory file will be written.
  Default: None.

• trajectory_object_id (str | None) – The object id to use when saving the trajectory to HDF5.
  Default: None.

• log_interval –
  Deprecated: from v2019.12, use trajectory_interval instead.

• spring_constant (PhysicalQuantity of type force/length) – The spring constant used for the NEB relaxation.
  Default: 5.0*(eV/Ang**2).

• spring_skew (float) – Skew the spring constant linearly depending on the image energies. For a value larger than 1.0, higher energy images are connected by a stronger spring. As a result images are spaced closer around the transition state, which can improve convergence of the NEB optimization. This parameter determines the maximum factor by which the spring constant is increased or decreased. The current transition state estimate and the lower of the endpoints are used as references.
  Default: 1.0.

• climbing_image (bool | PhysicalQuantity of type eV/Ang | Automatic | None) – Flag indicating if the climbing image algorithm should be used to find a transition state. If a force quantity is given, the climbing image algorithm is initially disabled and will be activated after the force threshold is reached. If Automatic is given, the threshold is set to 5 times the max_forces tolerance.
  Default: False.

• climbing_image_tolerance (PhysicalQuantity of type eV/Ang) – The tolerance with which to converge the climbing image. This can be different to the forces convergence tolerance of the overall NEB. If climbing_image is False this tolerance is ignored.
  Default: Same as the max_forces tolerance.

• preoptimization (bool) – Flag indicating if the endpoints should be optimized before the NEB optimization. Preoptimization is disabled when restarting. Note that the NEB must be re-interpolated after endpoint preoptimization. The interpolation is performed with the method specified when creating the NEB object.
  Default: False.

• optimizer_method (FIRE | LBFGS | VelocityProjectionOptimization) – The optimizer to use for optimizing the structure.
  Default: LBFGS.

• log_filename_prefix (Automatic | str | None) – The logging output from each image will be written to filenames starting with this value. If it is set to Automatic then the prefix will be the name of the calling python script. If it is set to None, then all output will be written to stdout.
  Default: Automatic.

• pre_step_hook (A callable function or method) – An optional user-defined function which will be called just before each optimization step. The signature of the function requires the arguments (step, neb_configuration). The return status is ignored. Unhandled exceptions will terminate the optimization.

• post_step_hook (A callable function or method) – An optional user-defined function which will be called just after each optimization step. The signature of the function requires the arguments (step, neb_configuration). The return status is ignored. Unhandled exceptions will terminate the optimization.

• optimize_cell (bool) – If it is None, it will be set automatically according to the NEB configuration. If the NEB configuration is composed of configurations with the type MoleculeConfiguration | SurfaceConfiguration | DeviceConfiguration, it is set False. If the NEB configuration is composed of configurations of the type BulkConfiguration and the lattice vectors of the configurations are different, it is set to True; otherwise it is False.

• disable_stress –
  Deprecated: from v2022.03, use optimize_cell parameter instead.

• target_stress (PhysicalQuantity of type pressure) – The target internal stress (tensor) of the system. Can be given as a single value in case of isotropic pressure, or as an internal stress vector in Voigt notation or as a 3x3-matrix.
  Default: 0*GPa.

• trajectory_interval (int | PhysicalQuantity of type time) – The resolution used in saving steps to a trajectory file. This can either be given as an integer (a value of 1 results in all steps being saved; a value of 2 results in every second step being saved; etc.) or as a time interval.
  Default: 1.

• restart_strategy (NoRestart | RestartFromTrajectory) – The restart mechanism based on trajectory data saved in a given file.
  Default: RestartFromTrajectory().

• write_raw_stress –
  Deprecated: from v2025.06.

• write_raw_forces –
  Deprecated: from v2025.06.

• surrogate_model (GaussianProcessNEBParameters | None) – Parameters for the surrogate model used to estimate the potential energy surface. The surrogate model is not compatible with pre_step_hook, and optimize_cell or target_stress arguments. Furthermore, it cannot be restarted.

Usage Examples

Find the reaction path for conversion of ethane to ethene, i.e.

\[\mathrm{C}_2\mathrm{H}_6 \rightarrow \mathrm{C}_2\mathrm{H}_4 +\mathrm{H}_2\]

# Find the reaction pathway
optimized_neb = OptimizeNudgedElasticBand(
        neb_configuration,
        max_forces=0.05*eV/Ang)

neb_c6h6.py

The two following examples demonstrate how to calculate the reaction path for rock-salt to wurtzite in CdSe. The first example is with a very small unit cell that results in a fully-concerted reaction mechanism, while the second is in a larger unit cell that first forms a line defect.

# Perform a VC-NEB optimization
neb_configuration = OptimizeNudgedElasticBand(
    neb_configuration,
    max_forces=0.01*eV/Ang,
    max_stress=0.01*GPa,
    max_steps=500,
    climbing_image=True,
    optimize_cell=True,
)

cdse_neb_concerted.py cdse_neb_line_defect.py

Variable Cell NEB

It is possible to study phase transitions using OptimizeNudgedElasticBand. When a NudgedElasticBand contains images that are bulk configurations with different lattice vectors, the generalized solid-state nudged elastic band method (GSS-NEB) [1] algorithm is used. The algorithm optimizes both the atomic and lattice degrees of freedom along the reaction pathway and can be used to determine the reaction mechanism and energy barrier of a phase change. The target_stress parameter can be set to apply pressure or arbitrary stress.

There are two important points to be aware of when using this technique:

1. In order to get an accurate estimate of the energy barrier one must use very large unit cells. This is because if a small unit cell is used, then all atoms will move during the phase transition. This is almost certainly not the favorable pathway in real systems. A large enough unit cell should permit the formation of a local nucleation process that occurs in either the initial or final phase of the material. After the local nucleation event the new crystal structure should then propagate across the unit cell. For more details, see for instance Sec. III.C. and Fig. 9 in Ref. [1].

2. Great care must be taken when constructing the initial and final images. Just as in the regular NEB algorithm, the atom indices in the initial and final structure must match. Atom #1 in the initial structure will be moved to atom #1 in the final structure. This means that you cannot simply take two different crystal structures with the same number of atoms and expect the resulting pathway to be correct.

Optimization with surrogate potential energy surface

Optimization of the NEB can be accelerated by using a surrogate potential energy surface (PES). The surrogate model is trained on-the-fly during the optimization based on the known true energies and forces determined by the reference calculator. The model is then used to predict the the next step in the optimization, thus significantly reducing the number of calls to the reference calculator. The surrogate model is updated periodically to ensure that the predictions are accurate. The present implementation is based on Gaussian process regression, a machine learning model that is well-suited for this purpose [2][3]. The surrogate optimization loop is carried out with the All Images Evaluated algorithm, which means that all images are optimized in parallel.

The surrogate model is enabled by passing a GaussianProcessNEBParameters object to the surrogate_model argument of the OptimizeNudgedElasticBand function.

Notes on Parallelization

For non-force field calculators, the number of processes used to calculate the energy and forces on each image is determined by the ParallelParameters object that has been set on the attached calculator. By default, the update is parallelized over images and then over the per-image energy/force calculation.

For force field calculators, which do not support ParallelParameters, the number of processes per image (PPI) is automatically determined. In the case where the total number of available MPI processes is equal to or smaller than the number of interior images the PPI is set to 1. This means that the number of images optimized in parallel corresponds to the number of available processes. If there are more available processes than interior images, the processes are distributed across those images. It is generally recommended not to use more MPI processes than interior images but rely on threading instead.

Care should be taken to ensure that the total number of MPI processes is divisible by the number of “interior” (moving) NEB images to achieve good load balancing. For example, if there are 7 images total, then 5 of them are interior images that move during the simulation and the total number of MPI processes should be a multiple of 5. The calculation will still run if the number of MPI processes is not a multiple of the number of interior images, however, a warning will be printed and the calculation of some images will be slower than others reducing the parallel efficiency of the overall calculation.

The output from the calculation will be logged to a separate file for each image whose name is the specified by the log_filename_prefix parameter with the image number appended to it. If the log files already exist, they will be appended to.

If preoptimization is selected and there are enough MPI processes for at least two images to be calculated in parallel, then the preoptimization will be performed in parallel and its output will be logged to files starting with the log_filename_prefix followed by “_preoptimization”.

NEB optimizations can be restarted from a previous run, e.g. if the simulation was interrupted or the optimization did not converge in the maximum number of steps. The detailed restart behavior is specified and documented in RestartFromTrajectory.

• See also NudgedElasticBand.

External Electric Field

The inclusion of a homogeneous external electric field in bulk periodic DFT calculations with periodic boundary conditions is supported via the object ElectricFieldConstraint. For other geometries (i.e. slab or molecules), it is also possible to add an electric field by means of metallic regions in the Poisson solver (see Poisson solvers).

Troubleshooting non-convergence

Unfortunately, NEB optimizations can run into convergence issues which cannot be completely avoided even when following best practice guidelines [4]. However, there are several strategies one can try to improve convergence:

• The convergence is very sensitive to the initial guess of the path. The simplest interpolation scheme is a LinearInterpolation between the provided initial and final states. While it is guaranteed to generate a path, it can lead to unphysical configurations when multiple atoms rotate around a common center. In such cases, it is recommended to choose a more sophisticated interpolation scheme, such as the ImageDependentPairPotential (IDPP) method and its sequential variant

• Related to the previous point, the endpoints should not only be well converged, but also chosen in such a way that they minimize the path length by avoiding unnecessary rotations and translations. A failed NEB optimization can often be salvaged by taking the two images that are closest to the transition state and optimizing them to the closest local minimum. Afterwards, the NEB can be restarted with those minima as the new endpoints.

• Increasing the number of images can also aid convergence. The more images are used the more accurate the forces will be that ensure the images stay evenly spaced and on the minimum energy path. However, one must be aware that the computational cost scales linearly with the number of images (see Notes on Parallelization).

• In addition to increasing the number of images one can also enable the spring skewing option. This will push the images towards higher energy regions, which will improve convergence as these regions are usually associated with more complex atomic rearrangements. In practice, one can achieve convergence faster and with fewer images by using spring skewing [5].

• In the NEB method, the spring constant is relatively arbitrary since it only serves to keep the images equidistant. However, the spring constant can have some impact on the convergence of the NEB calculation. Both too high and too low spring constants can lead to slower convergence. The default value preforms well in most cases.

• The LBFGS optimizer is the default for NEB optimizations since it is the most robust optimizer. If the initial path is far from the minimum energy path (MEP), convergence can be slow. Another optimizer, such as FIRE or VPO might be more efficient. One might also consider starting out with FIRE and then switching to LBFGS once the path has closed in on the MEP.

[1] (1,2)

Daniel Sheppard, Penghao Xiao, William Chemelewski, Duane D. Johnson, and Graeme Henkelman. A generalized solid-state nudged elastic band method. The Journal of Chemical Physics, 136(7):074103, 2012. doi:10.1063/1.3684549.

[2]

Olli-Pekka Koistinen, Freyja B. Dagbjartsdóttir, Vilhjálmur Ásgeirsson, Aki Vehtari, and Hannes Jónsson. Nudged elastic band calculations accelerated with Gaussian process regression. The Journal of Chemical Physics, 147(15):152720, 09 2017. URL: https://doi.org/10.1063/1.4986787, arXiv:https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4986787/15536366/152720\_1\_online.pdf, doi:10.1063/1.4986787.

[3]

Olli-Pekka Koistinen, Vilhjálmur Ásgeirsson, Aki Vehtari, and Hannes Jónsson. Nudged Elastic Band Calculations Accelerated with Gaussian Process Regression Based on Inverse Interatomic Distances. Journal of Chemical Theory and Computation, 15(12):6738–6751, 2019. URL: https://doi.org/10.1021/acs.jctc.9b00692, arXiv:https://doi.org/10.1021/acs.jctc.9b00692, doi:10.1021/acs.jctc.9b00692.

[4]

Andrew W. Ruttinger, Divya Sharma, and Paulette Clancy. Protocol for directing nudged elastic band calculations to the minimum energy pathway: nurturing errant calculations back to convergence. Journal of Chemical Theory and Computation, 18(5):2993–3005, 2022. URL: https://doi.org/10.1021/acs.jctc.1c00926, doi:10.1021/acs.jctc.1c00926.

[5]

G. Henkelman, B. P. Uberuaga, and H. Jonsson. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys., 2000.