Spin Orbit Coupling in lm and lmf (7.9)

Programs lm and lmf can include spin-orbit coupling to the scalar Dirac hamiltonian. At present lm has more functionality, e.g. it is a fully noncollinear code, with the ability to apply noncollinear external fields and also has the option to use the fully relativistic Dirac equations, where spin orbit is not included perturbatively, but lm is limited to the Atomic Spheres approximation. At present, lmf has the ability to handle SO coupling by adding λL·S to the scalar Dirac Hamiltonian. (In what follows, λ is implicit in wherever L·S is mentioned.) At present the density is restricted to spin-diagonal parts.

Turn on spin-orbit (SO) coupling with HAM_SO where SO is one of:

  1. SO=0 : No spin orbit coupling is included.
  2. SO=1 : Usual L·S coupling. The full L·S is added to the scalar Dirac hamiltonian. Only the spin-diagonal parts of the density are retained at present
  3. SO=2 : Lz·Sz part of L·S. The hamiltonian is still spin diagonal so the spin channels remain distinct.
  4. SO=3 : Lz·Sz + (L·SLz·Sz), with the last term included perturbatively in the eigenvalues only, as described below. This generates eigenvalues very similar to L·S for a given potential, but the eigenfunctions are generated from H0+Lz·Sz only. As a result the eigenfunctions remain spin diagonal. There is some effect on the density, but the approximation seems to be rather good. The error in the eigenfunctions is second order in L·S. This mode has advantages, particularly when used in conjunction with GW.
This Figure shows the valence bands of Fe near the Fermi level, calculated with SO coupling and without. Bands with SO coupling are shown as purple (majority spin) and red (minority spin). (Note the two spin channels are coupled, so neither has a sharply defined spin. Nevertheless except at crossing points the spins are nearly distinct. The projection onto spins is accomplished using color weights.)

Compare these bands to green and yellow dashed lines (majority and minority spins with SO=0). When SO=0 spins are strictly decoupled in both LDA and GW.

The SO=0 and SO=1 cases are very similar (Fe is a light element), but splittings can be seen. Note especially the splitting of the three states just above EF at H (EF=0 in the Figure), and the splitting of the t2g states at Γ near −2.3 eV.

The figure on the right enlarges the region EF to better show the splittings, especially at H.

The next figure compares SO=2 and SO=1 cases (the potential is kept fixed).
You can see that the splitting at H is well captured by the Lz·Sz part of L·S.

There is a small splitting at Γ near −0.2 eV which Lz·Sz doesn't get quite right. Also you can see that there where two bands of the same spin cross in the Lz·Sz case, they get split with the full L·S; see for example midway between H and N, near −0.2 eV.

In the next figure SO=3 is comared to SO=1.

The dashed lines almost perfectly overlay the solid ones, showing that the SO=3 approximation is nearly perfect, except where two degenerate levels of the same spin come together. Compare the region near −1 eV between Γ and H, or near −0.9 eV between Γ and P. In the latter case, one crossing is not well described but the agreement is subtantially improved relative to SO=2.

SO=3 also moves a little closer to the full SO=1 calculation of the magnetic moment and charge density. In the Table below, the density was made self-consistent with SO=1. In the RMS DQ line, a one-shot calculation was performed with SO=0,1,2,3, and the data shows the deviation in density relative to the self-consistent SO=1 density after the first iteration. The fully self-consistent magnetic moments are also shown for each approximation. Evidently the Lz·Sz density and moment are further removed from the SO=1 result than with no SO coupling at all, while the SO=3 case is pretty close to the full result.

RMS change in the density starting from a self-consistent density without spin orbit coupling
SO             0        1         2         3
RMS DQ         1.90e-5  0         3.13e-5   1.38e-5
M              2.24765  2.24861   2.24689   2.24814
Only a few checks have been made, but preliminary estimates of SO contribution to total energy seem to well described by SO=3. For example, in the top level directory run the Fe test case :
  fp/test/test.fp felz 4
This test gives the following total energies. (Note: there seems to be a term missing in the total energy when SO=3. This hasn't been worked out yet.)
 case             ehf      orb moment
 no SO       -2541.03436
 LzSz        -2541.03450    -0.0492
 L.S         -2541.03478    -0.0513
 L.S (pert)  -2541.03478    -0.0492

Perturbative inclusion of L+·S

Because the Fe potential is magnetic the the L+·S term is comparatively less important. It occurs as second order in perturbation theory whereas Lz·Sz is first order. But when the potential is nonmagnetic the L+·S term is comparable to Lz·Sz because ε+ε vanishes for the same bands. Thus the standard second order expression in perturbation theory, |λL+·S|2/(ε+ε), is not meaningful and should become ±λL+·S, which is is of the same order as the spin diagonal matrix element.

Trying to avoid pitfalls when |λL+·S| is large compared to |ε+−ε| is a bit subtle. Exact diagonalization has pitfalls because there is no general way to unambiguously partition the eigenvalues into spin-up and spin-down parts, in other words, retain an unambiguous one-to-one correspondence between the eigenvalues and eigenfunctions.

We solve the problem as follows. First H0+Lz·Sz is diagonalized exactly, and the L+·S parts are rotated to the basis of H0+Lz·Sz eigenfunctions. The residual L+·S still has only off-diagonal blocks in spin space. Call the 12 block V with matrix element Vij, the ith spin-up eigenvalue ε+i and the jth spin-down eigenvalue εj.

In second order perturbation theory, matrix element Vij contributes a shift δε+i=|Vij|2/(ε+i−εj) to ε+i. The difficulty, as noted, is that the denominator can vanish or be small. But the standard technique to resolve this case, exact diagonalization, is problematic because it is no longer possible to assign a particular eigenvalue to a spin channel.

Instead we choose to do exact diagonalization on a subset of the matrix V: we diagonalize exactly the coupling Vij between the two elements ε+i and εj only. Thus, coupling involving two states alone is computed to all orders, while third order and higher order couplings that connect three inequivalent states is omitted. This works very well in practice, but omits splittings when three different states all become degenerate. You can see instances of this in the Fe figures, above. The effects are minor, but readily seen. It would be possible to include such couplings to still higher order, but this has not been done.

When bands are not spin polarized: a case where SO coupling is large

To include SO coupling in the nonmagnetic case, you must make it magnetic by setting HAM_NSPIN=2. You must also set HAM_SO. The code will generate identical bands for spin-up and spin-down channels. You should not use SO=2 in such a case, even though it worked pretty well for Fe as we saw above.

To show how the different approximations work, consider Au -- a heavy element with large SO coupling. The figures on the left and right compare the following approximations:

SO          left    right
0 (none)    blue    blue
1 (L·S)     red     red
2 (Lz·Sz)    -      green 
3 (pert)    green   -
Circles are photoemission data.

Blue and red are very different: L·S now generates shifts on the order of 1 eV.

Comparing green and red in the right figure it is evident, Lz·Sz causes some splitting but it is a poor approximation to L·S.

The perturbative approach, on the other hand, is rather good, even for this heavy element (compare green and red in the left figure). The difference is negligible on the scale of LDA errors (compare to photoemission data). With (L·S) so large, the charge density is modified. As the following table shows the perturbative L·S approach captures only about 1/3 of the density shift: essentially the same shift as Lz·Sz gets.

RMS change in the density starting from a self-consistent density without spin orbit coupling
SO             0        1         2         3
RMS DQ         3.74e-4  0         2.58e-4   2.56e-4
Fortunately the perturbation in charge density is still small: the bands look essentially similar whether the potential is generated from the converged SO=1 density, or made self-consistent. All three approximations SO=1,2,3 work with lmf, for LDA or LDA+U calculations. The QSGW code at present requires the spin channels be kept separate and works with SO=2,3 only.