Unphysical/inconsistent result when changing SAXIS

[aedstrom]

Hi,

I attach a calculation including spin-polarization with spin-orbit coupling and constrained noncollinear magnetism, for a monolayer of Neel type antiferromagnet FePS3.

Running a calculation with e.g.
SAXIS = 0 0 1
and then restating from that calculation but rotating the SAXIS gives unphysical results.
In particular, when it is rotated to
SAXIS = 0 0 -1
One finds a different energy than that for SAXIS = 0 0 +1, inconsistent with time-reversal symmetry.

Moreover, if one instead starts the calculation from scratch but with
SAXIS = 0 0 -1
and then rotates it to
SAXIS = 0 0 1
one finds the opposite energies for the two cases, pointing again towards it being an unphysical, numerical inconsistency.

I have reproduced the exact same result with both VASP 5.4.4 and with VASP 6.4.3.

Below you see a plot of the energy as a function of the angle theta between SAXIS and the out-of-plane direction of the monolayer.
This should be 180 degree periodic as it is not.
Starting from the opposite direction of SAXIS, would invert the curve around theta=90.

E_of_th_FePS3.png

The calculations were done with LDAUTYPE = 2 and a U = 4 eV on the Fe d states, as is needed to get an AFM insulator.
Setting U=0, gives a metallic solution, but it also removes the above problem and the energy comes out the same for SAXIS = 0 0 1 and SAXIS = 0 0 -1.

[ahampel]

Hi,

thank you for reaching out to us via the official VASP forum.

I had a quick look at the calculations. Let me first stress that all symmetries are switched off during these runs as noted in the OUTCAR. I then noticed that in the end the magnetization is okay for both runs, but the orbital polarization is completely different for both runs and also the SOC term seems to be very different as the energy contributions are different. So somehow switching SAXIS here seems to change the systems orbital polarization. I will ask later in the team if someone as some comments about it and will get back to you.

Best,
Alex

[ahampel]

Hi,

I discussed a bit with the team and the problem seems to be quite technical an inheritly part of how LDA+U is working and not necessarily a bug in the code. In your system LDA+U is enforcing a certain orbital polarization hence populating a certain orbital with + or - m_l. The SOC term would like ideally to pair + and - ml orbitals to form the coupled SOC orbitals. So they both work a little bit against each other. I can confirm that running the calculation from scratch gives the same energy and the same magnetization. And I noted that running with SAXIS +1 or -1 gives an orbital polarization that switches from SAXIS +1 to -1 of the atom holding the large magnetic moment (Fe d). Thus it seems that the energy minimum with LDAU for a given SAXIS seems to force a specific orbital polarization. For 0 0 1 it is for me:

Code: Select all

orbital moment (z)

# of ion     p       d       tot
----------------------------------------
 1        -0.001   0.843   0.842
 2         0.001  -0.843  -0.842

and for 0 0 -1:

Code: Select all

orbital moment (z)                      
                                        
# of ion     p       d       tot        
----------------------------------------
 1         0.001  -0.843  -0.842        
 2        -0.001   0.843   0.842        

Now I hope that this is not a physical problem, but it is a consequence of this interplay of the SOC term and the LDAU term.

I also learned that the SAXIS term only sets the direction for the SOC term, it does internally not influence the LDAU term. This is also the problem when restarting from a WAVECAR, the interpretation of the orbitals is wrong now for LDAU and it works against the SOC term now differently. I am fairly certain that you get trapped in a certain local minimum now that you can't exit. I think this because the reported SOC matrix for the ion 1 is now suddenly very different with opposite sign compared to starting from scratch. It seems to be not able to change and the total energy is as you also find roughly 90 meV higher. I only run now +1 and -1 SAXIS and nothing in between but it fits to the plot you shared. This is probably not even fixable when going to super small angular changes of SAXIS. However, it should work when you start from scratch, which we or I should try, or is then hard to find to follow that minimum? I don't see why there would be a bug in the code only when reading a WAVECAR but I will leave that option open if this behavior turns out to be not just a local minimum but something really wrong.

Also it was suggested to get the anisotropy curve of the energy here to better not use constrained moments with a penalty function but rather just keep SAXIS the same and slightly rotate the magnetic moment. Or do you then converge back to the minimal solution? I guess there is a reason why you used constrained moments here...

Anyway, I hopefully find some time tomorrow to try these things myself and create an magnetic anisotropy curve for the system that correctly has 180 deg periodicity. But if you like you can maybe already comment on these things, i.e. possibility of starting all of the calculations from scratch and / or removing constrains and using MAGMOM instead of changing SAXIS.

Best regards,
Alex