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 $\theta=0.10$ and $\varphi=0.37$ (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 $9 \times 9 \times 5$ 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

$\begin{equation} E_\text{total} = K_0 + K_1 \text{sin}^2(\theta) + K_2 \text{sin}^4(\theta) + \dots \end{equation}$

In approximation we assume this expansion in terms of powers of sines of $\theta$ to end after the $K_2$ term. Use the total energies from calculations 1, 2, and 3 to set up a system of linear equations and determine the coefficients $K_0$ , $K_1$ , and $K_2$ . 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 $K_0$ term. Plot it. For this adapt the gnuplot script in the file plotEnergySurface.plt by plugging in $K_1$ and $K_2$ . 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.

$cap=Visualization of magnetocrystalline anisotropy of hcp Co obtained with the script above. The unit of the axes is Htr per unit cell.,width=0.8\textwidth$

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 $a=4.2426$ Angstrom and $c=12.3970$ 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 $\Gamma$ point are harder to see.

$cap=Band structure of BiTeCl without SOC.,width=0.8\textwidth$

$cap=Band structure of BiTeCl with SOC.,width=0.8\textwidth$

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 $+\vec{k}$ and those for the other spin at $-\vec{k}$ . 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/

Keys	Action
`?`	Open this help
`n`	Next page
`p`	Previous page
`s`	Search