Some NEB calculations are getting stuck

GULP
1

Hello,

I often get the output for the NEB optimisation run similar to this:

Start of LM-BFGS optimisation:

Iteration = 0 Force norm = 0.511337 RMSD = 0.00000000
Iteration = 1 Force norm = 0.511337 RMSD = 0.00000000
Iteration = 2 Force norm = 0.511337 RMSD = 0.00000000
Iteration = 3 Force norm = 0.511337 RMSD = 0.00000000
Iteration = 4 Force norm = 0.511337 RMSD = 0.00000000

**** Optimisation stopped as force is not dropping ****

I have tried changing the nebtangent recipe (1, 2 or 3), adding nodneb keyword and increasing the nebspring constant. These changes slightly affect the value of the “Force norm” but don’t make the NEB replicas move.

A typical .gin file would look like this:

cineb conv noautobond fix molmec nebspring 0.001 nebtangent 3

and would contain the atomic coordinates of a Metal-Organic Framework and a small guest molecule that I am trying to move through the pore channel.

I suspect that this problem might be cased by some of the NEB replicas having a relatively high energy and the subsequent gradient calculations returning large forces, though taking any single replica and optimising it on its own is never a problem when using

opti conv noautobond fix molmec

Some NEB calculations work without any problems, some get stuck after 4 iterations and some run full 1,000 iterations with RMSD still being 0.00000. This seems to be a numerical problem and not chemistry specific.

Could some adjustments be made perhaps on the LM-BFGS side to make it more stable and more forgiving to the initial conditions?

I would really appreciate any thoughts on this matter.

Many thanks,
Dmytro

2

Hi Dmytro,
Sorry that you’re having issues with NEB and I agree there can be fun getting things to converge for some cases. Fixing these could take a while, but I have a couple of suggestions:

  1. There is an alternative in the code to NEB call synchronous transit, which is a method that pre-dates NEB but works quite well. I’ve realised that this keyword isn’t documented, but if you add “sync” as a keyword to an NEB input then it should work (and I’ll update the documentation!). What the method does is use the initial and final states (like NEB with 2 beads) and moves one in the direction of the other until the energy of that state is higher than for the other configuration. Then the process switches and the other image is moved until the energy is again higher, etc. In this way the two images iterative toward convergence at the transition state.
  2. GULP is written as a lattice dynamics code, rather than MD, and so it doesn’t have to assume that you only have 1st derivatives available (2nd derivatives are available for most forcefields). Therefore you can use more sophisticated optimisers for transition states than NEB. The “transition” state keyword allows you to use the second derivative matrix with the RFO optimiser to search for a place where the Hessian has exactly 1 negative mode (or any other number if you like). So if you give it a reasonable guess for the transition state, or displace the system from the nearest minima toward the transition state, then you should be able to locate a saddle point this way as well.
    Because of the above, the best strategy in GULP would be to use something like “sync” (or NEB) to find a rough guess at where the transition state is and then you can refine to convergence using RFO.

Hope that helps.
Julian

PS I’ve just updated the source code version of 6.0 to include the “sync” keyword and added example79 which is a simple illustration of how the feature works.