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 C2H2 endpoints, following the same procedure as in earlier
tutorials, then save the optimised structures to disk in ChemShell
# 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, trustradius="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
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.