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