Reaction path optimisation

The nudged elastic band (NEB) method [12,13] is implemented in DL-FIND to optimise a reaction path between known reactant and product geometries. In a NEB optimisation, an initial path between the endpoints is guessed (either as a linear transit or via extra geometries specified by the user), and multiple geometry images are defined along it. The energies and gradients of each image are calculated independently and then an overall NEB gradient is calculated, consisting of a spring force acting on a local tangent between each of the images, and the true force acting perpendicular to the local tangent. By optimising the NEB gradient the images describe the minimum energy reaction path.

In the DL-FIND implementation of NEB, it is possible as in the dimer method to weight coordinates so that they take a greater or lesser part in the calculation of the spring force. By setting the weight to zero, no spring force is calculated for a given coordinate and so it is simply minimised rather than taking part in the NEB force. In this way, an environmental region can be defined that is relaxed independently of a NEB inner region.

The microiterative NEB optimisation has been implemented as an extension to this weighting procedure. As in the dimer case, the system is split into an inner region where the standard NEB gradient is calculated and an outer region where the weights are set to zero. In a macroiterative step the inner region moves according to the NEB gradient, which is followed by a microiterative loop where the outer region is relaxed fully. The outer region gradient is also defined over all images, but as there is no spring force it is effectively a series of independent minimisations (albeit indirectly linked due to the fixed inner regions). In the DL-FIND implementation however, the set of outer region gradients are collated and treated as one optimisation problem using a single L-BFGS optimiser, as this is more convenient than optimising each independently. In the microiterative loop, the RMS and maximum step and gradient criteria are calculated as normal, while the average of the energies of all the images is used as the energy convergence criterion.

The microiterative scheme must also be able to handle frozen images and the climbing image. Images are frozen after the force on them falls below a certain threshold. In the microiterative case, the outer region coordinates of this image - which should be fully relaxed after the previous step - are also frozen at this point. The spawning of a climbing image introduces one further set of outer region coordinates, which are simply minimised along with all the other outer region images.

ESP fitting is also more complicated with NEB. In order to perform the microiterative cycle correctly, an ESP fit must be carried out for every NEB image during the macroiterative step. These are stored in separate ChemShell fragments indexed by image number, and then retrieved as appropriate during the microiterative gradient calculations.

The microiterative NEB algorithm was tested using solvated glycine, optimising the reaction path between the minimised standard and zwitterionic geometries. The same QM setup was used as for the other algorithms, and multiple image ESP fitting was used as described above. Image freezing and a climbing image were both used in the test. Although, as before, the standard and microiterative optimisations are not directly comparable because the microiterative environment is minimised, the energy and structure of the optimised climbing image indicated that the microiterative optimisation had converged to the same endpoint as the standard optimisation.

It should be noted that the exclusion of the environment from the NEB spring forces carries the risk that the path taken by the environment from reactant to product will become discontinuous during the optimisation, as the minimisation of each environmental image is independent. One possible solution to this problem is to perform a series of NEB optimisations with restraint terms on the environmental degrees of freedom which are gradually relaxed, as suggested by Xie et al. [14]. However, this would substantially increase the overall cost of the NEB run. In our test cases we have not observed any problems with discontinuities in the environment, so we have not chosen to implement this procedure, but it would be quite straightforward to do so if it became clear over a wider set of test cases that it is necessary.

Tom Keal 2013-06-04