# Spin-orbit coupling (SOC)

## Objectives

- Setting up calculations with SOC.
- Understanding the different modes of SOC.
- Calculating the magnetocrystalline anisotropy energy (MAE)

## Introduction

In this tutorial we will perform calculations related to physics in which spin-orbit coupling (SOC) plays an important role. In calculations without explicit treatment of noncollinear magnetism Fleur treats SOC in so-called second variation. This means that a conventional Hamiltonian matrix, i.e., a Hamiltonian matrix with a scalar-relativistic description in the MT spheres, is set up and diagonalized. After diagonalization the eigenfunctions are used as a new basis for which a new Hamiltonian matrix with SOC is set up and diagonalized again to obtain the final energy eigenvalues and eigenfunctions.

This calculation procedure implies that the number of eigenfunctions calculated in the first step now
becomes a cutoff parameter defining the basis set size for the second step. This parameter can be set in
`/calculationSetup/cutoffs/@numbands`

. A value of zero indicates that a material-dependent default value is chosen
by the program. In principle SOC calculations have to be converged with respect to the numbands parameter. Due to time
limitations in this tutorial we will not do this, but you are free to perform convergence tests to get a feeling for the relevance of
this parameter. Since some of the investigated quantities are sensitive to very small energy differences a significant
quantitative dependence on the cutoff parameters is expected.

## Magnetocrystalline anisotropy of hcp Co

In perfect infinite ferromagnetic crystals the preferred magnetization direction is determined by spin-orbit coupling. In this example we will calculate the energy difference between different magnetization directions. Note that this typically is a very small quantity. Also be aware that magnetocrystalline anisotropy is not the only source of magnetic anisotropy. For example, in finite systems the shape of the sample also contributes to the magnetic anisotropy.

As example system we choose hcp Co. An inpgen input `inpCo.txt`

is available in the `Co`

subdirectory. Change into that directory and inspect the file.

```
cd Co ; cat inpCo.txt
```

Note the line "&soc 0.10 0.37 /". With this line we define a spin quantization axis (SQA) in terms of the two angles and (spherical coordinate system). These are not the angles we will use in the actual calculations. We choose these angles in the inpgen input because the SQA breaks the symmetry of the crystal. Here we want to have the freedom to choose any other SQA after generating the inp.xml file. So we need to define an axis that breaks as many symmetries as possible. Otherwise a later choice of another axis in the inp.xml file may be incompatible to the here detected symmetry operations.

The line "&kpt div1=9 div2=9 div3=5 /" defines a k-point mesh directly in the inpgen input. As no name is provided for this mesh and it is the only k-point mesh specified in this file it overrides the default.

Generate the Fleur input.

```
inpgen -f inpCo.txt
```

In the `inp.xml`

file the angles appear in `/calculationSetup/soc/@theta`

and `/calculationSetup/soc/@phi`

. The SOC
calculation is activated with the XML attribute `calculationSetup/soc/@l_soc`

. If there is a spin quantization axis
defined in the inpgen input spin-orbit coupling is activated by default.

To investigate different SQAs prepare separate directories for each SQA and copy all xml files to each of these directories.

- theta="0.0", phi="0.0" (z axis)
- theta="Pi/4.0", phi="0.0"
- theta="Pi/2.0", phi="0.0" (x axis)
- theta="Pi/2.0", phi="Pi/2.0" (y axis)

```
for d in th-0.00-ph-0.00 th-0.25-ph-0.00 th-0.50-ph-0.00 th-0.50-ph-0.50
do
mkdir $d
cp *.xml $d
done
```

Now we have to adopt in each of the directories the `inp.xml`

file to the respective SQA, e.g., in the `th-0.50-ph-0.50`

directory set

```
<soc l_soc="T" theta="0.50*Pi" phi="0.50*Pi" spav="F"/>
```

The following command will do the modification, you might want to check it by inspecting the files.

```
sed -e 's%<soc.*>%<soc l_soc="T" theta="0.00*Pi" phi="0.00*Pi" spav="F"/>%' inp.xml >th-0.00-ph-0.00/inp.xml
sed -e 's%<soc.*>%<soc l_soc="T" theta="0.25*Pi" phi="0.50*Pi" spav="F"/>%' inp.xml >th-0.25-ph-0.00/inp.xml
sed -e 's%<soc.*>%<soc l_soc="T" theta="0.50*Pi" phi="0.00*Pi" spav="F"/>%' inp.xml >th-0.50-ph-0.00/inp.xml
sed -e 's%<soc.*>%<soc l_soc="T" theta="0.50*Pi" phi="0.50*Pi" spav="F"/>%' inp.xml >th-0.50-ph-0.50/inp.xml
```

Now we are ready to run FLEUR in these directories

```
for d in th-0.00-ph-0.00 th-0.25-ph-0.00 th-0.50-ph-0.00 th-0.50-ph-0.50
do
cd $d
fleur_MPI
cd -
done
```

Extract the converged total energy for each of the calculations. hcp Co is an uniaxial system, i.e., the magnetocrystalline anisotropy in the x-y plane is negligible. Verify that the total energy difference between calculations 3 and 4 is much smaller than the difference to calculation 1. What is the "easy axis" (energetically most preferable), what is the "hard axis"?

```
for f in $(find . -name "out"); do grep -H "total energy" $f | tail -1; done
```

Depending on the unit cell there are different models to describe the magnetocrystalline anisotropy. For uniaxial systems the total energy in dependence of the magnetization direction is expressed by the equation

In approximation we assume this expansion in terms of powers of sines of to end after the term. Use the total
energies from calculations 1, 2, and 3 to set up a system of linear equations and determine the coefficients ,
, and . In the `calcCoefficients.ipynb`

you will find a Python script that you can modify to do this.

What are the values of the coefficients you calculated?

We want to visualize the anisotropy. For this we ignore the term. Plot it. **For this adapt the gnuplot script in the file plotEnergySurface.plt by plugging in and **. It will generate a file

`energySurface.png`

with the visualization.```
cd ..; gnuplot < plotEnergySurface.plt
```

The result of plotting the magnetic anisotropy is a file `surface.png`

that should be similar to the plot below (not identical due to slightly different parameters). The distance between the origin and a point on the surface is the energy advantage the easy axis has over the axis in the direction of the point on the surface.

Note that the procedure of determining (extended) Heisenberg model parameters is similar to what is sketched in this example for the magnetocrystalline anisotropy: One performs DFT calculations for different magnetic configurations and plugs the results of these calculations in equation systems where each equation relates to the respective configuration in the Heisenberg model.

## Rashba effect in BiTeCl

The Rashba effect, similarly to the Dresselhaus effect, is a band splitting due to spin-orbit coupling that appears in systems with broken inversion symmetry. While the Rashba effect appears in low dimensional systems (films, heterostructures) and uniaxial bulk systems, the Dresselhaus effect is observed in cubic bulk systems. In this exercise we calculate the Rashba splitting in bulk BiTeCl. Note that it also appears in other bismuth tellurohalides (BiTeX with X=Cl,Br,I).

BiTeCl features a hexagonal unit cell (use latsys='hP' for inpgen) with 6 atoms and the lattice parameters
Angstrom and Angstrom. The `inpBiTeCl.txt`

inpgen input in the `BiTeCl`

directory represents such a unit cell.
Change into that directory and inspect the file.

```
cd BiTeCl
cat inpBiTeCl.txt
```

We will calculate this system without spin polarization. In such a case the choice of the SQA is not important, but some SQA should be specified in the inpgen input. We choose one that does not break too many symmetries. Also note that we specified two different k-point sets in the inpgen input. Both feature relatively few points to reduce the computational demands. The second k-point set is a user-defined k-point path that deviates from the default path for this unit cell.

Generate the Fleur input.

```
inpgen -f inpBiTeCl.txt
```

Create the directories `withSOC`

and `withoutSOC`

and copy all xml files to these directories.

```
mkdir withoutSOC
mkdir withSOC
```

```
cp *.xml withoutSOC/
cp *.xml withSOC
```

In the `withoutSOC`

directory we now set `/calculationSetup/soc/l_soc`

to `"F"`

to disable SOC in the respective calculation.

```
sed -e 's%l_soc="T"%l_soc="F"%' inp.xml >withoutSOC/inp.xml
```

Now we will perform the self consistency using aiida-fleur in both directories. After that we also perform a bandstructure calculation.

```
These calculations will take some time :-).
```

```
cd withoutSOC ; aiida-fleur launch scf ; aiida-fleur launch band -P scf.json ;cd -
cd withSOC ; aiida-fleur launch scf ; aiida-fleur launch band -P scf.json ;cd -
```

If the calculations are finished, we can plot these bandstructures:

```
display < withoutSOC/bandstructure.png
display < withSOC/bandstructure.png
```

The results look similar to the plots below. However, due to their larger energy range the details close to the point are harder to see.

Compare the two plots. You should observe differences at the gamma point in the bands defining the band gap. For the case with SOC there actually is a spin-dependent band splitting where the states for one spin can be found at and those for the other spin at . We don't see the spin in this plot because we made a non-spinpolarized calculation. Otherwise this can also be highlighted with the weights that are stored in the `banddos.hdf`

file.

## Learn more

- https://www.flapw.de/MaX-7.0/documentation/socSecVar/

```
```