Big force difference in nearly identical calculations

[steen_lysgaard1]

Hi,

We have an issue with two calculations giving markedly different forces, even though only NSW has been changed in the INCAR. Two folders are attached with input and OUTCARs.

First, we did a singlepoint calculation with NSW=0. However, when the calculated forces were large we started another calculation with NSW=999. To our surprise, it ended after one step since the calculated forces were very small. What can be the cause of this?

The INCARs in question: (The rest of the input and OUTCARs are attached)

Code: Select all

 ENCUT = 520.000000
 SIGMA = 0.050000
 NUPDOWN = 18.000000
 EDIFF = 1.00e-06
 EDIFFG = -3.00e-02
 SYMPREC = 1.00e-04
 ALGO = all
 GGA = PE
 PREC = accurate
 IBRION = 2
 ISIF = 3
 ISMEAR = 0
 ISPIN = 2
 ISTART = 1
 ISYM = 0
 LMAXMIX = 4
 LORBIT = 12
 NELM = 400
 NELMIN = 8
 NSW = 0
 NCORE = 4
 LCHARG = .TRUE.
 LDAU = .TRUE.
 LWAVE = .FALSE.
 LREAL = auto
 LDAUL = 2 -1 -1 -1
 LDAUU = 5.300 0.000 0.000 0.000
 LDAUJ = 0.000 0.000 0.000 0.000
 MAGMOM = 2*5.0000 2*4.0000 22*0.0000

Thanks in advance!
Steen

[andreas.singraber]

Hello Steen!

I had a closer look at your results and could reproduce the behavior you described. My calculations also initially showed large forces for the singlepoint calculation. However, I repeated this run with different random seeds (INCAR tag RANDOM_SEED) and sometimes also the singlepoint calculation would end up in the same minimum as found in the "relax" run. So, I believe that there is no actual relation between the NSW tag changes and the contradicting results. Instead, I think that there are many close-by local minima and therefore the optimizer sometimes does not converge to the global minimum, depending on the starting conditions and maybe also rounding errors along the minimization (e.g. MPI-based reduction algorithms do not guarantee an order of summation). Actually, magnetic calculations are known to be notoriously difficult to converge and your system seems to fit into this category. Therefore, I suggest to use the procedure described here to get reliable results. I followed this recipe five times with different random seeds and could always obtain the desired energy minimum with the vanishing forces. Maybe it even suffices to just reduce the step size for the conjugate gradient algorithm by reducing the value for the TIME tag from 0.4 (default) to 0.05, as described in the link above. However, I have not tested whether this change on its own improves the reliability.

All the best,
Andreas Singraber

[andreas.singraber]

Hello!

I just had a discussion with my colleagues about your problem and I got a few more hints about what you can try to improve the reliability of results. Also, I need to correct myself, my assumption that systems like this are hard to converge is most likely wrong, there seem to be other issues with the setup causing the issues. I may also be wrong with the RANDOM_SEED remarks, as I am not sure if this is actually used to initialize the orbitals.

My colleagues gave me the following list of INCAR changes you can try:

• TIME = 0.05
• LASPH = .TRUE., this is recommended for LDA+U
• EDIFF = 1.0E-7 or even 1.0E-8
• ALGO = Damped, as an alternative to ALGO = All, the tag TIME is also important here.

I am currently running a few tests myself, I hope this helps...

Best,
Andreas Singraber

[steen_lysgaard1]

Hi Andreas,

Many thanks for your suggestions. We have implemented them all now (except ALGO=Damped, which couldn't converge), and the stability is a little better. However it still converges to a high force approximately one out of six times.

The issue, as you correctly point out, is that the density around the Fe atoms is a little different in the low and high force cases (different orbitals are occupied), leading to no energy difference but a large force difference. This is quite obvious from charge density difference plots.

We were interested in whether we could use the high force result in an ionic relaxation to arrive at a lower ground state energy / global minimum. However, every time we set NSW>0 we get the low force result.

Thus from now on we simply set NSW=1 (i.e. almost a single point calculation, but we are only interested in the energy and forces anyway)

Best regards,
Steen