Collinear spin and magnetic systems

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.

The bulk(:Fe) function from AtomsBuilder returns a BCC iron setup with a single iron atom inside the unit cell.

using AtomsBuilder
using PseudoPotentialData
using DFTK

bulk(:Fe)
FlexibleSystem(Fe, periodicity = TTT):
    cell_vectors      : [  -1.435    1.435    1.435;
                            1.435   -1.435    1.435;
                            1.435    1.435   -1.435]u"Å"

    Atom(Fe, [       0,        0,        0]u"Å")

First we consider a setup without spin polarization. To get the ground-state energy of this system we use an LDA model and rather moderate discretisation parameters.

Ecut  = 15         # kinetic energy cutoff in Hartree
kgrid = [3, 3, 3]  # k-point grid (Regular Monkhorst-Pack grid)
pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")

model_nospin  = model_DFT(bulk(:Fe); pseudopotentials, functionals=LDA(), temperature=0.01)
basis_nospin  = PlaneWaveBasis(model_nospin; kgrid, Ecut)
scfres_nospin = self_consistent_field(basis_nospin; tol=1e-4, mixing=KerkerDosMixing());
┌ Warning: Negative ρcore detected: -0.0002722115697908162
@ DFTK ~/work/DFTK.jl/DFTK.jl/src/terms/xc.jl:39
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -120.8141844411                    0.12    6.2    122ms
  2   -121.8094998534       -0.00       -0.55    3.0    206ms
  3   -121.8262230563       -1.78       -1.39    3.2   24.6ms
  4   -121.8263250251       -3.99       -2.59    1.0   39.1ms
  5   -121.8263297044       -5.33       -3.26    2.0   19.7ms
  6   -121.8263299667       -6.58       -3.71    2.0   20.8ms
  7   -121.8263299750       -8.08       -3.80    1.5   18.3ms
  8   -121.8263299785       -8.47       -4.33    1.0   16.6ms
scfres_nospin.energies
Energy breakdown (in Ha):
    Kinetic             47.8547283
    AtomicLocal         -92.6856069
    AtomicNonlocal      -2.4131245
    Ewald               -85.8893115
    PspCorrection       3.0096982 
    Hartree             24.6937171
    Xc                  -16.3733614
    Entropy             -0.0230693

    total               -121.826329978453

Since we did not specify any initial magnetic moment on the iron atom, DFTK will automatically assume that a calculation with only spin-paired electrons should be performed. As a result the obtained ground state features no spin-polarization.

Now we repeat the calculation, but give the iron atom an initial magnetic moment. For specifying the magnetic moment pass the desired excess of spin-up over spin-down electrons at each centre to the 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 1 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.

magnetic_moments = [4];
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.

model = model_DFT(bulk(:Fe); pseudopotentials, functionals=LDA(),
                  temperature=0.01, magnetic_moments)
basis = PlaneWaveBasis(model; Ecut, kgrid)
ρ0 = guess_density(basis, magnetic_moments)
scfres = self_consistent_field(basis, tol=1e-6; ρ=ρ0, mixing=KerkerDosMixing());
┌ Warning: Negative ρcore detected: -0.0001361057848954081
@ DFTK ~/work/DFTK.jl/DFTK.jl/src/terms/xc.jl:39
n     Energy            log10(ΔE)   log10(Δρ)   Magnet   |Magn|   Diag   Δtime
---   ---------------   ---------   ---------   ------   ------   ----   ------
  1   -120.8704363498                    0.06    1.690    1.691    5.4   63.4ms
  2   -121.8199385062       -0.02       -0.66    2.142    2.147    2.6   47.4ms
  3   -121.8337925540       -1.86       -1.45    2.205    2.217    3.0   64.0ms
  4   -121.8339852324       -3.72       -1.99    2.249    2.264    1.0   32.0ms
  5   -121.8340083253       -4.64       -2.65    2.276    2.293    1.1   33.3ms
  6   -121.8340112866       -5.53       -3.41    2.280    2.298    1.8   36.9ms
  7   -121.8340115157       -6.64       -3.66    2.280    2.298    2.2   41.4ms
  8   -121.8340115495       -7.47       -4.28    2.280    2.298    1.5   36.2ms
  9   -121.8340115491   +   -9.44       -4.32    2.280    2.298    1.2   34.8ms
 10   -121.8340115502       -8.95       -5.49    2.280    2.298    1.5   36.1ms
 11   -121.8340115503      -10.48       -5.62    2.280    2.298    2.9   55.8ms
 12   -121.8340115503      -12.62       -6.09    2.280    2.298    1.1   34.6ms
scfres.energies
Energy breakdown (in Ha):
    Kinetic             47.9532618
    AtomicLocal         -92.7969903
    AtomicNonlocal      -2.4298258
    Ewald               -85.8893115
    PspCorrection       3.0096982 
    Hartree             24.7610456
    Xc                  -16.4309357
    Entropy             -0.0109538

    total               -121.834011550268
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

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

Notice that with the small cutoffs we use to generate the online documentation the calculation is far from converged. With more realistic parameters a larger energy difference of about 0.1 Hartree is obtained.

The spin polarization in the magnetic case is visible if we consider the occupation of the spin-up and spin-down Kohn-Sham orbitals. Especially for the $d$-orbitals these differ rather drastically. For example for the first $k$-point:

iup   = 1
idown = iup + length(scfres.basis.kpoints) ÷ 2
@show scfres.occupation[iup][1:7]
@show scfres.occupation[idown][1:7];
(scfres.occupation[iup])[1:7] = [1.0, 1.0, 1.0, 1.0, 1.0, 0.9999917852729185, 0.9999917852729185]
(scfres.occupation[idown])[1:7] = [1.0, 1.0, 1.0, 1.0, 1.0, 0.986232364236331, 0.9862323632099121]

Similarly the eigenvalues differ

@show scfres.eigenvalues[iup][1:7]
@show scfres.eigenvalues[idown][1:7];
(scfres.eigenvalues[iup])[1:7] = [-3.355215892344667, -1.8664633300762623, -1.8664633234519032, -1.866463323451902, 0.13195509894054513, 0.41745736969727715, 0.4174573696972826]
(scfres.eigenvalues[idown])[1:7] = [-3.2931892739556847, -1.791772799751974, -1.7917727997519648, -1.791772783628266, 0.15434757340667193, 0.4918373939688468, 0.49183739472478444]
``k``-points in collinear calculations

For collinear calculations the 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.

We can observe the spin-polarization by looking at the density of states (DOS) around the Fermi level, where the spin-up and spin-down DOS differ.

using Plots
bands_666 = compute_bands(scfres, MonkhorstPack(6, 6, 6))  # Increase kgrid to get nicer DOS.
plot_dos(bands_666)
Example block output

Note that if same k-grid as SCF should be employed, a simple plot_dos(scfres) is sufficient. We can clearly see that the origin of this spin-polarization traces back to the 3D orbital contribution if we look at the corresponding projected density of states (PDOS).

plot_pdos(bands_666; iatom=1, label="3D")
Example block output

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

using Unitful
using UnitfulAtomic
bands_kpath = compute_bands(scfres; kline_density=6)
plot_bandstructure(bands_kpath)
Example block output