The convergence problem of NEB

[mwchang]

Dear Prof. Graeme

Thanks for your VTST code and publications of about NEB theory.

I have successfully implemented the code into VASP and also learned a lot of knowledge from the papers.




Recently, I encountered convergence problems of NEB calculations.

I used six images to find the reaction barriers and chose IBRION=3 and POTIM= 0.01 as an optimizer in the beginning. After the NEB force is smaller than 1eV/A, I changed the optimizer to IOPT=7.

When the NEB force is smaller than 0.1eV/A again, I switched IOPT=7 to LBFGS optimizer.

I hope the NEB force can be smaller than 0.01eV/A finally, however; I saw the NEB forces oscillate at level 0.04~0.15 eV/A.


I browsed the forum and read your paper saw you have mentioned an aggressive optimizer such as LBFGS(line) and IBRION=1 is unstable, so I changed the optimizer to FIRE. However, the method do not work.


Following is my INCAR file , could you give me some suggestions?

SYSTEM = HO+HO==> H2O+O on Pt(111)

ISTART = 0 !job : 0-new 1-cont 2-samecut
NWRITE = 3
ICHARG = 2 ! 1-read CHGCAR; 2-construct initial


ISYM = 0 ! 0 for broken Symm.

NPAR=2

ENCUT = 500.0 eV
EDIFF = 1E-5
EDIFFG = -0.01

NSW = 50 ! OPT cycle


IBRION = 3 ! Disable the optimizer of vasp
POTIM = 0.00 ! Disable the optimizer of vasp

IOPT = 7 ! 0-vasp, 1-LBFGS, 2-CG, 3-Quick Min, 4-Steepest, 7-Fast Relaxation
ICHAIN = 0 ! To run NEB
IMAGES = 6 ! Number of Images in the band
SPRING = -5 ! Spring force between images
LCLIMB = .FALSE. ! Use ci-NEB algorithm
LTANGENTOLD = .FALSE.
LDNEB = .FALSE.
NEBCELL = .FALSE.


ISMEAR = 0 ! Gaussian
SIGMA = 0.1 eV ! DOS related value

ISIF = 0 ! caculate forces on atoms, relax atoms
GGA = 91 !GGA-PW91
VOSKOWN=1
ISPIN = 2 !open spin-polarized


Thank you.

[graeme]

For LBFGS, there are two considerations. First, the forces must be accurate enough to estimate curvatures. You could try ediff=1e-6 to see if that gives you smooth convergence. Second, there is an inverse curvature, INVCURV, which is 0.01 by default. Dropping this to 0.001 will make the optimizer significantly more conservative when making an initial step.

[mwchang]

Thanks for your reply.

Actually, I have changed EDIFF value to 1E-7 when LBFGS is empolyied, but it did'nt work.

I will try to modify INVCURV to 0.001. Thanks your suggestion.

[mwchang]

Dear Prof. Graeme

Thanks for your previous suggestions, but my ciNEB calculations are still not converged.

For example:

This is the vef.pl result for image 01, as you can see the forces are still oscillating.

opt cycle force total energy

0 0.63030900 -145.934412 0
1 5.51628000 -145.772893 0.161519
2 0.43178700 -145.936477 -0.002065
3 0.47286000 -145.937322 -0.00291
4 0.47091500 -145.937718 -0.003306
5 0.33457600 -145.938485 -0.004073
6 0.35289000 -145.938611 -0.004199
7 0.54625000 -145.938286 -0.003874
8 0.64546000 -145.938624 -0.004212
9 0.49199900 -145.940069 -0.005657
10 0.21227700 -145.941297 -0.006885
11 0.26330000 -145.941394 -0.006982
12 0.70791200 -145.939655 -0.005243
13 1.05894100 -145.936763 -0.002351



Because you have mentioned if the angles near 90 degrees have a major problem,
I noted, in my calculation, the angle between images are at near 110 degree .


im01: NEB: distance to prev, next image, angle between 2.076460 0.986312 78.128055
im02: NEB: distance to prev, next image, angle between 0.986312 0.984101 115.674672
im03: NEB: distance to prev, next image, angle between 0.984101 1.020674 108.295474
im04: NEB: distance to prev, next image, angle between 1.020674 1.060602 128.359177
im05: NEB: distance to prev, next image, angle between 1.060602 1.079651 116.105476
im06: NEB: distance to prev, next image, angle between 1.079651 1.093867 129.759375



The result means that should I add some images ??


This is my first time doing NEB calculation and I quite lack experience in the work.

Could you give me some suggestions? I am thankful for your help.

[graeme]

Small angles between images usually indicates that there is poor resolution along the band. See if there are any intermediate minima developing. If there are, you can break up the band into separate calculations. If the band is just a long single-barrier path, you may need more images to sufficiently resolve the pathway. The NEB gives a peicewise linear approximation to the path, and if you are seeing angle approaching 90 degrees between images, the linear approximation is breaking down. But again, first just check the geometries along your path and the energy to see if there is anything strange or if there are any intermediate minima.

[mwchang]

Dear Prof. Graeme

I very appreciate your technical supports and quickly reply.
Yesterday, I added more images from six to twelve in the band.



However, I have some questions about the angle.

you said "Once angle approaching 90 degrees, the linear approximation is breaking down"

I am confused on the reason and how is the angle defined .

And another question is what are strange thins?

Could you give me some examples or related materials?



Sincerely

[graeme]

The angle is simply defined by cos(theta) = (v_i.v_(i+1)) where v_i and v_(i+1) are unit vectors along path segments i and i+1.

By strange, I am talking about a path which is not smooth (discontinuities in geometries), intermediate minima, or any other technical problems that might be present in the calculation.

[mwchang]

Dear Prof. Graeme

Thanks for your reply, again.

I learn a lot from you.

[mishralu]

Dear All,
I am trying to use force based optimizers- using IOPT =1 for LBFGS with these following settings-
ISTART = 0
ICHARG = 2 charge: 1-file 2-atom 10-const
ISPIN = 2
MAGMOM = 40*2.0 24*0
LASPH = T aspherical Exc in radial PAW
ICHAIN =0
IMAGES =5
IOPT =1
LCLIMB = .True.
Dispersion Interaction
LVDW = .TRUE. ## add a semi-empirical dispersion
Electronic Relaxation 1
ENCUT = 450.0 eV 36.75 Ry 6.06 a.u. 6.65 6.65 6.65*2*pi/ulx,y,z
EDIFF = 0.1E-07 stopping-criterion for ELM
Ionic Relaxation
EDIFFG = -0.01 stopping-criterion for IOM
NSW = 150 number of steps for IOM
NBLOCK = 1; KBLOCK = 500 inner block; outer block
IBRION = 3 ionic relax: 0-MD 1-quasi-New 2-CG
ISIF = 2 stress and relaxation of ions
POTIM = 0 time-step for ionic-motion
INCURV = 0.001
ISMEAR = -5
SIGMA = 0.1
GGA_COMPAT = .FALSE
----------------------------
However, after even 150 ionic steps contcar of all the images are same as the POSCARs and in all CONTCAR files, below the coordinates I am seeing Nan Nan inn place of forces....
Kindly correct me if I am doing any mistake.
Thanks,
Abhishek