Merge pull request #323 from JuliaMolSim/collinear_dos

DOS and bandstructure calculation for collinear spin

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"]

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))

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,23 +44,23 @@ 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.
 
 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)

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)

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

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.

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
 

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...)