NEB calculation on the defect surface of Ceria (111)

[ivast]

Dear All,

I am calculating the migration energy of the oxygen/vacancy on CeO2 surface, when its near neighbors are alkali metals (for example K, Li). The problem I have is that although the Li doped surface of Ceria (111) is converging nicely it has a problem with NEB
convergence. From the preliminary calculation it seems that Li case has two maxima close to start and final image in the NEB.
(see attachment). So I tried to split the NEB-path in two parts and I increased the number of images on the half-path.

The problem is that MEP path is changing considerably after some number of ionic iterations and it doesn't seem to converge.
For full MEP path I have done 295 ionic steps and forces were high as 2 eV/A (image 6). After 66 ionic iterations on the half-path
the biggest force is 0.7 eV/A (image 2). What I can conclude from the NEB movie and previous NEB calculation with Potassium instead of Lithium is that the Li is rather small and it is moving considerably together with a migrating oxygen. This possibly can increase locally the energy (create metastable states) and maybe NEB has then the problem in sampling the energy hypersurface (see attached movie).
The K case converged to 0.01 eV/A using FIRE and LBFGS optimizer in 179 ionic steps. Because ionic radius of K is rather large comparing to Ce+4 it didn't move much in the NEB calculation .
What I can change to increase the probability that in Li case the right MEP-path is found?
I have tried with reducing the MAXMOVE =0.02 in LBFGS (this helped in the case of K-doped ceria surface), increasing the number of images and using a different optimizers.

This my INCAR:

SYSTEM = CeO2-CeOVac-Li-surface
#LREAL=.FALSE.
# NEB stuff
IMAGES=6
ISTART=1
SPRING=-5
ISYM=0
ICHAIN=0
LCLIMB=.TRUE

LREAL = Auto
#LCHARG = .TRUE
#LWAVE = .TRUE
#LELF = .TRUE
PREC = Accurate
ISMEAR = 0
#NPAR = 1
# Restart from current
# WAVECAR and CHGCAR
ICHARG = 1
ISTART = 1
# LORBIT = 10
SIGMA = 0.01
ISIF = 2
# choose force-based minimizer
IBRION = 3
NSW = 90
# Keep POTIM small for
# better convergence
POTIM = 0
EDIFF = 1E-4
EDIFFG = -1E-2
NSIM = 4
NPAR = 4
LPLANE = .True.
NELM = 90
NELMIN = 2
NELMDL = 0
# IALGO = 48

IOPT=1
MAXMOVE=0.02
ILBFGSMEM=20
LGLOBAL=.TRUE.
INVCURV=0.01
LLINEOPT=.FALSE.
FDSTEP=2E-3

# Use FIRE optimizer for now
#IOPT = 7
#MAXMOVE = 0.2
#TIMESTEP = 0.1
#FTIMEMAX = 1.0
#FTIMEDEC = 0.5
#FTIMEINC = 1.1
#FALPHA = 0.1
#FALPHADEC = 0.99
#FNMIN = 5

I would be grateful for any comments or help.

Best regards,
Toni

[ivast]

Dear All,
I noticed that the spring forces for the Li case described in my previous post, are comparable or bigger
than real forces on the images. So I tried to reduce the spring forces by reducing SPRING constant from
-5 to -2 eV/A^2. The first iterations showed that spring forces are smaller than the real forces and the simulation seems to converge to particular MEP. But I would like to ask if this is really a good strategy in the case that NEB is converging slowly?

Thank you for any comments or help,
Toni