CI-NEB forces are fluctuating and calculations cannot converge

[youmukk]

Dear vtst developers,

I am working on an alloy system on which a carbon dioxide will absorb and then dissociate into a carbon monoxide and an oxygen atom. I am calculating the ernergy barrier but some problems appeared.

In the step1, I just used the one Gamma point and IBRION = 3;LCLIMB = .TRUE.;POTIM = 0.05 to calcualte in order to get resonable structures from the structures which were gotten from the commad nebmake.pl 3 .

In the next step, I increased the number of the K-points to a 2 2 1 Gamma point and used settings IOPT = 3;MAXMOVE = 0.15 to make the forces lower than 0.2 eV/Angstrom and it successfully converged.

Prombles appeared when I tried to use IOPT = 1 to get the forces lower than 0.05 eV/Angstrom, the forces are just fluctuating. Here I put the foreces on image 02. So I stopped it after 131 ionic steps. Files about this step are put in the direcotry named step3.
0 0.17411400 -536.945892 0
1 0.17491400 -536.946479 -0.00058726
2 0.18157700 -536.947340 -0.0014478
3 0.23383400 -536.952473 -0.00658155
4 5.83524900 -536.521717 0.424175
5 5.66401800 -536.617289 0.328603
6 5.28148000 -536.672328 0.273563
7 3.66198900 -536.845838 0.100053
8 1.67076400 -536.943098 0.00279359
9 0.14256800 -536.963727 -0.0178355
10 0.20828000 -536.963739 -0.0178468
11 0.15389100 -536.964086 -0.0181938
12 0.13894300 -536.964420 -0.0185285
13 0.16229700 -536.966445 -0.0205534
14 0.18500500 -536.969797 -0.023905
15 0.31098400 -536.980132 -0.0342407
16 0.30113600 -536.977357 -0.0314654
17 0.16741400 -536.975224 -0.0293327
18 0.20763400 -536.979553 -0.0336613
19 26.86745200 -535.103266 1.84263
20 3.70501400 -536.875373 0.0705183
21 2.27935300 -536.956139 -0.0102473
22 0.86716000 -536.996553 -0.0506613
23 0.56777000 -537.000369 -0.0544778
24 0.51003600 -537.008801 -0.0629093
25 1.31061800 -537.018274 -0.0723824
26 1.58790300 -537.058045 -0.112153
27 2.98917600 -536.959945 -0.0140531
28 1.65661700 -537.012482 -0.0665902
29 1.54790500 -537.021760 -0.0758685
30 2.37080800 -537.004216 -0.0583239
31 1.04223200 -537.029459 -0.0835669
32 1.04360100 -537.028095 -0.0822028
33 1.01029500 -537.023383 -0.0774914
34 1.00135300 -537.013120 -0.0672284
35 1.61405700 -536.941378 0.00451361
36 7.52282600 -536.674375 0.271517
37 4.92694000 -536.773387 0.172505
38 7.34990100 -535.857871 1.08802
39 5.59699700 -536.488304 0.457588
40 9.05165300 -536.561136 0.384756
41 2.47379100 -536.873097 0.0727943
42 1.00816500 -536.916857 0.0290347
43 1.10822900 -536.919897 0.0259943
44 1.64434800 -536.919908 0.0259836
45 1.48596400 -536.961691 -0.0157988
46 6.33228600 -535.874579 1.07131
47 5.93237900 -536.484412 0.46148
48 1.75246200 -536.968354 -0.0224624
49 0.97749500 -536.982581 -0.0366895
50 0.78722400 -536.990316 -0.0444246
51 0.85580000 -536.990435 -0.0445431
52 0.34432500 -536.997005 -0.0511135
53 1.60999500 -536.983088 -0.0371961
54 1.62099700 -536.994239 -0.0483469
55 0.42144400 -537.010217 -0.0643248
56 0.31848700 -537.012411 -0.0665194
57 0.33142700 -537.013126 -0.0672343
58 0.31549800 -537.014218 -0.068326
59 0.22065600 -537.016933 -0.0710417
60 0.23221800 -537.018897 -0.0730057
61 0.31312200 -537.028718 -0.082826
62 5.39701000 -536.896434 0.0494579
63 5.21725500 -536.769399 0.176493
64 3.38678800 -536.965633 -0.0197411
65 0.68564000 -537.065103 -0.119211
66 0.66585800 -537.071690 -0.125798
67 0.92260700 -537.071860 -0.125968
68 0.68149200 -537.072086 -0.126194
69 0.96284000 -537.062292 -0.1164
70 4.85307300 -536.752636 0.193256
71 1.10238000 -537.017566 -0.0716742
72 1.18508000 -537.021120 -0.0752279
73 0.62146700 -537.027092 -0.0812005
74 0.51316600 -537.028110 -0.0822183
75 0.46314400 -537.028740 -0.0828478
76 0.47061700 -537.031414 -0.0855226
77 1.24568100 -537.090752 -0.144861
78 8.86721300 -536.746751 0.19914
79 1.61452400 -537.104838 -0.158947
80 1.57125200 -537.107865 -0.161974
81 6.29947400 -536.920899 0.0249928
82 0.86209600 -537.123159 -0.177268
83 0.71585900 -537.124307 -0.178415
84 0.59657900 -537.124110 -0.178218
85 0.55851800 -537.122950 -0.177059
86 0.52348900 -537.118222 -0.17233
87 0.50971400 -537.094787 -0.148896
88 2.82655100 -537.006939 -0.0610477
89 1.08903100 -537.049724 -0.103833
90 0.80770400 -537.053725 -0.107834
91 0.55592500 -537.058013 -0.112121
92 0.56389000 -537.056991 -0.111099
93 0.45175200 -537.055452 -0.10956
94 0.42265900 -537.054185 -0.108293
95 0.73622200 -537.064381 -0.118489
96 6.83103900 -536.926757 0.0191344
97 3.98828200 -537.059122 -0.11323
98 3.11288500 -537.089389 -0.143497
99 1.20911500 -537.139649 -0.193757
100 0.77002000 -537.152870 -0.206978
101 0.88696100 -537.167163 -0.221272
102 0.93384200 -537.194362 -0.248471
103 1.62295300 -537.301073 -0.355182
104 5.82995300 -537.120238 -0.174347
105 2.65369900 -537.250961 -0.305069
106 2.10483900 -537.280244 -0.334352
107 1.71262800 -537.288985 -0.343094
108 3.50275500 -537.197254 -0.251362
109 7.45272000 -536.930127 0.0157652
110 6.45986600 -536.812900 0.132992
111 7.77922000 -535.814851 1.13104
112 6.34800400 -536.460276 0.485616
113 36.41654500 -534.173863 2.77203
114 4.72613700 -536.920468 0.0254242
115 2.43813300 -536.987179 -0.041287
116 1.18667300 -537.005653 -0.059761
117 0.99927300 -536.997890 -0.0519981
118 0.97033900 -536.967339 -0.0214475
119 1.87833500 -536.954447 -0.00855508
120 18.27968000 -535.879819 1.06607
121 3.43426700 -537.010553 -0.0646617
122 2.46872700 -537.047869 -0.101977
123 1.41885900 -537.104804 -0.158912
124 1.76554500 -537.115350 -0.169458
125 2.77090300 -537.089520 -0.143628
126 39.99645500 -533.592169 3.35372
127 3.41709500 -537.077202 -0.13131
128 1.89462200 -537.127342 -0.18145
129 0.84635600 -537.153958 -0.208066
130 1.06260800 -537.151974 -0.206082
131 2.09870300 -537.114693 -0.168801

In step4, I tried to calculate by using method IOPT = 2 because I read Professor Graeme`s comments about IBRION = 1/IOPT = 1 ,saying 'For IOPT=1, it is very important to set the INVCURV variable sufficiently small so as to ensure stable optimization.' However, forces were also spiking above 1 or even 10 eV/Ang.

I have attached the tar file of my calculations.I would be grateful if you could help me resolve this.

In addtion, I have done a Dimer calculation, although it did not converge, the vibrational analysis showed it only has one imaginary mode. I have attached the files about the Dimer calculation named dimer.tar.gz

I really do not know what to do now, can you give me some advice ?

Thanks!

[graeme]

In the attached calculation, I continued your step2 using the same optimizer (IOPT=3) and it converged fairly rapidly down to 0.1 eV/Ang. I will repeat this with IOPT=1 and find a suitable setting - I expect that you need to reduce INVCURV.

The problem with the dimer calculation is that you are using a binary there which does not have the VTST code linked into it. You can just grep the OUTCAR for VTST to check this. It is a little puzzling as your NEB calculation does, so for some reason you used a different binary for the dimer calculation. Anyway, I will also run a dimer calculation and make sure that it converges.

[youmukk]

Thank you, Professor Graeme! Can you give some advice about how to choose the different IOPT settings? Should I continue your calculation with the same setting, that is IOPT = 3, or maybe I need to try IOPT = 2 and EDIFF = 1E-7, for it is near the saddle point.

And I am also a little puzzled about the parameters settings with different optimizers, for example, when performing a calculation with IOPT = 3, how should I adjust the MAXMOVE and the TIMESTEP parameters if the forces don not converge? Should I go to the papers to see how it works, but for it seems a little hard for me.

[graeme]

Optimization is a bit of an art. But first, you can use any optimizer that will reach your force criterion. IOPT=3 is a safe and stable optimizer and if you can reach convergence with it (as we did) there is no need to use a different one.

That said, IOPT=3 and 7 are first-order optimizers in that they only use forces. Second order optimizers, including IOPT=1 and 2, make use of curvature information and will, in principle, converge faster, particularly near stationary points. You will need to increase the precision of the forces for these second order optimizers, as you mention, using (for example) EFIFF=1e-7. Note that these second-order optimizers are a little more finicky - there is almost always a tradeoff between stability and efficiency. If you have trouble with the second order methods, you can go back to IOPT=3.

In terms of the IOPT=3 settings, the primary control is with the TIMESTEP variable. If the value is too small, convergence will be slow. The rate of convergence will increase as the TIMESTEP value is increased until the method becomes unstable. A value of 0.1 is generally safe, but the actual value will depend upon the stiffest vibrational frequencies in your system. The MAXMOVE value can be safely set to something like 0.1 - 0.5 Ang and left alone - it is really just a sanity parameter to make sure that no optimizer moves the system more than a sensible trust-radius.

[youmukk]

Professor Graeme,
Thanks for your kindness and patience!
Now I think I know more about how to adjust the art of optimization. Simle :).

I have another two questions here.
The first one is that when we are performing a NEB calculation, the scripts neb.pl and vef.pl will give us, respectively, the current information of all images like forces and a list of the force and energy at each iteration. Is there a general rule which we can use to determine if a NEB calculation is good or bad. As far as I am concerned, I can check the structures to see if they go away and I can also check the forces and see if one of them goes too high like above 10 eV/Ang. But I am confused about the upper limit of force fluctuation, imagine a NEB calculation is running and is it safe that the forces fluctuate below 0.5 eV/Ang？Or if the force is higher than some value we should stop and go back the settings. Maybe I should just run about 200 steps and see the trend ?

The other question is about DIMER method. I performed a DIMER calculation on the same alloy but different adsorption site.
The settings are IOPT = 2;DRotMax = 4;EDIFF = 1E-8;EDIFFG = -0.05;NSW = 300. In the case of first run, call it run1, although the Force and Torque are high, the energy continues to go down. After the ionic steps reached 300, I copied CONTCAR to POSCAR and NEWMODECAR to MODECAR to continue the calculation but this time energy began to increase and the Curvature at last was positive.
Should I try IOPT = 3 or do something else? I have attached the files.