NEBs converging slowly

[spbatchu]

Hello Dr. Henkelman

I am running a few climbing NEB calculations on ZrO2 and TiO2. Although the NEB profiles look fine in terms of change in slope, I see that they are converging pretty slowly. I have tried both iopt=2 and iopt=3, but I don't see any improvement. I use a ediif of e-06 and edifg of -0.02 eV/Ang. I am unable to identify if there is any other factor that is causing the slow convergence.

Attached below are the folders. Your suggestions to improve the speed shall be really useful.

Thanks in advance

[graeme]

Goodness, looking at the date stamps on those slurm jobs, these have been epic calculations!

I see a few things. First, these calculations have both soft and stiff modes which makes optimization challenging. But also, you are getting just ~20 iterations per run which makes things slow-going. So while you might do a little better with a second-order optimizer, I don't think you will find anything along those lines that help even by a factor of 2-3. But I have a more significant suggestion for you:

I always recommend finding minima and transition states using modest computational settings. If I were to start calculations myself to find these barriers, I would restart with soft (O_s) potentials, change Zr_sv -> Zr, use a cutoff of about 240 eV, and start with gamma point. The kpoint sampling along gets you a factor of 10 speedup, ENCUT another factor of 3, and fewer electrons probably a factor of 4. And these all multiply, so that I honest expect that you do these calculations 100x faster to find very sensible minima and saddle geometries.

Another trick is to switch from a moderately converged NEB calculation to the dimer method (you can use the neb2dim.pl script to help) which allows you to find saddles precisely with a single image. That gets you a factor of 4-8 depending upon the number of images.

But this is the key point: it is very easy when you have a moderately-well converged geometry to simply up the precision and reconverge it. The other benefit of this strategy is that you get a measure of how sensitive each parameter is. If you go from a Zr to a Zr_sv potential, I would expect no difference - but you can tell by doing the test. It could take just a few iterations to reconverge and see if there is any difference (probably not!). With 1x1x1 -> 2x2x1 kpoints, I would expect some difference, but then almost nothing from 2x2x1 to 3x3x1. But again, this is what you can check with minimal computational work.

One last point: when you check for convergence, make sure to only check relative energies that you care about, such as binding energies and barriers. Do not try to converge absolute energies; that is very hard and they are meaningless. DFT calculations rely on cancelation of errors.

[spbatchu]

Thanks for the suggestions, Dr. Graeme. I will try them and let you know what I get.

I have a follow-up question, as you mention in the first paragraph that the calculations have both soft and stiff modes, how to check for the softness and stiffness of modes? I am unaware of the procedure and hence I am asking you.

Thanks in advance.

[graeme]

Formally, you would look at the eigenvalues of the Hessian matrix. But you can also just see it from the structure where you have both stiff C-H bonds in a molecule that is partly dissociating from the surface, so there are also going to be floppy modes associated with bending / torsions / rotation of the molecule. Codes that focus on molecules, such as Gaussian, typically work in internal coordinates to separate the stiff and soft modes. Vasp, which targets materials, works in Cartesian coordinates and so the soft and stiff modes are coupled. Quasi-Newton methods such as LBFGS aim to separate these modes, but they are still limited in their Cartesian representation.

[spbatchu]

Also, another question related to the usage of neb2dim.pl:

Just to clarify, the usage of this script is as follows right:

neb2dim.pl (image no.)

where image no. is the image number of the highest energy image, and the above command has to be used in the directory where my NEB is, right?

Thanks in advance

[graeme]

No, so we typically run things in the following way:

vfin.pl dir (copies the current files to directory dir and sets up a new calculation in the current directory)
cd dir
nebresults.pl (generates plots of the minimum energy path and calculates extrema along the path)
neb2dim.pl (uses the calculated estimate of the saddle location and tangent to the path to set up a dimer calculation)

[spbatchu]

I see. Thanks. I have just tried doing this by running the following and I face an issue:

vfin.pl new_folder
cd new_folder
nebresults.pl

After using the command nebresults.pl, I get the following:
No OUTCAR in 00
Unzipping the OUTCARs....

Then, I use the command neb2dim.pl, and I get the following:
NO EXTS.DAT IN THIS DIRECTORY; RUN nebresults.pl

and when I run nebresults.pl, the above repeats.

Am I missing something here?

[graeme]

Put the OUTCAR from your calculations of the minima in the appropriate folders (00 and NN); our scripts greps those for the energy of the endpoints which is needed for the interpolation.

[spbatchu]

Thanks Dr. Graeme

while using the nebresults.pl script in one of my cases, I have found the following:

Forces and Energy:
0 0.019946 -1428.686600 0.000000
1 0.149487 -1428.505800 0.180800
2 0.234107 -1428.723800 -0.037200
3 0.089871 -1428.953500 -0.266900
4 0.071976 -1429.017500 -0.330900
5 0.041199 -1429.066700 -0.380100
6 0.029083 -1429.149000 -0.462400
7 0.012219 -1429.202400 -0.515800
8 0.012344 -1429.219700 -0.533100

Extrema 1 found at image 0.047705 with energy: -0.000219
Extrema 2 found at image 1.186720 with energy: 0.199535

I am just trying to understand why two Extrema are predicted and also trying to understand the meaning and units of the numbers mentioned after the words 'image' and 'energy'

Thanks in advance

[graeme]

The algorithm to determine extrema simply goes between each pair of images and fits a cubic function to the energies and forces in that interval. In this case, there is a spurious minimum near the initial state with an energy of -0.000219. The position of 0.047705 is a scaled position along the path, which goes from 0 to 8 in this case. The approximation to the saddle is 18.7% from image 1 to image 2 and has an estimated barrier of 0.199535 with respect to the initial state (image 0).

[spbatchu]

Thanks for the explanation. That was helpful.