Spin-Orbit coupling¶

Author: Ramón Cuadrado (University of Southampton) & Nils Wittemeier (ICN2)

The main goal of this tutorial is to cover the basic notions to perform calculations using the spin and spin-orbit coupling (SOC) in Siesta.

Note

When performing a SOC calculation it is crucial to:

carefully converge the k–points sampling
carefully converge the mesh cutoff
converge the density matrix with a tolerance below 10^-5
use relativistic pseudopotentials (PPs) and use the non–linear core corrections when the PPs are built.

All the information regarding the basis set and all the parameters used in a common SIESTA calculation are valid for the SO.

In the next two examples the above mentioned parameters (mesh cutoff, etc) are not converged because we want to speed up the calculations to show an initial (and dirty) results, but remember, for a real calculation those values have to be converged.

An additional advice: remember to read the SOC section in the main Siesta manual!

Performing calculations with different spin configuration.¶

Hint

Please move to directory Example-1.

Note

The original fdf for this example can be found in example_1.fdf.original. We have made some changes (can you see them?) to speed the calculation while maintaining reasonable quality results.

In this example you will calculate the band structure, density-of-states (DOS), and Mulliken charges of bulk FePt using different spin options.

Getting started¶

To get started inspect the input file example_1.fdf file, and make sure that you are familiar with all options in the input file. Refer to the user manual to clarify any doubts. Alongside the input file you will find two pseudopotential files one for Fe and one for Pt (Fe_fept_SOC.psf, Pt_fept_SOC.psf)

In addition, to the usual flags we have included the following options:

to calculate and store and band structure information

BandLinesScale          ReciprocalLatticeVectors
%block BandLines
0.0000 0.0000 0.0000 Gamma
0.5000 0.0000 0.0000 X
0.5000 0.5000 0.0000 M
0.0000 0.0000 0.0000 Gamma
0.0000 0.0000 0.5000 Z
0.0000 0.5000 0.5000 R
0.5000 0.5000 0.5000 A
0.0000 0.0000 0.5000 Z
%endblock BandLines

to calculate (projected) density-of-states

%block PDOS.kgrid_Monkhorst_Pack
  30    0    0    0.0
   0   55    0    0.0
   0    0   55    0.0
%endblock PDOS.kgrid_Monkhorst_Pack
%block ProjectedDensityOfStates
    -10.00  5.00  0.02  1000  eV
%endblock ProjectedDensityOfStates

to calculate and print mulliken charges
```
WriteMullikenPop 1
```

Simulations¶

Run one SIESTA calculations for each possible options for the Spin flag: none, polarized, non-colinear, spin-orbit. Please, change the SystemLabel for each calculation to avoid overwriting output files or perform each calculation in a separate subdirectory.

Analysis¶

Visualize the bands structure, and compare the results. Do any of the calculations give the same band structure? Why?

Templates for plotting the bands/DOS with gnuplot can be found in the .gnp files.
Note

The default version of gnuplot on MareNostrum is quite old (v4.6). If you plan on running gnuplot on MareNostrum, we recommend loading a recent version first:
```
module load gnuplot/5.2.3
```
Note

If you prefer performing the analysis with sisl in a jupyter notebook and you have both installed, please feel free to give it go. To get started please, have a look at example_1.ipynb. Unfortunately, we can not support running the jupyter notebooks on MareNostrum.
Visualize the density-of-states (.DOS file).
Note
- For Spin none the .DOS file contains two columns: energy and DOS.
- For Spin polarized the .DOS file contains three columns: energy, DOS↑, DOS↓ (DOS for spin-up/spin-down electrons). The integral of DOS↑ (DOS↓) up to E_F will give the number of electrons with spin-up (spin-down). The total density of states is the sum of DOS↑ and DOS↓.
- For Spin non-colinear or Spin spin-orbit the .DOS file contains five columns: energy, DOS↑, DOS↓, Re{DOS↑↓} and Im{DOS↑↓}.
The integral of Re{DOS↑↓} up to E_F will give the magnetization along X axis, Mx, and the integral of Im DOS↑↓ up to E_F , My. The total density of states is the sum of DOS↑ and DOS↓.
Inspect the Mulliken charges printed in the output file (search for ‘Mulliken’). Are all atoms/orbitals equally spin-polarized? If not, electrons in which orbital do primarily contribute the magnetic properties?
More advanced: Visualize the projected density-of-states (pDOS)

Hint

In the interest of time, you might want to go through the examples below first before doing this.

Note

The .PDOS file can be processed as discussed in this how-to.

You can plot the PDOS↑, PDOS↓, PDOS↑↓, etc. for each atomic species or the total DOS (the sum of the relevant columns). It is worth noting that in general, for the transition metals the most relevant orbitals will be the d ones. You can play with the DOS/pDOS analysis programs to choose other PDOS curves.

Calculation of the magnetic anisotropy of Pt₂.¶

Hint

Please move to directory Example-2

Note

The original fdf for this example can be found in example_2z.fdf.original. We have made some changes (can you see them?) to speed the calculation while maintaining reasonable quality results.

In this example, you will calculate the magnetic anisotropy of a Pt₂ dimer. This means you will have to run multiple SIESTA calculations each with magnetic moments aligned along different directions and extract the total energy.

To get started inspect the input file example_2z.fdf file. Use this input file and the pseudopotential of Pt atom, Pt_pt2_SOC.psf (note that the pseudo for Pt is different from the one used in the first example), and calculate the total self-consistent energy Etot for three highly symmetric orientations of the initial magnetization, namely, X, Y and Z, keeping fixed the physical dimer (the (x,y,z) coordinates of each atom).

Each calculation has to be initialized modifying the DM.InitSpin, updating the orientation of each atomic magnetization. As a result one has to have three different energy values for each (θ,ϕ) pair.

Remember to rename each one of the .fdf files, change the SystemLabel flag inside for each magnetic orientation and do not re-use the density matrix (.DM file) obtained in other orientations. Reading the density matrix will cause SIESTA to ignore the DM.InitSpin. To calculate the total energy for each specific magnetic configuration you need to start from scratch.

Note

You can also try to change the orientation of the dimer (say put it along the Z or Y axes) and compute again the energy values associated to different spin orientations and compare the results. Some of the total energies (which?) have to be the same.

Note

What would happen if you used Spin non-colinear instead of Spin spin-orbit? Would you obtain the same curve?

Calculation of the magnetic anisotropy of FePt bulk¶

Hint

Please move to directory Example-3

Note

The original fdf for this example can be found in example_3z.fdf.original. We have made some changes (can you see them?) to speed the calculation while maintaining reasonable quality results.

In this example, you will calculate the magnetic anisotropy bulk FePt.

You will find the example_3z.fdf file and the PPs of Fe and Pt: Fe_fept_SOC.psf and Pt_fept_SOC.psf. Each calculation has to be initialized modifying the DM.InitSpin block, as above, updating the orientation of each atomic magnetization. As a result one has to have different energy values for each (θ,ϕ) pair.

A) Obtain the Etot for three highly symmetric orientations of the initial magnetization, namely, X, Y and Z, keeping fixed the atoms in their unit cell. Again do not use the same density matrix for same calculations. Allow SIESTA to calculate it from scratch.
B) Create a magnetization curve changing the angles at intervals of 20 degrees from the Z-axis initial orientation, to finalize along the X axis. Plot the resulting energy as a function of the varied angle. Check the orientation of the final magnetic moment in each output file. Does it coincide with the initial one?

If this set of calculations takes too much time to be performed completely, you can try to run only a few steps and see how the total energy evolves.

Note

If you do plot the magnetization curve, you might see that it has a minimum at around 60 degrees, which does not make too much sense, and it is quite shallow. In this case we went too far in the simplification of the parameters to gain speed. In particular, removing the polarization orbitals has a very significant effect. At the cost of a noticeably longer simulation time, you can restore those orbitals by using the block:

%Block PAO.Basis
Fe_fept_SOC   2
n=4   0   2  P
    0.0   0.0
n=3   2   2
    0.0   0.0
Pt_fept_SOC   2
n=6     0       2  P
    0.00000     0.00000
n=5     2       2
    0.00000     0.00000
%EndBlock PAO.Basis

You should now get a curve similar to that in Fig. 7:

Fig. 7 Magnetization curve from Z to the X axis¶

Spin-Orbit coupling¶

Performing calculations with different spin configuration.¶

Getting started¶

Simulations¶

Analysis¶

Calculation of the magnetic anisotropy of Pt2.¶

Calculation of the magnetic anisotropy of FePt bulk¶

Calculation of the magnetic anisotropy of Pt₂.¶