UT theoretical chemistry code forum

Dear vtst community,

I am using vtsttools to investigate the migration of oxygen in cubic ZrO2 but just met some problems. I think I need your kind help. There are basically two:

1. I just used one image to find out the saddle point and it converges well. The barrier is about 0.4 eV. But when I tried to manually move the oxygen off the saddle point a little bit alone the transition path (static calculation) I got energy increase! This does not make sense since cNEB "pushes" the oxygen opposite to the force alone transition path to the saddle point. Below is the test results (the displacement is from -0.4 to 0.4 angstrom)
OUTCAR.-0.4: energy without entropy= -894.43508677 energy(sigma->0) = -894.44102791
OUTCAR.-0.2: energy without entropy= -894.57379648 energy(sigma->0) = -894.57960596
OUTCAR.0.0: energy without entropy= -894.61576250 energy(sigma->0) = -894.62153739
OUTCAR.0.2: energy without entropy= -894.57379769 energy(sigma->0) = -894.57960717
OUTCAR.0.4: energy without entropy= -894.43509022 energy(sigma->0) = -894.44103134
As you can see the energy in OUTCAR.0.0 is the lowest.

2. So I thought maybe one image was not enough (From my understanding it shouldn't be the reason but just try). I increased the number of images to three. Then I have directories 00 01 02 03 04. Actually 02 has the POSCAR obtained from step 1. 01 and 03 were generated between 00 (initial) and 02, and between 02 and 04 (end), respectively, by script nebmake.pl. The calculation simply did not converge. Instead, it produced confusing results:
01/OSZICAR: 20 F= -.89470300E+03 E0= -.89467132E+03 d E =-.272434E-02
02/OSZICAR: 20 F= -.92416539E+03 E0= -.92412358E+03 d E =-.822467E-01
03/OSZICAR: 20 F= -.89437693E+03 E0= -.89434832E+03 d E =-.127992E-01
I just stopped the calculation at 20th step because there was no need to proceed. As you see, the energy of the middle one is far from expected as should be close the energy in OUTCAR.0.0. I do not know where is wrong. Every setting seems to be normal to me.

Any possible solutions or ideas are highly appreciated. Thank you so much.

My INCAR:

ISTART = 0
ICHARG = 2

PREC = Accurate
ENCUT = 480 eV
ALGO = N
NELM =60 ;
EDIFF = 1E-05
ISYM = 2
LMAXMIX= 4
NCORE = 8

ISIF = 2
EDIFFG =-1E-02
NSW = 60
IBRION = 3
POTIM = 0
IMAGES = 3
LCLIMB = .TRUE.
IOPT = 2

LREAL =.TRUE.
ISMEAR = 0 ; SIGMA = 0.08

Please post a .tar.gz file of the two calculations and I'll give advice. You can remove any CHG* and WAV* so that the files are not too large.

Dear Graeme,

Thank you very much for the reply. Please find the file I attached. Just two notes here.

Prob. 1 lies in the directory 02/test-force with the converged CONTCAR based on one-image cNEB calculation. The rest files correspond to Prob. 2.

Please let me know if any more information is needed.

H

You need to turn symmetry off (ISYM=0). The problem is that your 02 image detected symmetry and then it was broken during relaxation. With the geometry not having the symmetry of the wavefunction leads to the non-physical energies that you are seeing.

Thank you very much for your help. I will try it again.

H

[quote=graeme post_id=8900 time=1491923459 user_id=2]
You need to turn symmetry off (ISYM=0). The problem is that your 02 image detected symmetry and then it was broken during relaxation. With the geometry not having the symmetry of the wavefunction leads to the non-physical energies that you are seeing.
[/quote]

Dear Graeme,

I followed your suggestion to switch off the symmetry, the cubic structure transformed to another structure although I relaxed only the ionic positions, because the cubic one actually is not dynamically stable. I thus constrained several atomic layers (three 001 plane) using selective dynamics mode to prohibit the transformation, and even manually shift a little bit the migrating atom off the transition state (since from nebmake.pl output it was exactly at the middle of the cell, nominally the transition state). The neb calculation converged (ISYM=0), and the migrating atom was pushed back to the center. This was a good sign which indicated the transition state was indeed energetically the maximum. However, I tested the force on the migration atom at transition state again by manually shifting it off the saddle point along transition path (in 02/test-force as also described in the previous post), and the resulting energy still increases!! I am really confused by that. Would you please help me to check the files again? Thank you very much.

H

The configuration in OUTCAR.0.0 does not appear to be the same and in the 02/CONTCAR file, which I assume is your saddle point. There are also relatively high forces listed in OUTCAR.0.0.

[quote=graeme post_id=8937 time=1493151469 user_id=2]
The configuration in OUTCAR.0.0 does not appear to be the same and in the 02/CONTCAR file, which I assume is your saddle point. There are also relatively high forces listed in OUTCAR.0.0.
[/quote]

Dear Graeme,

Thank you very much for your quick response.

I think the position difference between OUTCAR.0.0 and 02/CONTCAR is due to that I used the below statement
zpos=`awk 'BEGIN{printf "%18.16f", '$disp'/10.3123638438128999 + 0.5}'`
in the job.sh script. $disp is the displacement I introduced to shift the migrating atom off the saddle point along the path (from -0.4 to 0.4 \AA). When $disp=0.0, the z-component is exactly 0.5, which is a little bit different from the original 0.5002580330854806 as in 02/CONTCAR (or POSCAR0). Since the displacement is relatively large, it is guaranteed that the displacements have overwhelmed this small difference and shifted the migrating atom off the real saddle point, and hence, if the converged structure is really the saddle point, the energy should drop. That is where I found hard to understand.

H

Double check the geometry in your test; the structure from your energy test is different from the CONTCAR in many degrees of freedom and it has a high force.

[quote=graeme post_id=8939 time=1493171071 user_id=2]
Double check the geometry in your test; the structure from your energy test is different from the CONTCAR in many degrees of freedom and it has a high force.
[/quote]

Dear Graeme,

I apologize for not checking it carefully before bothering you for that.

Yes, some forces are too high. I just compared the POSCARs generated by job.sh from POSCAR0 (identical to 02/CONTCAR). It is for sure that only the z component of the migrating atom had been modified by the script. I also transformed the POSCAR into Cartesian mode, and compared it with the positions in the corresponding OUTCAR. They are completely the same. So I think geometry is fine.

After that, I checked the forces in 02/OUTCAR, and I found that the forces did not converge (threshold is 0.01). I picked some positions here
2.52344 2.52344 6.46437 0.022702 0.022703 -0.041720 T T T
2.58478 7.72758 6.44918 -0.006882 0.006889 -0.001952 T T T
7.72758 2.58478 6.44918 0.006893 -0.006888 -0.001956 T T T
7.78893 7.78893 6.46437 -0.022703 -0.022703 -0.041724 T T T
0.00000 0.00000 3.86071 0.000002 0.000006 -0.000506 T T T
0.00000 5.15618 3.87023 -0.000004 -0.000008 -0.003460 T T T
5.15618 0.00000 3.87023 -0.000001 0.000005 -0.003462 T T T
5.15618 5.15618 5.25931 -0.000001 0.000001 -0.158547 T T T (the migrating atom)
0.00000 0.00000 6.45165 -0.000004 -0.000002 0.001997 T T T
0.00000 5.15618 6.44214 0.000003 0.000001 -0.002636 T T T
5.15618 0.00000 6.44214 -0.000005 -0.000004 -0.002642 T T T
2.52344 2.52344 3.84801 -0.020756 -0.020758 -0.033155 T T T
2.58477 7.72759 3.86317 0.005407 -0.005406 -0.001454 T T T
7.72759 2.58477 3.86317 -0.005416 0.005412 -0.001463 T T T
7.78892 7.78892 3.84801 0.020752 0.020757 -0.033158 T T T

So why did the calculation stop and report 'required accuracy is reached' while the forces had not been reduced to be small enough? I used quick min in this neb calculation. Does optimizer really matter so much? Besides, the z component in the coordinates of migrating atom is strange. 5.25931 is not consistent with CONTCAR, but the one in POSCAR (0.51 in fraction). Anywhere is wrong?

Best

H

Dear Graeme,

I switched to CG optimizer and continued the neb calculation. The calculation converged and the forces are OK now (Please see the attachment). I still tested the energies as I did in test-force/. The variation trend of energy due to displacement of migrating atom is the same, i.e., it increases with increasing displacement. Please refer to the 02/test-force/forces in which I picked the forces from OUTCAR.0.0 corresponding to atoms allowed to relax. They are small now. So, the problem remains. Really confusing.

H

Ah yes, my apologies, I understand now that you have a number of frozen atoms and those high forces are on the frozen atoms so they do not contribute to the convergence criterion.

Try making smaller displacements along your reaction coordinate. Then, also, try to minimize this displaced structure and see if it goes to your original minima, or back to the 'saddle' structure -- this is a way to check for an intermediate minimum. You can also use more NEB images to get a better resolution of the minimum energy path and again see if there is a intermediate minimum as well as a more accurate tangent to the path at the saddle.

Dear Graeme,

What you are recommending is exactly what I have tried to do.

The nebmake.pl generated 01 02 03, and the position of migrating atom in 02 was in the center of the cell due to the equivalence between 00/POSCAR and 04/POSCAR by symmetry. I then manually introduced displacement to the migrating atom along reaction path (in [001] direction), from z=0.5 to z=0.51. I used the OPT=3 (quick min) and calculation converged. The atom was also "pushed" back to the center of the cell, very close to z=0.5. This was why I called it "a good sign". The problem occurred when I tested the energy by introducing displacements in the same way but based on the converged saddle point structure (conventional static calculations). The energy should decrease due to the displacement, but it increased.

Since I used quick min, I then switched to CG and continued based on the previously converged structures. The calculation using CG also converged. Forces are fine (threshold even higher). However, testing energy had the same problem, as I attached the files.

Since the CG calculation was based on the previously converged structure, 02/POSCAR was too close to saddle point. I am trying again by manually shifting the atom off the saddle point and then see if it will be pushed back as quick min did, as you also recommended. I will report the results later. Thank you.

H

[quote=graeme post_id=8943 time=1493222352 user_id=2]
Ah yes, my apologies, I understand now that you have a number of frozen atoms and those high forces are on the frozen atoms so they do not contribute to the convergence criterion.

Try making smaller displacements along your reaction coordinate. Then, also, try to minimize this displaced structure and see if it goes to your original minima, or back to the 'saddle' structure -- this is a way to check for an intermediate minimum. You can also use more NEB images to get a better resolution of the minimum energy path and again see if there is a intermediate minimum as well as a more accurate tangent to the path at the saddle.
[/quote]

Dear Graeme,

I tried two ways:
(1) Introduce small displacement along reaction coordinate and use CG to do the neb calculation. The relaxation does not converge. The middle structure oscillates around the 'saddle' structure. However, when I increased spring constant from default to -7, it converged. According to my understanding of cNEB, the spring constant should not have effects on the point with highest energy. So why? The stranger thing is that when I manually displaced the migration atom along reaction coordinate (even as small as 0.02 \AA) using the 'saddle' structure and calculated the energy, the energy still increases compared to the original one. This conflicts with the cNEB results since force-based 'climbing' indeed happened. It is hard to understand.
(2) I used also more NEB images around the 'saddle' point. The neb calculation also converged, but testing energies using the 'saddle' structure get similar results to (1).

Could you find some possible reasons for that?

H