Reaction path optimisation

The nudged elastic band (NEB) method in DL-FIND can be used to optimise an unknown reaction path starting from the reactant and product structures.

NEB can be run on unoptimised endpoint structures, in which case the endpoints need to be optimised as well as the other images along the reaction path. However, it is usually more efficient to pre-optimise the endpoint structures so they can be frozen during the NEB run.

In the first step of our example, min_c2h2.nwchem.py we pre-optimise the C₂H₂ endpoints, following the same procedure as in earlier tutorials, then save the optimised structures to disk in ChemShell “punch” format:

# Save the optimised structures for the NEB run
rs.save("reactant_opt.pun", "pun")
ps.save("product_opt.pun", "pun")

In the second step, neb_c2h2.nwchem.py and its equivalents, the saved optimised endpoints are first read back in to ChemShell:

# Read in optimised reactant and product state geometries from disk
rs = Fragment(coords="reactant_opt.pun")
ps = Fragment(coords="product_opt.pun")

The QM theory is then set up as in previous examples, followed by a request for an NEB run with frozen endpoints by specifying neb="frozen" in the list of DL-FIND options:

# Run an NEB optimisation between the optimised endpoints
optPath = Opt(theory=neb_theory, frag2=ps, neb="frozen", nimages=8,
             coordinates="cartesian", nebk=0.01,
             maxstep=0.9, trust_radius="const")
optPath.run()

Note that we specify the product endpoint in the option frag2=ps. You can also specify one or more initial structures along the reaction path if you know them, which should be entered as a list in frag2 (with the product endpoint last in the list). This will accelerate the convergence of the method.

The number of NEB images should also be specified in the input, in this case 8 (reasonable values are 6-10). Two of these will be the frozen endpoints, and one will be a climbing image which is created when the rest of the path is sufficiently optimised and climbs towards the transition state. The five remaining images define the rest of the reaction path.

The recommended optimiser settings for use with NEB are generally algorithm=lbfgs (the default) and coordinates=cartesian.

Once the NEB calculation has run, the reaction path can be found in the file nebpath.xyz created by DL-FIND. A summary of the full path can be found in the ChemShell output file:

             Energy       F_tang    F_perp     Dist     Angle 1-3 Ang 1-2 Sum
Img    1    -76.8725400   0.00000   0.00000   0.70698    0.00    0.00   29.53   frozen
Img    2    -76.8627917  -0.00018   0.00013   0.72530   29.53   46.76   47.99   frozen
Img    3    -76.8354527  -0.00026   0.00012   0.75105   18.46   28.45   31.52   frozen
Img    4    -76.7965304   0.00022   0.00008   0.72900   13.06   41.92   49.34   frozen
Img    5    -76.7813373   0.00010   0.00021   0.71923   36.28   56.55   59.43   frozen
Img    6    -76.7933911   0.00001   0.00020   0.71786   23.15    0.00   23.15   frozen
Img    7    -76.7993514   0.00000   0.00000   0.00000    0.00    0.00    0.00   frozen
Cimg   5    -76.7805387  -0.00011   0.00014   0.14259    0.00    0.00    0.00

Here the Energy column lists the energies of each image along the path (where Cimg is the climbing image, i.e. the optimised transition state, which also states the regular NEB image it is closest to, in this case image 5). F_tang and F_perp describe forces on the image and should approach zero at convergence. Dist and the Angle entries describe the path itself, with distances and angles calculated over all the coordinates of neigbouring images. The Dist column can be used to plot energies vs. NEB reaction coordinate by summing all Dist values up to the given image, i.e.:

              R.C.        Energy
Img    1      0.00000    -76.8725400
Img    2      0.70698    -76.8627917
Img    3      1.43228    -76.8354527
Img    4      2.18333    -76.7965304
Img    5      2.91233    -76.7813373
Img    6      3.63156    -76.7933911
Img    7      4.34942    -76.7993514

The final converged energy given by DL-FIND is that of the climbing image.

If desired the climbing image result from NEB can be used as an input guess for further refinement using P-RFO or the dimer method. If you want a tightly converged transition state this is a much more efficient strategy than attempting to tighten convergence for the whole reaction path using NEB.