very different total energies with ISYM=-1 and ISYM=0?

[sophie_weber]

Hello,

I am trying to do a series of constrained magnetic calculations where I cant the magnetic moments of Ni in a vacuum-terminated slab of rock-salt NiO very slowly from their bulk, collinear in-plane directions in the in-plane perpendicular direction. The specific pattern of canting I'm performing should actually preserve the symmetry of the non-canted structure, but anyway since these calcuations are noncollinear and my understanding is one should always turn off symmetries in this case, I initially set ISYM=0.

I noticed afterwards however, that if I instead turn all symmetries, including time-reversal, the total energy changes by order of an meV, which for my purposes is a very big difference. This makes absolutely no sense to me. My understanding is that ISYM=0 assumes time-reversal as a symmetry, and if this is the case according to the POSCAR and MAGMOM, it reduces the k-mesh accordingly. My slab I'm pretty sure should have time-reversal plus a translation as a symmetry, so by default I'd assume that the k-mesh of ISYM=0 would be half that of the ISYM=-1 calculation and then I could maybe understand that you could end up with different energies. But checking both OUTCARS, I see that the number of kpoints used for each case is identical, so ISYM=0 did not change the k-mesh used. Moreover, with these constrained magnetic calculations I could imagine maybe converging to slightly different magnetic states, but if you look at the OUTCAR and OSZICARs for both calculations they are identical basically to the decimal place. Also, the penalty energies for both cases are order 10^-6 eV, so this should not be affecting anything (I do notice the penalty energy for ISYM=-1 is smaller, but still I'd be surprised if this matters since they are both so much smaller than the energy difference between the two calculations).

So the fact that I get such different energies for completely identical calculations, as far as I can tell, with the only difference being ISYM=0 or ISYM=-1, is extremely concerning. I'd really appreciate it if first, someone could explain how this can happen when the k-mesh and the final magnetic states are identical, and second, which setting is the more "correct" one to use. Interestingly, in the case where I cant out of plane, in which case I definitely am breaking the symmetry compared to the collinear state, the ISYM=-1 results make more "sense", whereas in the case where I'm canting in-plane as I describe above, the overall trend as a function of canting angle for the total energy makes more sense with ISYM=0. But this may just be a coincidence.

I attach the input files for one calculation (they are identical for both calculations, except ISYM is switched as described above) and the OUTCAR and OSZICARs for both ISYM settings. Any insight into this would be greatly appreciated. Thank you in advance!

[fabien_tran1]

Hi,

Noncollinear calculations are particularly prone to various problems of this kind, and this is sometimes not easy to solve the problem or even to know what is the reason. For the moment I have two questions.

1) The calculation ISYM=0 seems to have been restarted from a previously converged calculation, right? Would you get the same results if you start it from scratch, i.e. without preexisting WAVECAR? In other words, is the final total energy always the same, independently of how you start the calculation? I have myself experienced a strong dependency of the final total energy on the way the calculation is started, but this was mostly with MGGA and hybrid functionals.

2) What do you mean with "they are identical basically to the decimal place"?

[sophie_weber]

Hi Fabien,

Thanks for your response.
1) This is correct, it was restarted from a previously converged calculation. I haven't checked every single case, but for the few that I've checked, it doesn't seem to matter whether the calculations start from preconverged calculations or from scratch, the energies for ISYM=-1 and ISYM=0 are quite consistent (but different from each other). I am just using plain LDA+U

2) I mean that if you look at the converged magnetic moments of the atoms in the POSCAR (either in the OSZICAR or OUTCAR), they are identical both qualitatively and numerically between ISYM=0 and ISYM=-1, except for a few cases the magnitude of individual moments differs by 0.001 uB. It would shock me if this would lead to a ~meV energy difference.

Let me know if anything else is unclear!

[sophie_weber]

Shoot, I just noted however that for one calculation (a different compound, but same procedure of canting top and bottom spins on a vacuum-terminated slab), for the same setting ISYM=-1, changing ntasks-per-node and NCORE (halfing each, but keeping the same number of total cores) also results in substantially different energies. Have you run into such a dependence on parallelization settings?! In any case, it does seem like I need to be very carefully to use the exact same settings in each case to compare properly. I'll look into this more, but I'm still fairly certain that for the examples I gave you, the energy did not change based on reading in from another starting file.

In cases like this, how on earth does one choose the "correct" energies?

[fabien_tran1]

This is of course a problem.

For the problem ISYM=0/ISYM=-1 could you please repeat these two calculations with the following settings:
-ENCUT=300 to make the calculations faster.
-EDIFF=1E-4
-Reduce the k-mesh to make the calculations faster.
-Start them from scratch (no WAVECAR or charge files present).
-Without onsite U, i.e. just PBE.
-Same settings for the parallelization for both calculations.

Besides, do you think that it would be possible to construct a unit cell that is much smaller and should lead to the same problem?

[sophie_weber]

Sure, I submitted them now, they should be done by tomorrow. I'll let you know when I have results.

Best
Sophie

[sophie_weber]

Hello,

Attached are the input files, and the OUTCAR and OSZICAR for both cases, with the inputs you specified. As you can see the same thing happens, i.e. same k-mesh and everything else is used but the total energies differ by order meV.

In principle, in terms of the physics inferred it is not an issue if the absolute energies differ as long as trends are the same for both settings (though in my experience this is sometimes not the case, which make the problem worse). The main thing that bothers me in terms of reproducibility of the code is that whatever VASP is doing to print out substantially different energies, when everything else seems identical, is not rigorously understood by anyone, including the developers. I understand that noncollinear calculations are incredibly difficult to converge, but at the point that the outputs of the two calculations and the settings/final kmesh used appear identical but the ISYM flag leads to different total energies It feels difficult to trust any of the results completely.

Thanks a lot for your help, any further insight would be greatly appreciated!

[fabien_tran1]

Hi,

Indeed, these two calculations ISYM=-1 and 0 should in principle provide exactly the same results, in particular since you set LSORBIT=.TRUE., which means that there is no time-reversal symmetry (not implemented in the case of SOC) and therefore ISYM=-1 and 0 are the same. The sets of k-points in both calculations are exactly the same, the only difference being that they are ordered differently (see OUTCAR). This difference may however play a role because of the way the initial orbitals are generated (with a random generator, I don't know the details), leading to slightly different starting orbitals. Then, the hypothesis is that these two different sets of starting orbitals lead to different solutions, because of the possible flatness of the electronic energy surface (and/or the inability of the ALGO=All method to work properly).

A question about the calculations from your first post: how did you start the calculation OUT_isym_zero? Did you start it from the WAVECAR file of the calculation OUT_isym_minusone?

[sophie_weber]

Hello,

Thanks for the response. That is indeed concerning that they should in principle produce the same energies.

I'm sorry I don't remember exactly the settings for the initial files I sent you; I should have saved those directories. For sure I did not read in the wavecar of the exact ISYM=-1 calculation I sent you into the ISYM=0 calculation; I believe I did read in an ISYM=-1 WAVECAR however.

If it's important I can repeat the calculations, let me know.
You also the ALGO=All possible being an issue; should I repeat the reduced-size calculations with the default algorithm? The reason I switched was that the default minimization algorithn (Davidson I believe?) was taking many hundreds of iterations to converge, and ALGO=All was much quicker.