.. |C2H2| replace:: C\ :sub:`2`\ H\ :sub:`2`\ 

**************************
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
"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, 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 ``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.