Convergence issues with CI-NEB/NEB calculation in VASP

[nuha_sahir]

Hello everyone!
I am carrying out NEB calculation to find the transition state (TS) for the dissociation of H2. I have optimized the initial and final states and then used vaspkit (I used VTST and they generate the same images) to generate the intermediate images. There were some issues with the initial images in terms of the spacing between the H-atoms so I had to manually adjust them. I ran them using VASP, and no matter what settings I use, the calculation doesn't converge. It runs for the entire 300 steps (I tried 1000 steps as well) that I set and when I check the energy value it looks like it is oscillating. I have tried a few different suggestions that I found online, namely:
1. increasing the maximum number of ionic steps (NSW) and electronic steps (NELM)
2. changing the IBRION, PREC, ALGO, POTIM values
3. increasing the k-point mesh
4. trying to carry out the calculation without spin polarization
5. changing the convergence criteria (EDIFF and EDIFFG)
6. not using the climb feature, so NEB calculation instead of CI-NEB
I have double checked that the initial and final states are fully converged (they are). I don't know what else to try. I know that this should converge since this has been done in the literature. I just don't understand what I am doing wrong.
I have attached the INCAR file along with this post. For some more reference, I run the calculation in the cluster and VASP version is 6.2.0.

INCAR:

Code: Select all

SYSTEM = H2 dissociation
ISTART = 0 ! ISTART = 0 begin from scratch... ISTART = 1 if a WAVECAR exists
ICHARG = 2 ! ICHARG = 2 if ISTART = 0... ICHARG = 1 read the charge density from CHGCAR... ICHARG = 11 for restarting calculation DOS
ISIF = 2 ! 2: Relaxation only of internal parameters... 3: change of internal parameter, shape and volume simultaneously... ISIF = 0 for AIMD (IBRION = 0)
PREC = Normal ! PREC = Normal for standard precision
ALGO = Normal ! Default ALGO = Normal
!ISYM = 0 ! 0 for no symmetry
LREAL = Auto ! determines whether the projection operators are evaluated in real-space or in reciprocal space (Default LREAL = .FALSE. for reciprocal space)
IBRION = 2 ! -1 for NSW = 0... 2 for conjugate gradient algorithm... determines how the ions are updated and moved (IBRION = 0 for AIMD) (IBRION = 1 for DIIS, close to the local minimum) (IBRION = 5 for second derivatives in DOS)
POTIM = 0.15 ! 0.02 A stepwidth (default 0.5)... time step in molecular dynamics or the step width in ionic relaxations.
ENCUT = 500 ! cutoff energy for the plane-wave-basis set in eV. 1-1.3*Largest ENMAX on the POTCAR file
ISMEAR = 0 ! Gaussian smearing; determines how the partial occupancies fnk are set for each orbital
SIGMA = 0.05 ! determines the width of the smearing in eV (default 0.2)
EDIFF = 1E-6 !* global break condition for the electronic SC-loop (in eV) (Default 1e-4)
EDIFFG = -0.01 !* (Default: EDIFFG = EDIFFx10)... defines the break condition for the ionic relaxation loop (positive for energy, negative for forces)
NSW = 300 ! maximum number of ionic steps (default NSW = 0, for sc calculation) or steps in AIMD
NELM = 300 ! sets the maximum number of electronic SC (self-consistency) steps (Default 60)
NELMIN = 5 ! specifies the minimum number of electronic self-consistency steps for each ionic step; change in AIMD to a value between 4 and 8
IVDW = 11 ! specifies a vdW (dispersion) correction (Default IVDW = 0 no dispersion)... 1/10 for DFT-D2 Grimme, 11 for DFT-D3 Grimme, 12 for DFT-D3 Becke-Johnson, 13 for DFT-D4, 2/20 for TS...
ISPIN = 2 ! spin polarized calculation (default ISPIN = 1 non-spin-polarized)
MAGMOM = 60*0.7 0.7 -0.7 ! Default, NIONS*1.0 for ISPIN=2 initial magnetic moment for each atom, if and only if ICHARG=2, or if ICHARG=1 and the CHGCAR file contains no magnetisation density
LWAVE = .FALSE. ! do not write WAVECAR
LCHARG = .FALSE. ! do not write CHGCAR
NCORE = 4 ! determines the number of compute cores that work on an individual orbital (Default NCORE = 1)
SYMPREC = 1E-4
!-------------------------------------------------------------------------------------------------
IMAGES = 4 # Number of images that will be used
SPRING = -5 # Spring force (eV/A2) between images
LCLIMB = FALSE # Turns on the climbing image algorithm, making it CI-NEB

[marie-therese.huebsch]

Hi,
I am sure we can get to the bottom of it. Is it the electronic or the ionic convergence that is an issue? Can you please share a minimal reproducible example in accordance with the forum guidelines? And if you have the stdout, OUTCAR or OSZICAR of some of your trials with different ALGO, etc. that would be great.

At first glance, the only thing that looks suspicious to me is SYMPREC = 1E-4. If the actual ground state of the calculation is symmetry broken, this implies that VASP tries to go there every electronic step and then the symmetry broken part is symmetrized away.

For the example, please isolate the VASP calculation. We cannot assist with vaspkit.

Cheers,
Marie-Therese

[nuha_sahir]

Hello! Thank you so much for your reply.

The convergence issue is with the ionic steps. I see from the OSZICAR that the electronic steps converge easily, but the ionic steps don't. In the zip file that I have attached, I have added the minimal reproducible example (I have tried my best to make it minimal reproducible based on the guidelines, but I do apologize if I misunderstood anything), along with two sets of OSZICARs that I got using different algorithms and spin polarizations.

Thank you once again. Any help would be much, much appreciated!

[marie-therese.huebsch]

Hi,

thank you for preparing the minimal example!
Your system is a monolayer with an adsorbent. As VASP uses periodic boundary conditions one has to use tricks to model a surface:
Use a vacuum padding (like you did) and dipole corrections (missing!) such that the interaction between the layers is neglectable.

Basically, a dipole forms when charge transfers from or to the surface. This almost always happens except if the slab is completely symmetric (or sometimes for metals its is negligible). This dipole induces an electric field (i.e. gradient in the potential) which causes long-range interactions. Hence, the layers interact despite the vacuum padding being sizable. The solution is to introduce a step in the potential in the vacuum. Read more here: Electrostatic corrections.

Can you please check if this solves the issue? Let me know if the system converges!

Cheers,
Marie-Therese

[nuha_sahir]

Hello!

Sorry for the late reply, I had some difficulty getting the calculations running (technical problems from the clusters side).

I added the dipole corrections and the system (and the calculations) looks much better now, however, the calculations are still not converging. The electronic steps converge easily (like before), and the ionic steps are the problem. Another thing I noticed is that my system has a layer of carbon atoms (adsorbent) and H2 dissociates and adsorbs on it. So in all images, the spins of the H atoms are opposite. However, in one particular image (03), the spins are the same, and I am wondering if this is the problem.

I have attached the ZIP file with the POSCARs and the stdout. Again, thank you so much for all the help. I really hope we can figure out the problem and get to the bottom of it.

[michael_wolloch]

Dear Nuha Sahir,

sorry for the delayed response.

I am pretty sure the spins are not the issue unless specifically image 03 is not converging. The NEB images are only coupled via the energy and not directly via the electronic states. I doubt that the different spin states would change much, but of course, this is something you can test by running static calculations and trying to converge to different magnetic states using MAGMOM. (You could also constrain the magnetic moments if necessary.)

I would suggest these possible steps to try and achieve convergence:

1. Lower the number of intermediate structures from 4 to 3 or even only 1 and see if you can converge the NEB calculation then.

2. Decrease the value of SPRING to lower the spring constant between your images. You could try SPRING = -4 or -3. For your system I do not believe that the energy barrier is very high, thus a smaller spring force could reduce oscillations.

3. Try to pre-converge the images by applying a low dimensionality parameter via setting NFREE = 2 while still using IBRION = 1. Alternatively you can switch to a steepest descent minimization without line optimization for this task by setting IBRION = 3 and SMASS = 2.

More information about these steps can be found in the nudged elastic bands section on the wiki.

As a side note, please heed the warnings printed in your output and adapt your NPAR parameter for more efficient calculations.

Cheers, Michael

[nuha_sahir]

Hello Michael! Thank you for your reply!

I increased the number of images between the IS and the FS to 8, and that solved the issue with the spins. Increasing the number of images also helped with finding the TS. Since this was a data reproduction and the energy barrier is reported in the literature, the energy barrier I get from the "TS" image corresponds to the one in the literature. However, I STILL don't have ionic convergence. The forces between the H atoms - especially around the TS, are still above the EDIFFG criterion (-0.01, -0.05, -0.1) I am setting. I checked the OUTCAR of the images and the total force I see between the H atoms for some images is close to 2 eV/A. I tried implementing your suggestions and others that I found on the other posts, but they didn't help either. Here are some things I tried:
1) Decreased the number of intermediate images to 3, 2, and 1. Decreasing the images decreases the force compared to when I had 4 images, however, they are still high and there is no ionic convergence. And the energy barriers obtained from these aren't that great when compared to the literature.
2) Decreasing the value of the spring constant to -3. This did not help at all. The forces remained more or less the same. Same thing with IBRION = 1 and NFREE = 2, and IBRION = 3 and SMASS = 2.
3) Using slightly higher ENCUT (550 instead of 500), increased k-point mesh (4x4x1 instead of 2x2x1), and I even tried adding ADDGRID. None of these helped with decreasing the forces.
4) I tried a few VTST tools like using different IOPT optimizers and they did not help either.

I've noticed that the force does initially decrease significantly, however, after a while the force value stabilizes and then no matter how many iterations take place, they remain more or less the same (if they change, the change is very small). I am at a loss for what to do. I do know that the "TS" image I have is good in terms of the energy barrier I am getting. However, due to the fact, that there is no ionic convergence, I cannot guarantee that this is the TS and I cannot proceed with my other calculations since all of them require finding TS.

I am relatively new to VASP, so I don't know if I am doing something wrong or if it is my system itself that is the problem. I will appreciate ANY suggestion that could help in decreasing the forces. I have attached the INCAR and screenshots of two of the OUTCARs for your reference. Both screenshots are of the total forces acting on the atoms for the final iteration (1000th). One of the screenshots is of the first image, so the force is low. The second image is from the 6th image, and the H atoms in this image experience the highest force.

Again any help is much, much appreciated.

Nuha

[michael_wolloch]

Dear Nuha,

thanks for the update. I am sorry to hear that the proposed modifications did not result in convergence.
I also apologize that I did not realize sooner that you are using VTST tools to make a climbing image calculation, although it is even in the topic title .
Since we do not provide the climbing image code in our release version of VASP, I cannot really comment much on that, but in your case, I would try to use one of the optimizers they provide, for example LBFGS, and see if convergence improves then.
Do not forget to set IBRION=3 and POTIM=0 in case you select one of their optimizers!

I will try to run some simulations with a bit less accuracy (Γ-only, 400eV cutoff, only 2 images) in the meantime and see if I can figure something out.

Have a good weekend, Michael

[michael_wolloch]

Dear Nuha,

a quick update:

As stated in my previous post, I ran a NEB simulation (without CI or any VTST functionality) with reduced settings over the weekend.

Indeed I also did not converge in 500 ionic steps. However, when I restarted from the previous CONTCARs from fresh orbitals (ICHARG=2), I converged to 0.01 eV/Angstrom in 10 steps.

Please try that with your system as well and report how that goes.

Cheers, Michael

[nuha_sahir]

Hello Michael!

Thank you so much for running the simulation! I ran the NEB simulation with 10 images for an increased number of ionic steps (1000 steps) and monitored the change in the forces as well. The decrease is very small but it does decrease, and with a large enough number of ionic steps, I believe it will converge. However, since I had an idea of where the TS takes place (and I was short on time), I used it to run a dimer calculation, and it converged. Also, I used the fire optimizer from VTST (IOPT = 7) and the NEB calculation converged. So I finally have results after weeks!

Again thank you so much for all your help! Have a nice week!

Nuha

[michael_wolloch]

Dear Nuha,

glad to help and super glad that you have your results now.

I do not think that it is a good allocation of resources to converge with 1000 ionic steps, however. In that case, it is better to switch algorithms or methods as detailed in our conversation above.

If you have no other questions, please indicate that this issue is resolved and I can lock the topic,
Cheers, Michael