diff --git a/Project.toml b/Project.toml
index 930b849717..fc00392ccd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,7 +23,6 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PeriodicTable = "7b2266bf-644c-5ea3-82d8-af4bbd25a884"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -55,7 +54,6 @@ NLsolve = "4"
 Optim = "0.22, 1"
 OrderedCollections = "1"
 PeriodicTable = "1"
-Plots = "1"
 Polynomials = "1"
 Primes = "0.4, 0.5"
 ProgressMeter = "1"
@@ -73,8 +71,9 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua", "DoubleFloats", "IntervalArithmetic", "Random", "KrylovKit"]
+test = ["Test", "Aqua", "DoubleFloats", "IntervalArithmetic", "Plots", "Random", "KrylovKit"]
diff --git a/examples/cohen_bergstresser.jl b/examples/cohen_bergstresser.jl
index 84a8c666b7..5b1171c3a7 100644
--- a/examples/cohen_bergstresser.jl
+++ b/examples/cohen_bergstresser.jl
@@ -30,5 +30,6 @@ using Plots
 
 n_bands = 8
 ρ0 = guess_density(basis)  # Just dummy, has no meaning in this model
-p = plot_bandstructure(basis, ρ0, n_bands, εF=εF, kline_density=10)
+ρspin0 = nothing
+p = plot_bandstructure(basis, ρ0, ρspin0, n_bands, εF=εF, kline_density=10)
 ylims!(p, (-5, 6))
diff --git a/examples/collinear_magnetism.jl b/examples/collinear_magnetism.jl
index 4dc563e846..2b1cc20106 100644
--- a/examples/collinear_magnetism.jl
+++ b/examples/collinear_magnetism.jl
@@ -25,8 +25,8 @@ Ecut = 15          # kinetic energy cutoff in Hartree
 model_nospin = model_LDA(lattice, atoms, temperature=0.01)
 basis_nospin = PlaneWaveBasis(model_nospin, Ecut; kgrid=kgrid)
 
-scfres_nospin = self_consistent_field(basis_nospin, tol=1e-6, mixing=KerkerMixing())
-
+scfres_nospin = self_consistent_field(basis_nospin, tol=1e-6, mixing=KerkerMixing());
+#-
 scfres_nospin.energies
 
 # Since we did not specify any initial magnetic moment on the iron atom,
@@ -44,14 +44,14 @@ scfres_nospin.energies
 # 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, ]]
+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 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.
+#     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.
@@ -59,8 +59,8 @@ magnetic_moments = [Fe => [4, ]]
 model = model_LDA(lattice, atoms, magnetic_moments=magnetic_moments, temperature=0.01)
 basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)
 ρspin = guess_spin_density(basis, magnetic_moments)
-scfres = self_consistent_field(basis, tol=1e-10, ρspin=ρspin, mixing=KerkerMixing())
-
+scfres = self_consistent_field(basis, tol=1e-6, ρspin=ρspin, mixing=KerkerMixing());
+#-
 scfres.energies
 
 # !!! note "Model and magnetic moments"
@@ -72,7 +72,8 @@ scfres.energies
 # 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("Magnetic case:    ", scfres.energies.total)
+println("Difference:       ", scfres.energies.total - scfres_nospin.energies.total);
 # 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
@@ -85,15 +86,24 @@ println("Magnetic case:    ", scfres.energies.total);
 
 iup   = 1
 idown = iup + length(scfres.basis.kpoints) ÷ 2
-@show scfres.occupation[iup]
-@show scfres.occupation[idown];
+@show scfres.occupation[iup][1:7]
+@show scfres.occupation[idown][1:7];
 
 # Similarly the eigenvalues differ
-@show scfres.eigenvalues[iup]
-@show scfres.eigenvalues[idown];
+@show scfres.eigenvalues[iup][1:7]
+@show scfres.eigenvalues[idown][1:7];
 
 # !!! 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.
+plot_bandstructure(scfres, kline_density=3, unit=:eV)
diff --git a/examples/metallic_systems.jl b/examples/metallic_systems.jl
index e360678e7b..7aebf574ee 100644
--- a/examples/metallic_systems.jl
+++ b/examples/metallic_systems.jl
@@ -51,9 +51,4 @@ scfres.energies
 
 # The fact that magnesium is a metal is confirmed
 # by plotting the density of states around the Fermi level.
-
-εs = range(minimum(minimum(scfres.eigenvalues)) - .5,
-           maximum(maximum(scfres.eigenvalues)) + .5, length=1000)
-Ds = DOS.(εs, Ref(basis), Ref(scfres.eigenvalues))
-q = plot(εs, Ds, label="DOS")
-vline!(q, [scfres.εF], label="εF")
+plot_dos(scfres)
diff --git a/src/DFTK.jl b/src/DFTK.jl
index d35e2a2a5d..95e0e03109 100644
--- a/src/DFTK.jl
+++ b/src/DFTK.jl
@@ -156,13 +156,13 @@ include("external/pymatgen.jl")
 
 export high_symmetry_kpath
 export compute_bands
-export plot_band_data
 export plot_bandstructure
 include("postprocess/band_structure.jl")
 
 export DOS
 export LDOS
 export NOS
+export plot_dos
 include("postprocess/DOS.jl")
 export compute_χ0
 export apply_χ0
@@ -184,6 +184,7 @@ function __init__()
     @require DoubleFloats="497a8b3b-efae-58df-a0af-a86822472b78" begin
         !isdefined(DFTK, :GENERIC_FFT_LOADED) && include("fft_generic.jl")
     end
+    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plotting.jl")
 end
 
 end # module DFTK
diff --git a/src/PlaneWaveBasis.jl b/src/PlaneWaveBasis.jl
index 9ae4f7a573..0d8bf4377b 100644
--- a/src/PlaneWaveBasis.jl
+++ b/src/PlaneWaveBasis.jl
@@ -209,8 +209,10 @@ Creates a new basis identical to `basis`, but with a different set of kpoints
 """
 function PlaneWaveBasis(basis::PlaneWaveBasis, kcoords::AbstractVector,
                         ksymops::AbstractVector, symmetries=nothing)
+    # TODO This constructor does *not* keep the non-variational property
+    #      of the input basis!
     PlaneWaveBasis(basis.model, basis.Ecut, kcoords, ksymops, symmetries;
-                   fft_size=basis.fft_size)
+                   fft_size=basis.fft_size, variational=true)
 end
 
 
@@ -296,6 +298,15 @@ function index_G_vectors(basis::PlaneWaveBasis, kpoint::Kpoint,
     get(kpoint.mapping_inv, idx_linear, nothing)
 end
 
+"""
+Return the index range of ``k``-points that have a particular spin component.
+"""
+function krange_spin(basis::PlaneWaveBasis, spin::Integer)
+    n_spin = basis.model.n_spin_components
+    @assert 1 ≤ spin ≤ n_spin
+    spinlength = div(length(basis.kpoints), n_spin)
+    (1 + (spin - 1) * spinlength):(spin * spinlength)
+end
 
 #
 # Perform (i)FFTs.
diff --git a/src/external/pymatgen.jl b/src/external/pymatgen.jl
index b5a65c6720..bb144f5d08 100644
--- a/src/external/pymatgen.jl
+++ b/src/external/pymatgen.jl
@@ -37,20 +37,27 @@ function pymatgen_bandstructure(basis, λ, εF, klabels)
     Spin = pyimport("pymatgen.electronic_structure.core").Spin
     bandstructure = pyimport("pymatgen.electronic_structure.bandstructure")
 
-    # This assumes no spin polarization
-    @assert basis.model.spin_polarization in (:none, :spinless)
+    @assert basis.model.spin_polarization in (:none, :spinless, :collinear)
+    n_spin   = basis.model.n_spin_components
+    n_kcoord = div(length(basis.kpoints), n_spin)
 
-    eigenvals_spin_up = Matrix{eltype(λ[1])}(undef, length(λ[1]), length(basis.kpoints))
-    for (ik, λs) in enumerate(λ)
-        eigenvals_spin_up[:, ik] = λs
+    mg_spins = (Spin.up, Spin.down)
+    mg_eigenvals = Dict()
+    for σ in 1:n_spin
+        eigenvals = Matrix{eltype(λ[1])}(undef, length(λ[1]), n_kcoord)
+        for (ito, ik) in enumerate(krange_spin(basis, σ))
+            eigenvals[:, ito] = λ[ik]
+        end
+        mg_eigenvals[mg_spins[σ]] = eigenvals
     end
-    eigenvals = Dict(Spin.up => eigenvals_spin_up)
 
-    kcoords = [kpt.coordinate for kpt in basis.kpoints]
+    # Convert coordinates to Cartesian
+    kcoords = [basis.kpoints[ik].coordinate_cart for ik in 1:n_kcoord]
+    klabels = Dict(lal => basis.model.recip_lattice * vec for (lal, vec) in pairs(klabels))
     pylattice = pymatgen_lattice(basis.model.lattice)
     bandstructure.BandStructureSymmLine(
-        kcoords, eigenvals, pylattice.reciprocal_lattice, εF,
-        labels_dict=klabels, coords_are_cartesian=true
+        kcoords, mg_eigenvals, pylattice.reciprocal_lattice, εF,
+        labels_dict=klabels, coords_are_cartesian=true  # Issues if this is false!
     )
 end
 
diff --git a/src/plotting.jl b/src/plotting.jl
new file mode 100644
index 0000000000..e2f96386a3
--- /dev/null
+++ b/src/plotting.jl
@@ -0,0 +1,88 @@
+import Plots
+
+# This is needed to flag that the plots-dependent code has been loaded
+const PLOTS_LOADED = true
+
+"""
+Plot the trace of an SCF, i.e. the absolute error of the total energy at
+each iteration versus the converged energy in a semilog plot. By default
+a new plot canvas is generated, but an existing one can be passed and reused
+along with `kwargs` for the call to `plot!`.
+"""
+function ScfPlotTrace(plt=Plots.plot(yaxis=:log); kwargs...)
+    energies = Float64[]
+    function callback(info)
+        if info.stage == :finalize
+            minenergy = minimum(energies[max(1, end-5):end])
+            error = abs.(energies .- minenergy)
+            error[error .== 0] .= NaN
+            extra = ifelse(:mark in keys(kwargs), (), (mark=:x, ))
+            Plots.plot!(plt, error; extra..., kwargs...)
+            display(plt)
+        else
+            push!(energies, info.energies.total)
+        end
+    end
+end
+
+
+function plot_band_data(band_data; εF=nothing,
+                        klabels=Dict{String, Vector{Float64}}(), unit=:eV, kwargs...)
+    eshift = isnothing(εF) ? 0.0 : εF
+    data = prepare_band_data(band_data, klabels=klabels)
+
+    # For each branch, plot all bands, spins and errors
+    p = Plots.plot(xlabel="wave vector")
+    for ibranch = 1:data.n_branches
+        kdistances = data.kdistances[ibranch]
+        for spin in data.spins, iband = 1:data.n_bands
+            yerror = nothing
+            if hasproperty(data, :λerror)
+                yerror = data.λerror[ibranch][spin][iband, :] ./ unit_to_au(unit)
+            end
+            energies = (data.λ[ibranch][spin][iband, :] .- eshift) ./ unit_to_au(unit)
+
+            color = (spin == :up) ? :blue : :red
+            Plots.plot!(p, kdistances, energies; color=color, label="", yerror=yerror,
+                        kwargs...)
+        end
+    end
+
+    # X-range: 0 to last kdistance value
+    Plots.xlims!(p, (0, data.kdistances[end][end]))
+    Plots.xticks!(p, data.ticks["distance"],
+                  [replace(l, raw"$\mid$" => " | ") for l in data.ticks["label"]])
+
+    ylims = [-4, 4]
+    !isnothing(εF) && is_metal(band_data, εF) && (ylims = [-10, 10])
+    ylims = round.(ylims * units.eV ./ unit_to_au(unit), sigdigits=2)
+    if isnothing(εF)
+        Plots.ylabel!(p, "eigenvalues  ($(string(unit))")
+    else
+        Plots.ylabel!(p, "eigenvalues - ε_f  ($(string(unit)))")
+        Plots.ylims!(p, ylims...)
+    end
+
+    p
+end
+
+
+function plot_dos(basis, eigenvalues; εF=nothing)
+    n_spin = basis.model.n_spin_components
+    εs = range(minimum(minimum(eigenvalues)) - .5,
+               maximum(maximum(eigenvalues)) + .5, length=1000)
+
+    p = Plots.plot()
+    spinlabels = spin_components(basis.model)
+    colors = [:blue, :red]
+    for σ in 1:n_spin
+        D = DOS.(εs, Ref(basis), Ref(eigenvalues), spins=(σ, ))
+        label = n_spin > 1 ? "DOS $(spinlabels[σ]) spin" : "DOS"
+        Plots.plot!(p, εs, D, label=label, color=colors[σ])
+    end
+    if !isnothing(εF)
+        Plots.vline!(p, [εF], label="εF", color=:green, lw=1.5)
+    end
+    p
+end
+plot_dos(scfres; kwargs...) = plot_dos(scfres.basis, scfres.eigenvalues; εF=scfres.εF, kwargs...)
diff --git a/src/postprocess/DOS.jl b/src/postprocess/DOS.jl
index 5d23831d0f..d6e162e0c1 100644
--- a/src/postprocess/DOS.jl
+++ b/src/postprocess/DOS.jl
@@ -23,12 +23,12 @@ and temperature `T`. It increases with both `T` and better sampling of the BZ wi
 Fermi surface.
 """
 function NOS(ε, basis, eigenvalues; smearing=basis.model.smearing,
-             temperature=basis.model.temperature)
-    @assert basis.model.spin_polarization in (:none, :spinless)
+             temperature=basis.model.temperature, spins=1:basis.model.n_spin_components)
     N = zero(ε)
     if (temperature == 0) || smearing isa Smearing.None
         error("NOS only supports finite temperature")
     end
+    @assert basis.model.spin_polarization in (:none, :spinless, :collinear)
 
     # Note the differences to the DOS and LDOS functions: We are not counting states
     # per BZ volume (like in DOS), but absolute number of states. Therefore n_symeqk
@@ -38,10 +38,10 @@ function NOS(ε, basis, eigenvalues; smearing=basis.model.smearing,
     #
     # To explicitly show the similarity with DOS and the temperature dependence we employ
     # -f'((εik - ε)/temperature) = temperature * ( d/dε f_τ(εik - ε') )|_{ε' = ε}
-    for ik = 1:length(eigenvalues)
+    for σ in spins, ik = krange_spin(basis, σ)
         n_symeqk = length(basis.ksymops[ik])  # Number of symmetry-equivalent k-Points
-        for iband = 1:length(eigenvalues[ik])
-            enred = (eigenvalues[ik][iband] - ε) / temperature
+        for (iband, εnk) in enumerate(eigenvalues[ik])
+            enred = (εnk - ε) / temperature
             N -= n_symeqk * Smearing.occupation_derivative(smearing, enred)
         end
     end
@@ -53,16 +53,17 @@ end
 Total density of states at energy ε
 """
 function DOS(ε, basis, eigenvalues; smearing=basis.model.smearing,
-             temperature=basis.model.temperature)
-    @assert basis.model.spin_polarization in (:none, :spinless)
-    filled_occ = filled_occupation(basis.model)
-    D = zero(ε)
+             temperature=basis.model.temperature, spins=1:basis.model.n_spin_components)
     if (temperature == 0) || smearing isa Smearing.None
         error("DOS only supports finite temperature")
     end
-    for ik = 1:length(eigenvalues)
-        for iband = 1:length(eigenvalues[ik])
-            enred = (eigenvalues[ik][iband] - ε) / temperature
+    @assert basis.model.spin_polarization in (:none, :spinless, :collinear)
+    filled_occ = filled_occupation(basis.model)
+
+    D = zero(ε)
+    for σ in spins, ik = krange_spin(basis, σ)
+        for (iband, εnk) in enumerate(eigenvalues[ik])
+            enred = (εnk - ε) / temperature
             D -= (filled_occ * basis.kweights[ik] / temperature
                   * Smearing.occupation_derivative(smearing, enred))
         end
@@ -74,16 +75,17 @@ end
 Local density of states, in real space
 """
 function LDOS(ε, basis, eigenvalues, ψ; smearing=basis.model.smearing,
-              temperature=basis.model.temperature)
-    @assert basis.model.spin_polarization in (:none, :spinless)
-    filled_occ = filled_occupation(basis.model)
+              temperature=basis.model.temperature, spins=1:basis.model.n_spin_components)
     if (temperature == 0) || smearing isa Smearing.None
         error("LDOS only supports finite temperature")
     end
+    @assert basis.model.spin_polarization in (:none, :spinless, :collinear)
+    filled_occ = filled_occupation(basis.model)
+
     weights = deepcopy(eigenvalues)
-    for ik = 1:length(eigenvalues)
-        for iband = 1:length(eigenvalues[ik])
-            enred = (eigenvalues[ik][iband] - ε) / temperature
+    for (ik, εk) in enumerate(eigenvalues)
+        for (iband, εnk) in enumerate(εk)
+            enred = (εnk - ε) / temperature
             weights[ik][iband] = (-filled_occ / temperature
                                   * Smearing.occupation_derivative(smearing, enred))
         end
@@ -93,5 +95,18 @@ function LDOS(ε, basis, eigenvalues, ψ; smearing=basis.model.smearing,
     # weights (as "occupations") at each kpoint. Note, that this automatically puts in the
     # required symmetrization with respect to kpoints and BZ symmetry
     ldostot, ldosspin = compute_density(basis, ψ, weights)
-    ldostot.real
+
+    # TODO This is not great, make compute_density more flexible ...
+    if basis.model.spin_polarization == :collinear
+        ρs = [(ldostot.real + ldosspin.real) / 2, (ldostot.real - ldosspin.real) / 2]
+    else
+        @assert isnothing(ldosspin)
+        ρs = [ldostot.real]
+    end
+    return sum(ρs[iσ] for iσ in spins)
 end
+
+"""
+Plot the density of states over a reasonable range
+"""
+function plot_dos end
diff --git a/src/postprocess/band_structure.jl b/src/postprocess/band_structure.jl
index 5adaa8fbd3..e1ab6669b6 100644
--- a/src/postprocess/band_structure.jl
+++ b/src/postprocess/band_structure.jl
@@ -1,5 +1,4 @@
 using PyCall
-import Plots
 
 # Functionality for computing band structures, mostly using pymatgen
 
@@ -19,7 +18,7 @@ function high_symmetry_kpath(model; kline_density=20)
     (kcoords=kcoords, klabels=labels_dict, kpath=symm_kpath.kpath["path"])
 end
 
-@timing function compute_bands(basis, ρ, kcoords, n_bands;
+@timing function compute_bands(basis, ρ, ρspin, kcoords, n_bands;
                                eigensolver=lobpcg_hyper,
                                tol=1e-3,
                                show_progress=true,
@@ -30,7 +29,7 @@ end
     myrationalize(x::T) where {T <: AbstractFloat} = rationalize(x, tol=10eps(T))
     myrationalize(x) = x
     bs_basis = PlaneWaveBasis(basis, [myrationalize.(k) for k in kcoords], ksymops)
-    ham = Hamiltonian(bs_basis; ρ=ρ)
+    ham = Hamiltonian(bs_basis; ρ=ρ, ρspin=ρspin)
 
     band_data = diagonalize_all_kblocks(eigensolver, ham, n_bands + 3;
                                         n_conv_check=n_bands,
@@ -47,19 +46,21 @@ function prepare_band_data(band_data; datakeys=[:λ, :λerror],
     # from the band_data object into nicely plottable branches
     # This is a bit of abuse of the routines in pymatgen, but it works ...
     plotter = pyimport("pymatgen.electronic_structure.plotter")
+    basis = band_data.basis
 
-    ret = Dict{Symbol, Any}(:basis => band_data.basis)
+    ret = Dict{Symbol, Any}(:basis => basis)
     for key in datakeys
         hasproperty(band_data, key) || continue
 
         # Compute dummy "Fermi level" for pymatgen to be happy
         allfinite = [filter(isfinite, x) for x in band_data[key]]
         eshift = sum(sum, allfinite) / sum(length, allfinite)
-        bs = pymatgen_bandstructure(band_data.basis, band_data[key], eshift, klabels)
+        bs = pymatgen_bandstructure(basis, band_data[key], eshift, klabels)
         data = plotter.BSPlotter(bs).bs_plot_data(zero_to_efermi=false)
 
         # Check number of k-Points agrees
-        @assert length(band_data.basis.kpoints) == sum(length, data["distances"])
+        n_kcoords = div(length(basis.kpoints), basis.model.n_spin_components)
+        @assert n_kcoords == sum(length, data["distances"])
 
         ret[:spins] = [:up]
         spinmap = [("1", :up)]
@@ -88,7 +89,7 @@ i.e. where bands are cut by the Fermi level.
 """
 function is_metal(band_data, εF, tol=1e-4)
     # This assumes no spin polarization
-    @assert band_data.basis.model.spin_polarization in (:none, :spinless)
+    @assert band_data.basis.model.spin_polarization in (:none, :spinless, :collinear)
 
     n_bands = length(band_data.λ[1])
     n_kpoints = length(band_data.λ)
@@ -100,47 +101,6 @@ function is_metal(band_data, εF, tol=1e-4)
     false
 end
 
-
-function plot_band_data(band_data; εF=nothing,
-                        klabels=Dict{String, Vector{Float64}}(), unit=:eV, kwargs...)
-    eshift = isnothing(εF) ? 0.0 : εF
-    data = prepare_band_data(band_data, klabels=klabels)
-
-    # For each branch, plot all bands, spins and errors
-    p = Plots.plot(xlabel="wave vector")
-    for ibranch = 1:data.n_branches
-        kdistances = data.kdistances[ibranch]
-        for spin in data.spins, iband = 1:data.n_bands
-            yerror = nothing
-            if hasproperty(data, :λerror)
-                yerror = data.λerror[ibranch][spin][iband, :] ./ unit_to_au(unit)
-            end
-            energies = (data.λ[ibranch][spin][iband, :] .- eshift) ./ unit_to_au(unit)
-
-            color = (spin == :up) ? :blue : :red
-            Plots.plot!(p, kdistances, energies; color=color, label="", yerror=yerror,
-                        kwargs...)
-        end
-    end
-
-    # X-range: 0 to last kdistance value
-    Plots.xlims!(p, (0, data.kdistances[end][end]))
-    Plots.xticks!(p, data.ticks["distance"],
-                  [replace(l, raw"$\mid$" => " | ") for l in data.ticks["label"]])
-
-    ylims = [-4, 4]
-    !isnothing(εF) && is_metal(band_data, εF) && (ylims = [-10, 10])
-    ylims = round.(ylims * units.eV ./ unit_to_au(unit), sigdigits=2)
-    if isnothing(εF)
-        Plots.ylabel!(p, "eigenvalues  ($(string(unit))")
-    else
-        Plots.ylabel!(p, "eigenvalues - ε_f  ($(string(unit)))")
-        Plots.ylims!(p, ylims...)
-    end
-
-    p
-end
-
 function detexify_kpoint(string)
     # For some reason Julia doesn't support this naively: https://github.com/JuliaLang/julia/issues/29849
     replacements = ("\\Gamma" => "Γ",
@@ -159,23 +119,28 @@ If this value is absent and an `scfres` is used to start the calculation a defau
 `n_bands_scf + 5sqrt(n_bands_scf)` is used. Unlike the rest of DFTK bands energies
 are plotted in `:eV` unless a different `unit` is selected.
 """
-function plot_bandstructure(basis, ρ, n_bands;
+function plot_bandstructure(basis, ρ, ρspin, n_bands;
                             εF=nothing, kline_density=20, unit=:eV, kwargs...)
+    if !isdefined(DFTK, :PLOTS_LOADED)
+        error("Plots not loaded. Run 'using Plots' before calling plot_bandstructure.")
+    end
+
     # Band structure calculation along high-symmetry path
     kcoords, klabels, kpath = high_symmetry_kpath(basis.model; kline_density=kline_density)
     println("Computing bands along kpath:")
     println("       ", join(join.(detexify_kpoint.(kpath), " -> "), "  and  "))
-    band_data = compute_bands(basis, ρ, kcoords, n_bands; kwargs...)
+    band_data = compute_bands(basis, ρ, ρspin, kcoords, n_bands; kwargs...)
 
     plotargs = ()
     if kline_density ≤ 10
         plotargs = (markersize=2, markershape=:circle)
     end
+
     plot_band_data(band_data; εF=εF, klabels=klabels, unit=unit, plotargs...)
 end
 function plot_bandstructure(scfres; n_bands=nothing, kwargs...)
     # Convenience wrapper for scfres named tuples
     n_bands_scf = size(scfres.occupation[1], 2)
     isnothing(n_bands) && (n_bands = ceil(Int, n_bands_scf + 5sqrt(n_bands_scf)))
-    plot_bandstructure(scfres.basis, scfres.ρ, n_bands; εF=scfres.εF, kwargs...)
+    plot_bandstructure(scfres.basis, scfres.ρ, scfres.ρspin, n_bands; εF=scfres.εF, kwargs...)
 end
diff --git a/src/scf/scf_callbacks.jl b/src/scf/scf_callbacks.jl
index 2ba08cb5c3..21155d9ae6 100644
--- a/src/scf/scf_callbacks.jl
+++ b/src/scf/scf_callbacks.jl
@@ -1,24 +1,5 @@
-"""
-Plot the trace of an SCF, i.e. the absolute error of the total energy at
-each iteration versus the converged energy in a semilog plot. By default
-a new plot canvas is generated, but an existing one can be passed and reused
-along with `kwargs` for the call to `plot!`.
-"""
-function ScfPlotTrace(plt=plot(yaxis=:log); kwargs...)
-    energies = Float64[]
-    function callback(info)
-        if info.stage == :finalize
-            minenergy = minimum(energies[max(1, end-5):end])
-            error = abs.(energies .- minenergy)
-            error[error .== 0] .= NaN
-            extra = ifelse(:mark in keys(kwargs), (), (mark=:x, ))
-            plot!(plt, error; extra..., kwargs...)
-            display(plt)
-        else
-            push!(energies, info.energies.total)
-        end
-    end
-end
+# For ScfPlotTrace() see DFTK.jl/src/plotting.jl, which is conditionally loaded upon
+# Plots.jl is included.
 
 """
 Default callback function for `self_consistent_field`, which prints a convergence table
diff --git a/src/scf/self_consistent_field.jl b/src/scf/self_consistent_field.jl
index d5892de0fa..10bd28a100 100644
--- a/src/scf/self_consistent_field.jl
+++ b/src/scf/self_consistent_field.jl
@@ -1,4 +1,3 @@
-using Plots
 include("scf_callbacks.jl")
 
 function default_n_bands(model)
diff --git a/test/compute_bands.jl b/test/compute_bands.jl
index ff23e8d674..9de1fc2a05 100644
--- a/test/compute_bands.jl
+++ b/test/compute_bands.jl
@@ -140,13 +140,14 @@ end
 
     # Build Hamiltonian just from SAD guess
     ρ0 = guess_density(basis, [spec => testcase.positions])
-    ham = Hamiltonian(basis; ρ=ρ0)
+    ρspin0 = nothing
+    ham = Hamiltonian(basis; ρ=ρ0, ρspin=ρspin0)
 
     # Check that plain diagonalization and compute_bands agree
     eigres = diagonalize_all_kblocks(lobpcg_hyper, ham, n_bands + 3, n_conv_check=n_bands,
                                      tol=1e-5)
 
-    band_data = compute_bands(basis, ρ0, [k.coordinate for k in basis.kpoints], n_bands)
+    band_data = compute_bands(basis, ρ0, ρspin0, [k.coordinate for k in basis.kpoints], n_bands)
     for ik in 1:length(basis.kpoints)
         @test eigres.λ[ik][1:n_bands] ≈ band_data.λ[ik] atol=1e-5
     end
@@ -196,14 +197,16 @@ end
         @test ret.λerror[3][:up][iband, :] == ret.λ[3][:up][iband, :] ./ 100
     end
 
+    B = model.recip_lattice
     ref_distances = zeros(3, 3)  # row idx is kpoint, col idx is branch,
     ikpt = 1
     for ibr in 1:3
         ibr != 1 && (ref_distances[1, ibr] = ref_distances[end, ibr-1])
         ikpt += 1
         for ik in 2:3
-            ref_distances[ik, ibr] = (ref_distances[ik-1, ibr]
-                                      + norm(kcoords[ikpt-1] - kcoords[ikpt]))
+            ref_distances[ik, ibr] = (
+                ref_distances[ik-1, ibr] + norm(B * (kcoords[ikpt-1] - kcoords[ikpt]))
+            )
             ikpt += 1
         end
     end