H2O drifting away during NEB

[Soon Ho Kwon]

Hello, prof. Graeme

I'm struggling with proton transfer from neutral water to Cu(111). (Volmer reaction)
I have 3 H2O in the system and the final state would be 2H2O(l)+OH*+H*.
During the NEB, H2O drift away from the surface making high forces as follows.

Step max_force barrier image delta_barrier
0 3.05694 eV/A 1.66845 eV 03 1.67
1 2.55418 eV/A 1.57188 eV 03 -0.0966
2 2.22045 eV/A 1.57375 eV 03 0.00187
3 1.55696 eV/A 1.53141 eV 03 -0.0423
4 6.27581 eV/A 1.49808 eV 03 -0.0333
5 1.83588 eV/A 1.50562 eV 03 0.00754
6 2.50277 eV/A 1.50235 eV 03 -0.00327
7 2.49879 eV/A 1.45564 eV 03 -0.0467
8 4.09923 eV/A 1.42517 eV 03 -0.0305
9 3.45421 eV/A 1.40589 eV 03 -0.0193
10 2.17794 eV/A 1.39411 eV 03 -0.0118
11 0.98926 eV/A 1.36637 eV 03 -0.0277
12 1.30494 eV/A 1.35983 eV 03 -0.00655
13 0.96488 eV/A 1.35131 eV 03 -0.00851
14 0.78059 eV/A 1.34058 eV 03 -0.0107
15 3.02258 eV/A 1.37906 eV 03 0.0385
16 1.15891 eV/A 1.34412 eV 03 -0.0349
17 0.84033 eV/A 1.33545 eV 03 -0.00868
18 0.75162 eV/A 1.31875 eV 03 -0.0167
19 0.74229 eV/A 1.3175 eV 03 -0.00125
20 0.7003 eV/A 1.30286 eV 03 -0.0146
21 0.8811 eV/A 1.29382 eV 03 -0.00904
22 0.81971 eV/A 1.28425 eV 03 -0.00957
23 1.01024 eV/A 1.26323 eV 03 -0.021
24 1.84678 eV/A 1.25703 eV 03 -0.0062
25 1.21607 eV/A 1.25321 eV 03 -0.00381
26 0.93729 eV/A 1.22658 eV 03 -0.0266
27 0.885 eV/A 1.18146 eV 03 -0.0451
28 0.84607 eV/A 1.17657 eV 03 -0.00489
29 0.99158 eV/A 1.17014 eV 03 -0.00643
30 0.9178 eV/A 1.16502 eV 03 -0.00512
31 0.75951 eV/A 1.1532 eV 03 -0.0118
32 1.15428 eV/A 1.14522 eV 03 -0.00798
33 1.03127 eV/A 1.14526 eV 03 3.83e-05
34 0.92635 eV/A 1.13295 eV 03 -0.0123
35 0.85759 eV/A 1.12702 eV 03 -0.00593
36 0.70499 eV/A 1.11811 eV 03 -0.00891
37 0.70231 eV/A 1.11296 eV 03 -0.00516
38 0.64163 eV/A 1.10212 eV 03 -0.0108
39 0.74924 eV/A 1.09918 eV 03 -0.00293
40 0.44481 eV/A 1.08917 eV 03 -0.01
41 0.57803 eV/A 1.08614 eV 03 -0.00303
42 0.39372 eV/A 1.07898 eV 03 -0.00717
43 0.49728 eV/A 1.0749 eV 03 -0.00408
44 0.4172 eV/A 1.06939 eV 03 -0.0055
45 0.52307 eV/A 1.06724 eV 03 -0.00215
46 0.52626 eV/A 1.06256 eV 03 -0.00468
47 0.60925 eV/A 1.05816 eV 03 -0.0044
48 0.61046 eV/A 1.0532 eV 03 -0.00496
49 0.61396 eV/A 1.0511 eV 03 -0.00211
50 0.54554 eV/A 1.04708 eV 03 -0.00401
51 0.55641 eV/A 1.04688 eV 03 -0.000207
52 0.43906 eV/A 1.04191 eV 03 -0.00496
53 0.44366 eV/A 1.03933 eV 03 -0.00258
54 0.40398 eV/A 1.03688 eV 03 -0.00245
55 0.58089 eV/A 1.03762 eV 03 0.000738
56 0.37569 eV/A 1.03341 eV 03 -0.00421
57 0.3964 eV/A 1.0307 eV 03 -0.00271
58 0.36282 eV/A 1.02942 eV 03 -0.00128
59 0.45649 eV/A 1.02735 eV 03 -0.00206
60 0.29281 eV/A 1.02591 eV 03 -0.00144
61 0.35017 eV/A 1.02343 eV 03 -0.00248
62 0.27517 eV/A 1.02118 eV 03 -0.00225
63 0.29774 eV/A 1.01955 eV 03 -0.00163
64 0.29806 eV/A 1.01636 eV 03 -0.00319
65 0.40621 eV/A 1.01347 eV 03 -0.00289
66 0.36251 eV/A 1.01016 eV 03 -0.00332
67 0.43622 eV/A 1.00464 eV 03 -0.00552
68 0.45509 eV/A 1.00119 eV 03 -0.00345
69 0.63143 eV/A 0.99717 eV 03 -0.00402
70 0.5466 eV/A 0.99434 eV 03 -0.00283
71 0.63207 eV/A 0.98774 eV 03 -0.0066
72 1.34444 eV/A 0.97984 eV 03 -0.0079
73 0.69862 eV/A 0.98446 eV 03 0.00462
74 0.68705 eV/A 0.98131 eV 03 -0.00315
75 0.66958 eV/A 0.97749 eV 03 -0.00382
76 0.55903 eV/A 0.9765 eV 03 -0.000994
77 0.90746 eV/A 0.97769 eV 03 0.00119
78 0.43758 eV/A 0.96874 eV 03 -0.00895
79 0.43727 eV/A 0.96628 eV 03 -0.00246
80 0.3384 eV/A 0.96567 eV 03 -0.000612
81 1.63943 eV/A 0.96833 eV 02 0.00266
82 1.45558 eV/A 1.0287 eV 02 0.0604
83 1.70772 eV/A 1.05007 eV 02 0.0214
84 1.73478 eV/A 1.0574 eV 02 0.00732
85 1.8025 eV/A 1.07822 eV 02 0.0208
86 1.93782 eV/A 1.11749 eV 02 0.0393
87 2.1825 eV/A 1.18997 eV 02 0.0725
88 2.02991 eV/A 1.14479 eV 02 -0.0452
89 1.72578 eV/A 1.14054 eV 02 -0.00425
90 0.95223 eV/A 1.13709 eV 02 -0.00344
91 1.03883 eV/A 1.16826 eV 02 0.0312
92 0.93624 eV/A 1.16958 eV 02 0.00132
93 0.7085 eV/A 1.17541 eV 02 0.00583
94 1.12863 eV/A 1.21601 eV 02 0.0406
95 0.88914 eV/A 1.20862 eV 02 -0.00739
96 1.31276 eV/A 1.19668 eV 02 -0.0119
97 2.27431 eV/A 1.2042 eV 02 0.00753
98 1.89144 eV/A 1.19571 eV 02 -0.0085
99 1.71654 eV/A 1.18288 eV 02 -0.0128
100 1.22759 eV/A 1.15605 eV 02 -0.0268
101 1.88538 eV/A 1.15776 eV 02 0.00171
102 3.57721 eV/A 1.41577 eV 02 0.258
103 2.62545 eV/A 1.21486 eV 02 -0.201
104 2.50091 eV/A 1.20896 eV 02 -0.0059
105 2.1522 eV/A 1.21538 eV 02 0.00642
106 2.20083 eV/A 1.42999 eV 03 0.215
107 4.53661 eV/A 2.49761 eV 03 1.07
108 3.17172 eV/A 1.76435 eV 03 -0.733
109 3.24476 eV/A 1.88985 eV 03 0.125
110 3.87588 eV/A 2.42489 eV 03 0.535
111 11.92817 eV/A 4.28011 eV 03 1.86
112 4.57053 eV/A 2.84807 eV 03 -1.43
113 4.22857 eV/A 3.04511 eV 03 0.197
114 4.63462 eV/A 3.71301 eV 03 0.668
115 4.85503 eV/A 4.04243 eV 03 0.329
116 4.70873 eV/A 3.80435 eV 03 -0.238

I set -4 negative charge of the system on purpose.
(I attached the input & output files.)

It will be so great if I can get any advice to adjust those calculations.
Thank you,

Sincerely,
Soonho

[graeme]

Thanks for uplaoding the files. I see the problem and will adjust parameters to find a solution.

[Soon Ho Kwon]

Thank you so much for your prompt response.
I'm also trying to see if # of images (to reduce the distance below 1 A) or LBFGS optimizer can resolve this problem.

I hope you have a lovely day!

Sincerely,
Soonho

[graeme]

I did my standard calculation for a problematic NEB and the results seem ok (attached). I used somewhat cheaper settings, which I think are appropriate for your system. In the first calculation (run1) I accidentally used built-in vasp CG because you specified this (IBRION=2) earlier in the INCAR file. This is actually the only optimizer that you should not use for the NEB because it relies on the force being consistent with the energy, which is not true for the NEB. This can lead to early termination, as it did in run1.

The subsequent run with IBRION/IOPT=3 converged reasonably well within 200 steps. I see a primary barrier of 1.18 eV and then a small barrier associated with H diffusion to the final state site that you specified. I used 8 images because I was suspicious that this would be a two-step process, as it was. You should be able to easily increase the precision and reconverge, as needed.

[Soon Ho Kwon]

First of all, Thank you for your help.

I noticed that you relaxed the structure at the neutral net charge of the system for the initial state. Since I tried to mimic the large negative potential environment, I increased the total number of electron on purpose which evokes water desorption and makes it H pointing down. but, unfortunately, that might be one of the factors for the abrupt force changes during NEB.

I also tried neb with 8 images, and IBRION=3, IOPT=1 which showed there were two additional local minimums. (I might be able to assume the first maximum would be the highest barrier and further converge it with high criteria using the dimer method) But I don't understand why the forces are so fluctuating and how I can remove the variance. (please look at the list_barrier file)

Any suggestions and advice will be very helpful.
I appreciate your help.
Have a lovely day!

Sincerely,
Soonho

[graeme]

I did a conservative relaxation of the NEB for the charged system and it looks ok (attached).

You really have to be careful with IBRION=1 / IOPT=1. For IOPT=1, it is very important to set the INVCURV variable sufficiently small so as to ensure stable optimization. In your case, forces were spiking above 1 or even 10 eV/Ang, which is not ok.

Here, I used IOPT=3 and was able to get reasonable convergence (<0.1 eV/Ang) in 100 iterations. I think you should be able to use a similar approach to converge further and/or with more precise settings.

[Soon Ho Kwon]

Dear, prof. Graeme.

Your results obviously show that it's not a problem of the pathway but the optimizer I used. I might need to use QM (or FIRE) and test those variables.

Thank you so much for your help.
Have a wonderful day & weekend!

Sincerely,
Soonho