Divergence of NEB calculation

[nokimann]

I tested twice a straightforward example of the chemical reaction

2 C2H4 ---> C4H8

with different INCARs

Here is my first divergent INCAR

SYSTEM = NEB
NCORE = 4
ISTART = 0
ICHARG = 2
ENCUT = 550.0
EDIFF = 0.1E-02
NSW = 250
IBRION = 1
ISIF = 2
ISYM = 1
POTIM = 0.5000
ISMEAR = -5
IALGO = 38
LWAVE = .FALSE.
LCHARG = .FALSE.
LVTOT = .FALSE.
SPRING = -5
IMAGES = 4

And my second (qusai?-)convergent INCAR, because of NSW(?), it just calculated one loop
So the result of it is reasonable quite well, compared to a result of activation graph which I saw in another tutorial

SYSTEM = NEB
NCORE = 4
ISTART = 0
ICHARG = 2
ENCUT = 550.0
EDIFF = 0.1E-02
NSW = 0 <-------------------- only different this
IBRION = 1
ISIF = 2
ISYM = 1
POTIM = 0.5000
ISMEAR = -5
IALGO = 38
LWAVE = .FALSE.
LCHARG = .FALSE.
LVTOT = .FALSE.
SPRING = -5
IMAGES = 4

When I run with 1st INCAR, it showed very divergent CONTCAR images at the end.
Even though I set EDIFF as 0.1E-02, it couldn't converge. (Maybe this is the INCAR problem, not the VTST and NEB?)
I want to get the more precise value less than 0.1E-05 or 0.1E-06 of EDIFF.
What should I change inside INCAR or what do I have to do?

And I'm searching for the tutorials on NEB calculation, but most browsers might be out of date.
If you know any tutorials or advice, could you let me know them?

[nokimann]

I read a few posts. They mostly said "adjust POTIM and IBRION."
But I want to know how I can decide the value of POTIM and IBRION for all cases?
Does personal experience just govern it?

[graeme]

There is some art to optimization, but there are also some sensible guidelines. When forces are high, use a conservative optimizer such as IBRION=3 and a small time step < 0.1. Once the forces are below 0.2 eV/Ang, you can switch to a second order optimizer (e.g. IBRION=1).

[nokimann]

You said 'forces,' but what did you mean by 'forces'?
+ How can I check that and before or after the NEB calculation?
And 'small time step' means the TIME-tag?

[graeme]

You can see the current NEB forces using our script nebef.pl.

The time step is controlled by POTIM is you use the built-in vasp-optimizers and by the TIMESTEP tag if you use our vtst optimizers.

[nokimann]

I'm still suffering from the divergence.
I guessed the reason is that the two C2H4 want to be separated from one another until all C2H4 has the same distance, considering the adjacent cells.

And, then I expand the cell size to lower the interaction between nearest cells, but it takes too long time to finish up.
In other words, it's difficult to check the results.

Do you have any advice?...

Note 1: When I tried just 'relaxation' steps, it was convergent indeed. Troubles only happened on IMAGE folders.
Note 2: The force value is small (0.01~0.02)

[graeme]

If you can post a .tar.gz file of the calculation, I can provide more specific comments.
You can remove any large CHG* and WAV* files.

[nokimann]

Preparing the example for you.
By the way, I wonder the difference between inserting OUTCARs into 00 and 04 (if IMAGES = 4) folders and not.

[graeme]

If you include OUTCAR files in the endpoint directories, our script will read them to get the endpoint energies, for example when minimum energy paths are plotted.

[nokimann]

I've prepared four calculations.
0_NEB is the 1st one, and 3_NEB is the last one.

Below is the description of what I did.

0_NEB: the first initial and final trial POSCAR, POTIM=0.01
not convergent for the calculation
images seem to be divergent. (See each CONTCAR in the image folder)

1_NEB: the initial trial POSCAR shifted from 0_NEB to have nearly equilibrium intervals between C2H4s
not convergent for the calculation
images seem to be divergent

2_NEB: POTIM increased as 0.1
convergent for the calculation
but the result images (CONTCARs) seem to be divergent

3_NEB: POTIM increased more as 0.25
not convergent both the calculation and the images

[graeme]

A few points.

- make sure that you relax with a force criterion (ediffg < 0)
-- the relaxed endpoints are quite different in the attached path because of this

- when you have problems converging, use IBRION =3 and a small POTIM - all of you calculations used second order optimizers with unstable parameters

- Since the reactant molecules moved apart in the initial state, I needed more images to resolve the path

- turn of symmetry; you can see that the optimized path has broken the symmetry of the molecules. You would not be able to see this with symmetry turned on.

- start with low cutoffs for finding minimum energy pathways; it is easy to reconverge a path with more accurate setting when you have it.

- There is no need for any k-point sampling with isolated molecules; use the gamma point and a gamma point vasp binary.

[nokimann]

These are your comments. And I added queries for each comment.

- Make sure that you relax with a force criterion (EDIFFG < 0)
-- The relaxed endpoints are quite different in the attached path because of this

I searched the meaning of this Tag, and it means "If EDIFFG is negative it has a different meaning: In this case the relaxation
will stop if all forces are smaller than | EDIFFG |."
And force is defined by ionic positions, wavefunctions, and partial occupancies, isn't it?
Anyway, what I'm wondering is where all forces come from and how we can check?
I found the info. like the line below. So this line indicates the forces?

POSITION TOTAL-FORCE (eV/Angst)

If EDIFFG > 0, this means it stops when the change in ionic step is smaller than this Tag.


- When you have problems converging, use IBRION =3, and a small POTIM - all of you calculations used second-order optimizers with unstable parameters

I misunderstood about this. Checking the forces should be conducted on the results of NEB, not after relaxation, shouldn't it?


- Since the reactant molecules moved apart in the initial state, I needed more images to resolve the path

Ok.


- Turn off symmetry; you can see that the optimized path has broken the symmetry of the molecules. You would not be able to see this with symmetry turned on.

I tried it and also tried it with changing IALGO from 48 to 38. It gave me much faster convergent. What do you think about this?
And it has nearly symmetric results in CONTCAR. I might be caused by the smaller SC-loops compared to IALGO=48?


- Start with low cutoffs for finding minimum energy pathways; it is easy to reconverge a path with more accurate setting when you have it.

I've seen that ENCUT is usually x1.3~1.5 to the highest ENMAX in POTCARs. And it mostly guarantees the convergence in relaxation steps.
But you said 'low cutoffs.' Can you explain about this in NEB if you can?

- There is no need for any k-point sampling with isolated molecules; use the gamma point and a gamma point vasp binary.

I selected k-point as 'G' with 2 2 1. Did you mean by other things?
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Really happy for your advice. And it really works well. But I need time to understand your INCAR Tags.

[graeme]

Yes, EDIFFG<0 is a force criterion - you should always use this.

Using symmetry for a symmetric structure will be fast, but as you can see from my calculation, the lowest energy path breaks that symmetry. Very few NEB calculations will have symmetry - just turn it off for these calculations.

In terms of the cutoff and all other parameters, start with inexpensive parameters when searching the potential landscape and finding transition states. Once you have found the reaction mechanisms, it is easy to reconverge with more accurate parameters and check for convergence.

For an isolated molecule you only need the gamma point (1 1 1) and you can also run a gamma point version of vasp which should increase the speed by about a factor of 2.

[nokimann]

Thank you, professor.