Forces cannot converge in NEB calculations

[minminmin]

Dear Graeme,
I am getting into trouble in simulating the hydrogen evolution reaction (Tafel reaction) on a Pt(111) surface covered by a layer of water. Forces and energies are fully converged for both the initial state and the final state structures. However, the transition state calculation cannot converge because forces on atoms are too large.
I have tried a two-step calculation. First, I used IBRION = 3, IOPT=0, POTIM = 0.01, and EDIFFG = -1 to relax the elastic band and atoms (see the folder “first-step”). Although the CHAIN-FORCE could be converged fast, forces on atoms cannot get smaller than 1 eV/Å even after 100 ionic steps. Then I copy the obtained CONTCARs to POSCARs to start a second calculation with IBRION = 3, IOPT=1, POTIM = 0, and INVCURV = 0.001 (see the folder “second-step”). However, the convergence gets even worser (forces > 10 eV/Å)) and the geometries of images get rather unrealistic. Use of IOPT=3 in the second step cannot help to converge.
Any suggestions that can help us are appreciated.

INCAR for the first calculation:
LCLIMB = .T.
IOPT = 0
SPRING = -5
IMAGES = 7
IBRION = 3
POTIM = 0.01
EDIFFG = -1

INCAR for the second calculation:
LCLIMB = .T.
IOPT = 1
INVCURV = 0.001
ICHAIN = 0
SPRING = -5
IMAGES = 7
IBRION = 3
POTIM = 0
EDIFFG = -1

[graeme]

I did my standard calculation for problematic cases. I went to minimal precision settings and re-relaxed the endpoints. Then I did an 8 image band with IOPT=3 TIMESTEP=0.1. It seems to be converging just fine, to a barrier of about 0.77 eV with residual forces on the order of 0.03 eV/Ang. From this, you should be able to increase the precision settings and make sure that it is converged.

[anirudhnk]

Hello,

How do we re-run an NEB calculation? Like I understand for normal vasp calculations, we can copy the CONTCAR to POSCAR and include the WAVECAR. But for NEB, there would be many CONTCARS generated right? Do we simply modify the INCAR and re submit?
Thank you.

[graeme]

We have a script in the vtstscript package to help with this. For any vasp calculation, regular, dimer, or neb, you can use the vfin.pl script. This takes an input of a directory name. The script will then copy the calculation results to the specified directory as well as copying all the CONTCAR -> POSCAR files in the original directory so that it is ready for a restart. A nice thing about this system is that you can then do whatever analysis you want to in the subdirectory with your completed calculation without messing up the subsequent one.

[anirudhnk]

Got it, i shall take a look at the script. Thank you very much!

[anirudhnk]

Hello Prof. Graeme, i just wanted to ask if we would have to put the OUTCAR files in the 00 and NN directories before using the vfin.pl script?
Thank you.

[graeme]

Yes, that is better. It is not required, but if you want to use subsequent neb-related scripts for determining saddle energies or plotting the minimum energy path, they need to know the energy of the endpoints and they get that from the OUTCAR files.

[Bhubnesh Lama]

Hi All,
I am new to the neb calculation. I am looking at the MEP for motion of Li-ion in Li3Fe4P4O16. I optimize the IS and FS using EDIFFG = -0.01 and also I am using LDA+U. My INCAR for the neb calculation is:

system = LFP
ENCUT = 520 eV
LMAXMIX = 4 # L channel mixing for the mixing charge and wavefunction
PREC = Accurate

#Others variable
ADDGRID = .TRUE.
EDIFF = 1.0E-6
SYMPREC = 1.0E-4
LORBIT = 11
LREAL = Auto

NELM = 200

# Few more
ISMEAR = 0 # Gaussian Smearing of the Fermi function
SIGMA = 0.05 # Width of the the gaussian smearing ( -2 read the occupancies formt the incar or potcar file and keep it constant)

LDAU = .TRUE
LDAUTYPE = 2
LDAUL = -1 2 -1 -1
LDAUU = 0.0 5.0 0.0 0.0
LDAUJ = 0.0 0.0 0.0 0.0
LDAUPRINT = 2

# Parallization
ALGO = Normal
WEIMIN = 0 # Maximum Wight for the band to be considered empty.
ISPIN = 2
MAGMOM = 3*0.0 4*5.0 4*0.0 16*0.0

#LDIPOL = TRUE
#IDIPOL = 3
NPAR = 4
AMIN = 0.01
EDIFFG = -0.01
NSW = 200
ISIF = 2
IBRION = 3
#LNEBCELL = .TRUE. # Flag to turn on SS-NEB Used with ISIF=3 and IOPT=3
POTIM = 0.0
IOPT = 3 # QM (Quick-Min) force-based optimizers are used and good for a high force far from minimum

ICHAIN =0
IMAGES = 5
SPRING = -5
LCLIMB =.TRUE
TIMESTEP = 0.001

my forces per atom in the OUTCAR in each atom in the images are >0.01 eV and I am having difficulty in decreasing the force per atom to EDIFFG level. And my activation energy is larger (~10 eV) than in the literature (~0.4 eV).
How can I resolve these issues?

[Bhubnesh Lama]

Bhubnesh Lama wrote:
> Hi All,
> I am new to the neb calculation. I am looking at the MEP for motion of Li-ion in
> Li3Fe4P4O16. I optimize the IS and FS using EDIFFG = -0.01 and also I am using
> LDA+U. My INCAR for the neb calculation is:
>
> system = LFP
> ENCUT = 520 eV
> LMAXMIX = 4 # L channel mixing for the mixing charge and wavefunction
> PREC = Accurate
>
> #Others variable
> ADDGRID = .TRUE.
> EDIFF = 1.0E-6
> SYMPREC = 1.0E-4
> LORBIT = 11
> LREAL = Auto
>
> NELM = 200
>
> # Few more
> ISMEAR = 0 # Gaussian Smearing of the Fermi function
> SIGMA = 0.05 # Width of the the gaussian smearing ( -2 read the occupancies formt
> the incar or potcar file and keep it constant)
>
> LDAU = .TRUE
> LDAUTYPE = 2
> LDAUL = -1 2 -1 -1
> LDAUU = 0.0 5.0 0.0 0.0
> LDAUJ = 0.0 0.0 0.0 0.0
> LDAUPRINT = 2
>
> # Parallization
> ALGO = Normal
> WEIMIN = 0 # Maximum Wight for the band to be considered empty.
> ISPIN = 2
> MAGMOM = 3*0.0 4*5.0 4*0.0 16*0.0
>
> #LDIPOL = TRUE
> #IDIPOL = 3
> NPAR = 4
> AMIN = 0.01
> EDIFFG = -0.01
> NSW = 200
> ISIF = 2
> IBRION = 3
> #LNEBCELL = .TRUE. # Flag to turn on SS-NEB Used with ISIF=3 and IOPT=3
> POTIM = 0.0
> IOPT = 3 # QM (Quick-Min) force-based optimizers are used and good for a high force
> far from minimum
>
> ICHAIN =0
> IMAGES = 5
> SPRING = -5
> LCLIMB =.TRUE
> TIMESTEP = 0.001
>
> my forces per atom in the OUTCAR in each atom in the images are >0.01 eV and I am
> having difficulty in decreasing the force per atom to EDIFFG level. And my activation
> energy is larger (~10 eV) than in the literature (~0.4 eV).
> How can I resolve these issues?

[graeme]

This looks like a case of misordering atoms between the initial and final states. You can find many similar examples reported in this forum.

[Bhubnesh Lama]

Thank you, Prof. Graeme
I ordered the atoms between initial and final states and did my calculation and it gave me reasonable force and energy.