In this example we consider iron in the BCC phase. To show that this material is ferromagnetic we will model it once allowing collinear spin polarization and once without and compare the resulting SCF energies. In particular the ground state can only be found if collinear spins are allowed.

First we setup BCC iron without spin polarization using a single iron atom inside the unit cell.

In [1]:

```
using DFTK
a = 5.42352 # Bohr
lattice = a / 2 * [[-1 1 1];
[ 1 -1 1];
[ 1 1 -1]]
atoms = [ElementPsp(:Fe, psp=load_psp("hgh/lda/Fe-q8.hgh"))]
positions = [zeros(3)];
```

To get the ground-state energy we use an LDA model and rather moderate discretisation parameters.

In [2]:

```
kgrid = [3, 3, 3] # k-point grid (Regular Monkhorst-Pack grid)
Ecut = 15 # kinetic energy cutoff in Hartree
model_nospin = model_LDA(lattice, atoms, positions, temperature=0.01)
basis_nospin = PlaneWaveBasis(model_nospin; kgrid, Ecut)
scfres_nospin = self_consistent_field(basis_nospin; tol=1e-4, mixing=KerkerDosMixing());
```

In [3]:

```
scfres_nospin.energies
```

Out[3]:

Energy breakdown (in Ha): Kinetic 15.9208690 AtomicLocal -5.0693382 AtomicNonlocal -5.2202436 Ewald -21.4723040 PspCorrection 1.8758831 Hartree 0.7793484 Xc -3.4467515 Entropy -0.0182879 total -16.650824697571

`Model`

and the guess density functions.
In this case we seek the state with as many spin-parallel
$d$-electrons as possible. In our pseudopotential model the 8 valence
electrons are 2 pair of $s$-electrons, 1 pair of $d$-electrons
and 4 unpaired $d$-electrons giving a desired magnetic moment of `4`

at the iron centre.
The structure (i.e. pair mapping and order) of the `magnetic_moments`

array needs to agree
with the `atoms`

array and `0`

magnetic moments need to be specified as well.

In [4]:

```
magnetic_moments = [4];
```

!!! tip "Units of the magnetisation and magnetic moments in DFTK" Unlike all other quantities magnetisation and magnetic moments in DFTK are given in units of the Bohr magneton $μ_B$, which in atomic units has the value $\frac{1}{2}$. Since $μ_B$ is (roughly) the magnetic moment of a single electron the advantage is that one can directly think of these quantities as the excess of spin-up electrons or spin-up electron density.

We repeat the calculation using the same model as before. DFTK now detects the non-zero moment and switches to a collinear calculation.

In [5]:

```
model = model_LDA(lattice, atoms, positions; magnetic_moments, temperature=0.01)
basis = PlaneWaveBasis(model; Ecut, kgrid)
ρ0 = guess_density(basis, magnetic_moments)
scfres = self_consistent_field(basis, tol=1e-4; ρ=ρ0, mixing=KerkerDosMixing());
```

In [6]:

```
scfres.energies
```

Out[6]:

Energy breakdown (in Ha): Kinetic 16.2946033 AtomicLocal -5.2226744 AtomicNonlocal -5.4099625 Ewald -21.4723040 PspCorrection 1.8758831 Hartree 0.8191819 Xc -3.5406692 Entropy -0.0131622 total -16.669104174398

!!! note "Model and magnetic moments"
DFTK does not store the `magnetic_moments`

inside the `Model`

, but only uses them
to determine the lattice symmetries. This step was taken to keep `Model`

(which contains the physical model) independent of the details of the numerical details
such as the initial guess for the spin density.

In direct comparison we notice the first, spin-paired calculation to be a little higher in energy

In [7]:

```
println("No magnetization: ", scfres_nospin.energies.total)
println("Magnetic case: ", scfres.energies.total)
println("Difference: ", scfres.energies.total - scfres_nospin.energies.total);
```

In [8]:

```
iup = 1
idown = iup + length(scfres.basis.kpoints) ÷ 2
@show scfres.occupation[iup][1:7]
@show scfres.occupation[idown][1:7];
```

Similarly the eigenvalues differ

In [9]:

```
@show scfres.eigenvalues[iup][1:7]
@show scfres.eigenvalues[idown][1:7];
```

`kpoints`

field of the `PlaneWaveBasis`

object contains
each $k$-point coordinate twice, once associated with spin-up and once with down-down.
The list first contains all spin-up $k$-points and then all spin-down $k$-points,
such that `iup`

and `idown`

index the same $k$-point, but differing spins.

In [10]:

```
using Plots
plot_dos(scfres)
```

Out[10]:

Similarly the band structure shows clear differences between both spin components.

In [11]:

```
using Unitful
using UnitfulAtomic
plot_bandstructure(scfres; kline_density=6)
```

Computing bands along kpath: Γ -> H -> N -> Γ -> P -> H and P -> N Diagonalising Hamiltonian kblocks: 100%|████████████████| Time: 0:00:01

Out[11]: