Skip to main content

Tips for the NEB Method

  • Updated

The NEB (Nudged Elastic Band) method, a calculation technique for reaction pathway analysis, is introduced as an Example within Matlantis. However, depending on the system, the calculation may not converge or complete using the example code as is. This article introduces tips and settings to try when NEB calculations do not converge.

 

1. Devise the Creation of the Initial Path

Most failures in NEB calculations are due to the inability to create an appropriate initial path. Here, we introduce four methods for creating initial paths. Please select the appropriate method according to your system.

  • Creating Initial Paths Using Linear Interpolation

    This is a common method that linearly interpolates from the initial state to the final state. In ASE, it can be used as follows.

    from ase.mep.neb import NEB
neb = NEB(images, k=0.1, climb=False, allow_shared_calculator=False, parallel=True)
neb.interpolate()

    However, in organic reactions involving twisting, such as the isomerization reaction from the cis to trans form of stilbene, atoms may get too close, making it impossible to create an appropriate initial path. (Reference: Below are the SMILES for stilbene. You can confirm that linear interpolation does not work well in practice.)

    Stilbene (cis form): c1ccccc1/C=C\c1ccccc1

    Stilbene (trans form): c1ccccc1/C=C/c1ccccc1

     

  • Creating Initial Paths Using IDPP (Image Dependent Pair Potential)

    This interpolation method uses the distance matrix between atoms. As a result, it can generate physically more reasonable initial paths than linear interpolation. For more details, please refer to the paper.

    In ASE, you can execute it simply by specifying "idpp" as an argument to interpolate. In the stilbene example above, using IDPP also allows creating an appropriate initial path.

    neb = NEB(images, k=0.1, climb=False, allow_shared_calculator=False, parallel=True)
neb.interpolate("idpp")

    *Note: IDPP tends to work well for isolated molecular systems but may not generate appropriate paths for crystalline systems.

     

  • Creating Initial Paths Using Restraint Scan

    Restraint Scan is a function to continuously change molecular structures and obtain different structures. It is similar to the scan implemented in Gaussian. By using this function, you can sequentially change simple structures to obtain more complex structures.

    In the stilbene example, you can create a reasonable initial path by rotating the dihedral angle formed by the carbon atoms of the double bond and the ipso carbon of the phenyl group by 180°.

    from matlantis_features.utils.calculators import pfp_estimator_fn
from matlantis_features.features.reaction.rest_scan import (
    RestScanFeature,
    RestScanFeatureResult,
    ScanDihedral,
    ScanDistance,
    calculate_dihedral,
)

# Please adjust the indices as needed
current_dihedral = calculate_dihedral(cis_atoms, [5, 6, 7, 8])
scans = [
    [ScanDihedral(indices=[5, 6, 7, 8], direction=1, destination=current_dihedral + 180.0, exceed=1.0, fmax=0.05)],
]

scan_feature = RestScanFeature(
    scans=scans,
    trajectory="output/scan.traj",
    estimator_fn=pfp_estimator_fn(model_version="v8.0.0", calc_mode="PBE"), 
    show_progress_bar=True, show_logger=True, logger_interval=10
)

scan_result = scan_feature(cis_atoms, fmax=0.05, n_run=1000)

    For detailed explanations and usage, please see here.

     

  • Manual Creation of Initial Paths

    If the above methods do not produce an appropriate initial path, manually manipulating atomic coordinates can be more reliable and faster. With ASE alone, you can directly adjust atomic coordinates as shown below.

    atoms.positions[0] += [1., 1., 1.] # Move atom at index 0 by (1,1,1)

    However, if many atoms need adjustment, repeatedly modifying coordinates by code and checking structures is inefficient. Therefore, we recommend using the GUI feature of pfcc_extras that we provide. Running the code below will launch the GUI.

    from pfcc_extras import show_gui
show_gui(images, representations=["ball+stick"], replace_structure=True)

    If images contains multiple structures, you can select the target structure using the lever at the bottom left. After selecting a structure, input the indices of the atoms you want to edit (multiple indices allowed) into "Selected atoms," and press buttons like "X+" to visually change their coordinates.

 

2. Review NEB Method Settings

  • Perform Regular NEB First, Then Switch to Climbing Image NEB

    Climbing Image NEB (CI-NEB) is a method to find transition states with higher accuracy. However, calculations may not converge if CI-NEB is used from the start. We recommend a two-step calculation procedure: first optimize a rough reaction path using regular NEB, then use the obtained path as the initial value for CI-NEB.

     

  • Divide the Reaction Path

    Multi-step reactions involving stable intermediates are not inherently the direct target of the NEB method. The NEB method is designed to estimate transition states by interpolating between predefined initial and final states. Therefore, attempting to calculate a multi-step reaction involving stable intermediates all at once may result in a failure to converge. In such cases, please try splitting the path by reaction step—such as 'initial state to intermediate state' and 'intermediate state to final state'—and calculate each separately.

     

[Conclusion] 

If calculations still do not converge after trying these initial path creation and setting adjustments, please consider using other reaction pathway analysis methods such as the Reaction String method without insisting on the NEB method.

Related to

Related articles