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.
using DFTK
a = 5.42352 # Bohr
lattice = a / 2 * [[-1 1 1];
[ 1 -1 1];
[ 1 1 -1]]
Fe = ElementPsp(:Fe, psp=load_psp("hgh/lda/Fe-q8.hgh"))
atoms = [Fe => [zeros(3)]];
To get the ground-state energy we use an LDA model and rather moderate discretisation parameters.
kgrid = [3, 3, 3] # k-point grid (Regular Monkhorst-Pack grid)
Ecut = 15 # kinetic energy cutoff in Hartree
model_nospin = model_LDA(lattice, atoms, temperature=0.01)
basis_nospin = PlaneWaveBasis(model_nospin; kgrid, Ecut)
scfres_nospin = self_consistent_field(basis_nospin, tol=1e-6, mixing=KerkerMixing());
n Free energy Eₙ-Eₙ₋₁ ρout-ρin α Diag --- --------------- --------- -------- ---- ---- 1 -16.64709508327 NaN 3.28e-01 0.80 5.0 2 -16.64775614862 -6.61e-04 7.80e-02 0.80 1.0 3 -16.64782704334 -7.09e-05 3.22e-03 0.80 2.0 4 -16.64783311599 -6.07e-06 1.74e-03 0.80 5.8 5 -16.64783427919 -1.16e-06 2.47e-04 0.80 2.0 6 -16.64783431093 -3.17e-08 2.04e-05 0.80 2.0
scfres_nospin.energies
Energy breakdown (in Ha): Kinetic 15.9223941 AtomicLocal -5.0701158 AtomicNonlocal -5.2211381 Ewald -21.4723040 PspCorrection 1.8758831 Hartree 0.7795705 Xc -3.4438356 Entropy -0.0182884 total -16.647834310935
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 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.
magnetic_moments = [Fe => [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.
model = model_LDA(lattice, atoms, magnetic_moments=magnetic_moments, temperature=0.01)
basis = PlaneWaveBasis(model; Ecut, kgrid)
ρ0 = guess_density(basis, magnetic_moments)
scfres = self_consistent_field(basis, tol=1e-6; ρ=ρ0, mixing=KerkerMixing());
n Free energy Eₙ-Eₙ₋₁ ρout-ρin Magnet α Diag --- --------------- --------- -------- ------ ---- ---- 1 -16.65996145048 NaN 3.11e-01 2.620 0.80 4.4 2 -16.66625978381 -6.30e-03 6.98e-02 2.428 0.80 1.0 3 -16.66694942782 -6.90e-04 9.41e-03 2.349 0.80 2.8 4 -16.66698695966 -3.75e-05 4.32e-03 2.320 0.80 3.1 5 -16.66699445764 -7.50e-06 1.78e-03 2.309 0.80 2.1 6 -16.66699630484 -1.85e-06 5.57e-04 2.301 0.80 1.5 7 -16.66699651941 -2.15e-07 4.50e-05 2.297 0.80 1.5
scfres.energies
Energy breakdown (in Ha): Kinetic 16.3014939 AtomicLocal -5.2260153 AtomicNonlocal -5.4135655 Ewald -21.4723040 PspCorrection 1.8758831 Hartree 0.8201272 Xc -3.5395383 Entropy -0.0130777 total -16.666996519408
!!! 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
println("No magnetization: ", scfres_nospin.energies.total)
println("Magnetic case: ", scfres.energies.total)
println("Difference: ", scfres.energies.total - scfres_nospin.energies.total);
No magnetization: -16.64783431093486 Magnetic case: -16.666996519407565 Difference: -0.019162208472703668
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, 0.9999988143651348, 0.999998814361033, 0.9999988143563308, 0.9585782956963069, 0.9585777518185276, 1.1544588297395824e-29] (scfres.occupation[idown])[1:7] = [1.0, 0.8324269079406773, 0.8324246161110305, 0.832423305282138, 7.827584139052771e-6, 7.827511182882219e-6, 1.5494575518861517e-32]
Similarly the eigenvalues differ
@show scfres.eigenvalues[iup][1:7]
@show scfres.eigenvalues[idown][1:7];
(scfres.eigenvalues[iup])[1:7] = [-0.06926952258662326, 0.3568787838005498, 0.3568788183965499, 0.3568788580570714, 0.4619146312725758, 0.4619147682481193, 1.1596444537839434] (scfres.eigenvalues[idown])[1:7] = [-0.03040251104243435, 0.477301835769149, 0.47730200006612405, 0.47730209403618556, 0.6109095817708741, 0.6109096749759809, 1.2257792744271487]
!!! note "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
plot_dos(scfres)
Similarly the band structure shows clear differences between both spin components.
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:04