Reaction path optimisation with DL-FIND
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.
In the example neb.chm, the pre-optimised reactant (rs.c) and product (ps.c) fragments are defined in the input file.
# Reactant state
c_create coords=rs.c {
coordinates
c 7.85748346687591e-02 -1.88972666351032e-06 -2.46312642501926e-01
c -7.01226542208765e-01 3.77945332702064e-06 1.74531186215820e+00
h -2.48123189618235e+00 1.88972666351032e-06 2.42431899825741e+00
h 1.37622556805467e+00 -3.77945332702064e-06 -1.77654337472603e+00
}
# Product state
c_create coords=ps.c {
coordinates
c -4.49883447328575e-01 -0.00000000000000e+00 5.02474540374067e-01
c 6.74900760059402e-01 9.44863331755160e-06 2.69416629662669e+00
h -2.51752220702840e+00 -3.77945332702064e-06 2.65293057110223e-01
h 5.64846858629888e-01 -3.77945332702064e-06 -1.31515905092333e+00
}
The NEB method can be selected using the keyword neb. We recommend that you pre-optimise your starting structures to speed up the convergence of the NEB run, in which case neb=frozen should be set to tell DL-FIND that the endpoints do not require further optimisation.
# Optimise the reaction path
dl-find coords=rs.c coords2= { ps.c } \
theory= mndo: { accuracy=veryhigh hamiltonian=am1 } \
list_option=full \
neb=frozen nimage=8 optimiser=lbfgs \
coordinates=cartesian \
result=ts_neb.c
The reactant should be specified using the coords option, while the product should be specified using coords2. 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 coords2 (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 optimiser=lbfgs coordinates=cartesian. For the full list of options that can be selected when using the nudged elastic band method, please consult the DL-FIND manual pages.
Once the reaction path is converged, the final climbing image structure ts_neb.c is exported as an XYZ file. This is the transition state found by the NEB method.
The path itself 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:
NEB Report
Energy F_tang F_perp Dist Angle 1-3 Ang 1-2 Sum
Img 1 0.1347777 0.00000 0.00000 0.71541 0.00 0.00 53.58 frozen
Img 2 0.1452152 0.00077 0.00020 0.63886 53.58 60.30 62.89 frozen
Img 3 0.1918941 0.00073 0.00015 0.56573 9.31 25.20 26.49 frozen
Img 4 0.2147686 0.00028 0.00013 0.53769 17.18 31.87 33.29 frozen
Img 5 0.2026017 -0.00004 0.00015 0.54140 16.11 42.38 47.86 frozen
Img 6 0.1858755 -0.00014 0.00013 0.55539 31.75 0.00 31.75 frozen
Img 7 0.1798594 0.00000 0.00000 0.00000 0.00 0.00 0.00 frozen
Cimg 4 0.2148332 0.00004 0.00005 0.03044 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 4). 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 0.1347777
Img 2 0.71541 0.1452152
Img 3 1.35427 0.1918941
Img 4 1.92000 0.2147686
Img 5 2.45769 0.2026017
Img 6 2.99909 0.1858755
Img 7 3.55448 0.1798594
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.