Hellmann–Feynman forces do not converge

[jeffrey_zom]

Hello,

I am doing a geometric convergence calculation for an electron polaron in a LiVO3 supercell (80 atoms, 545 electrons)
I have found that after a few dozen ionic steps, both the average force and the maximum force stay nearly constant.
However, the maximum force of ca. 0.025 eV/A is still above the convergence threshold I put of 0.01 eV/A.
I've done exactly the same calculations using the same settings in the INCAR for MVO3 (M=Cs/Rb/K/Na) and had no trouble.

Things I have tried:
- Increase the FFT grid density to ROPT = 3*-1E-4 (from -2.5E-4 normally)
This significantly reduced the total drift force, but did not cause lowering of the forces, see figure OUTCAR_1.png (purple is default ROPT for Accurate, yellow is increased)
Note: In these figures, the dots represent the maximum force, the straight line the average, the striped line the drift force. The red line is the convergence criterion.

- Changed IBRION from 1 to 2.
IBRION 2 performed much worse. The average and maximum forces increased again. see figure OUTCAR_2.png (purple is IBRION 2, yellow is IBRION 1)

- Increased EDIFF from 1E-6 to 1E-8
This did not change the forces (not shown)


I have noticed that the forces for most atoms stay constant (see OUTCAR_3.png in the .zip). I am puzzeled, as a non-zero force should result in some displacement of the atom, thereby changing the total force. A non-changing force I think suggests that either the displacement is negligibly small to result in no change in force, or that it is in some potential well and constantly jumps to the other side with the same force but different vector components, oscillating between these positions. I have also found that the atoms for which these forces stay high (and constant) are those that at the very border of my unit cell, i.e. some fractional coordinate is near equal to 0 or 1.

Could someone give advice on how to still achieve the convergence criterion?

The POSCAR of the relaxed supercell (so far as possible) is attached, as well as one of the INCAR/POTCAR/KPOINTS used (IBRION=1, ROPT = 3*-1E-4).

[henrique_miranda]

Very interesting question. And thank you for sharing the results of your careful testing.

I would like to suggest, if you can, to try and use the reciprocal space projection scheme LREAL=.FALSE.
This might be more accurate for computing forces.

Another INCAR tag that might help to reduce the noise in your forces is ADDGRID=.TRUE.

[jeffrey_zom]

Hi Miranda,

Thank you for your quick response.

I am currently running 2 additional calculations using LREAL=.FALSE. and ADDGRID=.TRUE., as you suggested.
I've not updated the other incar settings, but perhaps a combination of adaptions are needed (e.g. both lowering of the value of ROPT and addition of support grid).

If I read your response correctly, you think the noise, partially expressed in the result of a drift force, is the problem for failure of convergence here?
I will post an update here once the additional calculations have run for some time.

[henrique_miranda]

Note that when you set LREAL=False then ROPT is not used anymore (so in a sense it's one less parameter to tune).

I don't know what is the reason that you get a max force that never disappears. But we will try to find out:
It can be due to the projection scheme (real or reciprocal space) or due to insufficient ENCUT, but looking at your INCAR I doubt that (you use ENCUT=600 while ENMAX=400 from the POTCAR).
ADDGRID may or may not help, but it is probably worth a try.

Another possible source of error is the additional charge in your system (NELECT), you might need to do some electrostatics correction.
You can check if this is the case by not setting NELECT and seeing if in that case you reach the ionic minimization criteria.

The reason the ionic minimization does not converge is that the forces are not getting to be below the EDIFFG criteria, as you show in your plots.
But as you say, they should when you move the ions along their directions.
The reason they don't is probably related to numeric noise, I would say.

[jeffrey_zom]

Hi,

A quick update on the calculations:

I've run 3 additional calculations using (separetly) LREAL=.FALSE., EDIFF=1E-9 combined with ROPT=3*-1E-4, ADDGRID=.TRUE. .
Those 3 calculations in addition to the 2 performed earlier (IBRION=2 and ROPT=1E-4) are plotted in the attached figure.

For clearity:
3: IBRION = 2
4: ROPT = 1E-4
5: LREAL =.FALSE.
6: EDIFF = 1E-9, ROPT=3*-1E-4
7: ADDGRID =.TRUE.

Dots: Max force
Line: Average force
Dashed: Drift force

Unfortanetly, all of these changes to the INCAR did not result in an improvement of the convergence.
While the drift force does seem to change signficantly, the average and maximum forces do not.

Do you have any idea on other things I could try?
Electrostatics corrections (LDIPOL = .TRUE.?) are only available for cubic cells if I understand correctly, which is not the case here.

[henrique_miranda]

Ok, thank you for doing these tests.

I think we can safely conclude that the issue is not numerical noise in the forces.
Both LREAL=Auto or LREAL=FALSE yield very similar results, as is to be expected if the calculation is sufficiently converged.
ADDGRID seems to worsen the results, that is why we recommend users to test before using it.

Now, I would say that the problem is due to the need for electrostatics corrections.
But to be sure about that, it would perhaps be good if you could first run the calculation without the NELECT tag (meaning a neutral system).

[henrique_miranda]

Another INCAR tag that might be leading to noise in your forces is the ISMEAR=-5 which uses the tetrahedron method.
Perhaps you can try running the calculation with ISMEAR=0.

[jeffrey_zom]

Hi Miranda,

I've run additional tests, following your suggestions:

1) I've done calculations using ISMEAR=0 for the electron polaron, continueing from the almost converged POSCAR.

8: ISMEAR=0, EDIFF=1E-6,
9: ISMEAR=0, EDIFF=1E-9, ROPT=3*1E-4

The results in terms of forces against ionic step are shown in OUTCAR_1.png.
As you can see, the smearing method does not lead to improved convergence.

2) I've also looked if the additional electron is the bottleneck for convergence.

I have removed NELECT from the INCAR so that the system is neutral again.
I used the initial settings for these calculations (ISMEAR=-5, EDIFF=1E-6)
See OUTCAR_2.png

3: I've used the POSCAR of the almost converged geometry for the electron-polaron, therefore, the starting forces for some atoms were relatively large.
In that case, the geometry could not relax to the ground-state uncharged geometry, so that leads me to believe that the convergence problem does not have to do with the added electron.
4: Confusingly enough, when I use the exact same INCAR, only changing the POSCAR to that of the supercell constructed of the neutral fully converged (EDIFFG<0.01) unit cell, convergence is reached after 2 ionic steps. Ofcourse, this fast convergence makes sense from the perspective that this supercell constructed of the neutral unit cell should be identical. It does not make sense that a slightly different initial POSCAR (which has the exact same cell shape/size) can not reach convergence.

Any suggestions on further tests to solve the convergence problem are very welcome.

[henrique_miranda]

Firstly: thank you for your prompt responses and careful testing following my suggestions.

So far, we have excluded that this problem is due to:
1. sufficient ENCUT or using ADDGRID
2. the real or reciprocal projection scheme onto the PAW projectors
3. smearing method (ISMEAR=-5)
4. or due to the inclusion or not of electrostatic corrections (the problem remains in the neutral system)

At this point, I am just as puzzled as you, and unfortunately I don't have any further suggestion.
I will discuss with my colleagues, perhaps someone else has another idea.
I will try and play around with the input files you shared and let you know what I find.

[jeffrey_zom]

Hi Miranda,

Thank you for your help so far.
I will try to do some further testing, If anything interesting comes from these test, I will post an update here.
It would be great if one of your colleagues has an idea of the origin of the problem and/or how to solve it.

For if you want to do testing:
I have added the POSCAR for the (almost) converged electron-polaron,
the neutral supercell constructed from the converged unit cell,
and the initial POSCAR for the electron polaron (4 V-O bonds surrounding a single V atom are extended by 0.1A),

Kind regards,

Jeffrey

[henrique_miranda]

Thank you for sharing the files.
This calculation takes quite a while to complete, which makes it difficult to test slightly different INCARs to track down what the issue might be.
I tried some of your INCAR files locally, but before that I did a few changes
1. Gamma only calculation (just performance)
2. ENCUT=520 (just performance)
3. ISMEAR=0 (already tested)
4. EDIFF=1E-8 (if EDIFF is too large, then it is difficult to reach EDIFFG, but you already tested this before)
5. ALGO=Normal (results might differ)
6. remove NPAR=4 and KPAR=4 and set NCORE=4 (just performance)
7. remove ICHARG and ISTART (probably not harmful, but just to be sure)

Like this, the calculation runs faster.
In my case, I was able to converge to the threshold that you specified.
Could you try one of your problematic calculations with these changes and let me know if that works for you?

[jeffrey_zom]

Hi Miranda,

If I understood correctly, you changed KPOINTS, ENCUT and NCORE to increase the speed of the calculation, right?
(I used ENCUT = 650 as ENMAX = 499 (Li_sv)* 1.3 ~= 650)
And you changed ISMEAR from -5 to 0 because the tetrahedron method is not available for less than 4 kpoints, correct?
So the only significant change you made is to use ALGO=Normal?


I don't know if you used the additional electron.

I have performed two similar calculations (with NELECT = 545):
1) One that should be identical to yours. This calculation does not appear to converge, see OUTCAR_1.png. I also see that it results in a different CONTCAR, namely, the electron is now delocalised over the system, the different (polaron) geometry surrounding vanadium atom 18 is removed.

2) One that changes only ALGO to Normal. Now the geometry stays correct, but again convergence does not seem to improve, see OUTCAR_2.png

Would you mind sharing the INCAR used for the calculation and the relaxed CONTCAR you obtained? I want to check if I have the same settings, and when convergence is reached using these settings whether an electron polaron is calculated or a delocalised electron.

Kind regards,

Jeffrey

[henrique_miranda]

Yes, changing ENCUT, KPOINTS and NCORE is just to speed to up the calculations and get an idea of where the problem might be more quickly without wasting too much computational resources.
Once you know the problem and fix it, you should increase ENCUT and KPOINTS again for production calculations.

So, I changed your INCAR files a bit and ran the calculations locally.
I got the ionic relaxation of your calculation with reduced ENCUT and gamma-only to converge.
Turns out I made a mistake in advising you to change ALGO from Fast to Normal since it is very unlikely that that was the problem.
Sorry for this.

Attached you find a set of input files (step2_ibrion1.zip) that I ran locally and converged below the threshold in the forces that you specified.
The POSCAR however, is not the one that you have sent but an intermediate step that I obtained by using IBRION=2.

I think this is the crucial point: if you are close to convergence IBRION=1 works well but not if you are far:
https://www.vasp.at/wiki/index.php/IBRI ... (RMM-DIIS).

I think you can fix this problem you reported by either:
1. run the calculation with IBRION=2 (ibrion2.zip): the forces indeed look more noisy during the minimization than IBRION=1, but you should monitor the total energy as well, if the total energy is not changing then the ionic minimization is not happening. I think this is what is happening in your calculations where the max force remains constant with IBRION=1.
2. run the calculations with IBRION=1 but first run a few ionic minimization steps with IBRION=2. You might also try setting NFREE=10 and lowering POTIM, but in my testing it did not help much.

These input files I am sharing are just meant to demonstrate the behavior of the ionic minimizer.
For production calculations, you should revert to the INCAR settings you shared initially.

Hope this helps!

[jeffrey_zom]

High Miranda,

I think you are right in that my attention had be to much focused on the forces, as I use a force criterion, rather than the convergence of the energy.
As I showed in my very first post, I used IBRION=2, and based on the forces, it seemed to perform much worse, resulting in erraticly changing forces that were much higher.
Interestingly enough, IBRION=2 in the end did seem to be solution.

For completeness of the discussion and in hopes of helping others in the future facing similar convergence problems,
I have provided some figures in the attachment comparing the different calculations.

OUTCAR_1.png
1) IBRION=2, initially increases forces and generates much more erratic forces, but is eventually able to converge.
3) IBRION=1, does not appear to converge.
5) IBRION=1, NFREE=10, initially shows the same behaviour as without the use of NFREE, but is also adventually able to converge.

OSZICAR_1.png
This show the difference in energy of the newly obtained geometry after an ionic step as function of the ionic steps.
Here too I thought that IBRION=2 was performing much worse, as ΔE was initially going up, not down, while with IBRION=1, ΔE was much lower initially, but also failed to converge.

OSZICAR_2.png
This show the total energy as function of the ionic steps.
This way of looking at the convergence seems to show best how IBRION=2 is the proper choice.


It is somewhat scary to see how both IBRION=2 and IBRION=1, NFREE=10 reach the force criterion adventually, and also converge in terms of ΔE, but result in a significantly different energy (~20 meV)
I should also mention that in my case, the provided POSCAR was very close to the converged solution, as the maximum forces were around 0.02 eV/A, and only for some atoms out of 80. This leads to believe IBRION=1 should be the proper choice. NELECT should remain in the INCAR however, as the geometry of the converged solution without the addition of an additional electron differs significantly, as the additional electron becomes localised and creates a local distortion. Removing NELECT results in very high initial forces.

Thank you so much for you help,

Kind regards,

Jeffrey

[henrique_miranda]

Hi Jeffrey,

I am glad I could help in the end.
Thank you for sharing these figures, they are very instructive for other users that might stumble upon the same problem.
They show that IBRION=2 is clearly better in this case since it is able to quickly decrease the total energy, while IBRION=1 at some point seems to not be updating the positions anymore or doing so very very slowly.

Regarding the structure you find with IBRION=1 and NFREE=10 is this some meta-stable configuration?
It is very odd that the ionic minimization stops there.
It probably warrants some investigation on our side.

You should be able to run these calculations when setting NELECT.
My suggestion to remove it was just to try to figure out if that was causing problems.

You are right that changing ENCUT, KPOINTS and NELECT will mean that one gets large initial forces of course.

To get accurate polaron binding energies, you probably need to make a convergence study with the size of the supercell.
There is where the electrostatic corrections might help by converging the energy differences with smaller supercell sizes but if you include them or not should have no influence in whether you are able to relax the structure.