Convergence and rotational steps after CENTCAR is written

[sebagg]

Dear all,
I am having some problems understanding the structure of a dimer calculation, in particular about the convergence criteria, when the CENTCAR is written and specially why there seem to be rotational steps after the last CENTCAR (point 4 below). In the attached example the last step only says " Dimer: Projection", previous X rotational steps appear as " Dimer: Itr 13 Rot X \n Dimer: Rotation", before that the central step is "Dimer: Itr 13 Rot 0 \n Dimer: Central Point" and before that there is another " Dimer: Projection", so my questions are:

1) When is the CENTCAR file written? in "Projection" or "Central Point"? Is there some indication in the OUTCAR of when the CENTCAR file is written?
The entries in DIMCAR correspond to "Projection steps"? Does the atomic forces ("TOTAL-FORCE (eV/Angst)" block) directly above a "Dimer: Projection" file correspond to the coordinates written in the converged CENTCAR file?

2) Other of my calculations that end up without rotations in the last steps only say "Dimer: Central Point", does it indicates the positions and atomic forces in the last "TOTAL-FORCE (eV/Angst)" block directly correspond to the CENTCAR file? CONTCAR and CENTCAR files are the same in those cases.

3) According to the webpage, the dimer converges when the FORCE in DIMCAR is below the selected EDIFFG. This is not the case in the attached example but it ended saying it reached convergence. All atomic forces in the last "TOTAL-FORCE (eV/Angst)" are below EDIFFG as if that was the criteria used to claim convergence. What is used as criteria at the end?

4) I know from this forum to always use the CENTCAR file to continue dimers or for frequencies of a converged dimmer and it worked fine for the attached transition state but something I don't understand is that the CENTCAR file was written minutes before the CONTCAR, they have different geometries. According to XDATCAR the CENTCAR was the 68th geometry and CONTCAR is the last one(74th). It seems like the last geometry written was the CONTCAR file and if this correspond to the last "Dimmer: Projection" after Iteration 13, the CONTCAR geometry and not the CENTCAR geometry would have the lower atomic (and tangential) forces. Why is that?

5) are WAVECAR and CHGCAR along with every step (with the default flags in the INCAR)? in my example where it seems there are steps after the CENTCAR file is written, the WAVECAR would not correspond to the CENTCAR. In this case the CENTCAR should not be used with the WAVECAR for a frequency calculation.

Thanks and best regards,
Sebastian.

[graeme]

Hi Sebastian. I'll do my best to answer your questions.

1. The CENTCAR is written in "Dimer: Central Point"

2. Yes, that is exactly correct. When the dimer translations, we start by calculating the central point and write it to the CENTAR (in the vtstcode). But it is also written (by vasp) to the CONTCAR which will always have the latest point calculated. If the dimer does rotations, the CONTCAR will follow those, but the CENTCAR will remain at the center of the dimer.

3. I don't get this question. When I look at the forces at the end of your calculation I see:
FORCES: max atom, RMS 0.031290 0.007365
FORCES: max atom, RMS 0.042804 0.011165
FORCES: max atom, RMS 0.019824 0.005786
where convergence is reached (max atom force < 0.02 eV/Ang) in the final iteration.

4. I believe that step 14 is not a real step - it is just the vtstcode figuring that the calculation is converged. Oh, actually, I think this will answer your previous question properly. So in step 13, the CENTCAR was written and the force was below 0.02 eV/Ang. But the dimer method also rotates to lower the torque below 1 eV/Ang^2. In your INCAR you have the maximum rotations set to 2, so it did those and then decided that it was converged. In principle the rotations were not required, but that is why the CONTCAR is different from the CENTCAR - it was from step 13. I do not believe any calculations were done in step 14.

5. Yes, that is right, if there is rotation on the last step, the WAVECAR and CHGCAR will not be consistent with the CENTCAR.

Graeme

[sebagg]

Dear Prof. Greame, thanks for your quick answer, it helped me a lot with some of my questions but I still struggle to understand the final steps in my calculation. I would greatly appreciate if you could help me clarify some last things:

The final geometries(#X), relevant "Dimer" messages and projected forces are:

(#69) Dimer: Itr 13 Rot 0, Dimer: Central Point, Projection -- FORCES: max atom, RMS 0.029270 0.007231 <---- Here CENTCAR is written
(#70) Dimer: Itr 13 Rot 1, Dimer: Projection, Trial Rotation --- FORCES: max atom, RMS 0.031221 0.007335
(#71) Dimer: Itr 13 Rot 1 , Dimer: Rotation, Projection --------- FORCES: max atom, RMS 0.042654 0.011097
(#72) Dimer: Itr 13 Rot 2 , Dimer: Projection,Trial Rotation --- FORCES: max atom, RMS 0.031290 0.007365
(#73) Dimer: Itr 13 Rot 2 , Dimer: Rotation, Projection, RotNum > RotMax ---- FORCES: max atom, RMS 0.042804 0.011165
(#74) Dimer: Projection --------------------------------------------------------------------------- FORCES: max atom, RMS 0.019824 0.005786

1) If I understand it correctly, CENTCAR is written for the geometry #69 (Dimer: Central Point), geometries #70-#71 and #72-#73 are the dimer rotations used to adjust the dimer direction and #74 would be the new Central Point (and CENTCAR) if the convergence criteria was not meet. Is this correct?

2) The final geometry #74 is the first to meet the convergence criteria (EDIFFG<-0.02) and should be the best approximation to the transition state, but it correspond to the CONTCAR geometry while the CENTCAR geometry is still from geometry #69 and has larger "FORCES: max atom". Why is that?

3) This also means that the last "energy(sigma->0)" is not the energy of the central image. For energy barrier calculations I should be using the "energy(sigma->0)" right above the last "Dimer: Central Point" instead of the energy from the last one image since it may correspond to a rotation step, right?

Thanks and best regards,
Sebastian

[graeme]

Hm, that doesn't see quite right. Let me look into this more carefully and I'll let you know.

[graeme]

Ok, so what I said previously about the the CENTCAR vs CONTCAR was not correct. Below is a table of the calculation. The CENTCAR was last written at vasp iteration 69, dimer step 13. Then, after the final translation after dimer step 13 rotation 2, the force dropped below 0.02 eV/Ang and the calculation was terminated. Thus, the final CONTCAR is actually the best estimate of the saddle and this structure is consistent with the CHG* and WAV*. I think that I can fix the situation by writing the CENTCAR when we detecct convergence. Thank you for getting me to look into this more carefully!

------------
vasp vasp dimer dimer dimer
itr force step, rot task force torque comment

1 0.046393 1, 0 Dimer: Central Point 0.173853 write CENTCAR
2 0.046589 1, 1 Dimer: Trial Rotation 0.913984
3 0.118663 1, 1 Dimer: Rotation FN<FNMax so stop rotation
4 0.040270 OPT: CG step
5 0.088043 2, 0 Dimer: Central Point 0.162017 write CENTCAR
6 0.088535 2, 1 Dimer: Trial Rotation 0.986362
7 0.120989 2, 1 Dimer: Rotation FN<FNMax so stop rotation
8 0.031425 OPT: CG step
9 0.087184 3, 0 Dimer: Central Point 0.167795 write CENTCAR
10 0.088987 3, 1 Dimer: Trial Rotation 1.132475
11 0.104930 3, 1 Dimer: Rotation FN>FNMax so rotation
12 0.089367 3, 2 Dimer: Trial Rotation 1.04071
13 0.174345 3, 2 Dimer: Rotation RotNum > RotMax, stop rotating
14 0.077705 OPT: CG step
…
69 0.029270 13, 0 Dimer: Central Point 0.059627 write CENTCAR
70 0.031221 13, 1 Dimer: Trial Rotation 1.28922
71 0.042654 13, 1 Dimer: Rotation FN>FNMax so rotation
72 0.031290 13, 2 Dimer: Trial Rotation 1.343055
73 0.042804 13, 2 Dimer: Rotation RotNum > RotMax, stop rotating
74 0.019824 14 F<|EDIFFG|, stop
-----------

<note that the attached data is formatted better>

[graeme]

If you want to try a small patch, you can change the following in chain.F, around line 370 depending upon the version:

IF (ICHAIN==2) THEN
CALL dimer_fin()
ENDIF

to

IF (ICHAIN==2) THEN
CALL dimer_fin(posion,a,b)
ENDIF

and in dimer.F around line 195:

SUBROUTINE Dimer_Fin(POSION, LATT_A, LATT_B)
IF (IU6>=0) WRITE(IUdim,'(I5,F14.5,7x,A3,4x,F14.5,7x,A3,11x,A3,4x)') &
& Itr,F0r,'---',U0,'---','---'
END SUBROUTINE Dimer_Fin

to

SUBROUTINE Dimer_Fin(POSION, LATT_A, LATT_B)
REAL(q),DIMENSION(3,NIONS) :: POSION
REAL(q),DIMENSION(3,3) :: LATT_A, LATT_B

IF (IU6>=0) WRITE(IUdim,'(I5,F14.5,7x,A3,4x,F14.5,7x,A3,11x,A3,4x)') &
& Itr,F0r,'---',U0,'---','---'
! Write CENTCAR
POSION = R0
CALL KARDIR(NIONS, POSION, LATT_B)
IF (IU6>=0) THEN
OPEN(UNIT=IUcent, FILE='CENTCAR', STATUS='UNKNOWN')
CALL OUTCENT(IUcent, .TRUE., DINFO%SZNAM2, DINFO, 1._q, A, DINFO%NTYP, &
& DINFO%NITYP, DINFO%LSDYN, DINFO%NIONS, POSION, DINFO%LSFOR)
CLOSE(IUcent)
END IF
END SUBROUTINE Dimer_Fin

[sebagg]

Prof. Graeme, thank you very much for looking into this, for taking the time to explain it in detail my example and for the patch. I see that the small patch is already included in the main folder of the last downloadable version (vtstcode-184.tgz), and I will be talking to my cluster administrators to try it out.

Thanks again, and best regards,
Sebastian

[graeme]

Ok, but I should warn you that that patch is not for these changes. I could make a patch if you want. It's just I saw you were using 5.4.x and the current code is for vasp.6.2.x. I hope that with such simple changes, someone could just make them in the source code.