Merge branch 'master' into kgrid

examples/anyons.jl

@@ -0,0 +1,34 @@
+# # Anyons
+
+# We solve the almost-bosonic anyon model of https://arxiv.org/pdf/1901.10739.pdf
+
+using DFTK
+using StaticArrays
+using Plots
+
+# Unit cell. Having one of the lattice vectors as zero means a 2D system
+a = 14
+lattice = a .* [[1 0 0.]; [0 1 0]; [0 0 0]];
+
+# Confining scalar potential
+pot(x, y, z) = ((x - a/2)^2 + (y - a/2)^2)
+
+# Parameters
+Ecut = 50
+n_electrons = 1
+β = 5
+
+terms = [Kinetic(2),
+         ExternalFromReal(X -> pot(X...)),
+         Anyonic(1, β)
+]
+model = Model(lattice; n_electrons=n_electrons,
+              terms=terms, spin_polarization=:spinless)  # "spinless electrons"
+basis = PlaneWaveBasis(model, Ecut, kgrid=(1, 1, 1);
+                       fft_size=DFTK.compute_fft_size(lattice, Ecut, supersampling=1.2))
+scfres = direct_minimization(basis, tol=1e-14)  # Reduce tol for production
+E = scfres.energies.total
+s = 2
+E11 = π/2 * (2(s+1)/s)^((s+2)/s) * (s/(s+2))^(2(s+1)/s) * E^((s+2)/s) / β
+println("e(1,1) / (2π)= ", E11 / (2π))
+display(heatmap(scfres.ρ.real[:, :, 1], c=:blues))

src/eigen/preconditioners.jl

@@ -33,7 +33,8 @@ mutable struct PreconditionerTPA{T <: Real}
 end
 
 function PreconditionerTPA(basis::PlaneWaveBasis{T}, kpt::Kpoint; default_shift=1) where T
-    kin = Vector{T}([sum(abs2, G + kpt.coordinate_cart)
+    scaling = only([t for t in basis.model.term_types if t isa Kinetic]).scaling_factor
+    kin = Vector{T}([scaling * sum(abs2, G + kpt.coordinate_cart)
                      for G in G_vectors_cart(kpt)] ./ 2)
     @assert basis.model.spin_polarization in (:none, :collinear, :spinless)
     PreconditionerTPA{T}(basis, kpt, kin, nothing, default_shift)

src/terms/anyonic.jl

@@ -2,14 +2,63 @@
 # This term does not contain the kinetic energy, which must be added separately
 # /!\ They have no 1/2 factor in front of the kinetic energy,
 #     so for consistency the added kinetic energy must have a scaling of 2
-# Energy = <u, ((-i∇ + βA)^2 + V) u>
+# Energy = <u, ((-i hbar ∇ + βA)^2 + V) u>
 # where ∇∧A = 2π ρ, ∇⋅A = 0 => A = x^⟂/|x|² ∗ ρ
-# H = (-i∇ + βA)² + V - 2β x^⟂/|x|² ∗ (βAρ + J)
-#   = -Δ + 2β (-i∇)⋅A + β²|A|^2 - 2β x^⟂/|x|² ∗ (βAρ + J)
+# H = (-i hbar ∇ + βA)² + V - 2β x^⟂/|x|² ∗ (βAρ + J)
+#   = - hbar^2 Δ + 2β hbar (-i∇)⋅A + β²|A|^2 - 2β x^⟂/|x|² ∗ (βAρ + J)
 # where only the first three terms "count for the energy", and where
-# J = 1/(2i) (u* ∇u - u ∇u*)
+# J = 1/(2i) hbar (u* ∇u - u ∇u*)
+
+# for numerical reasons, we solve the equation ∇∧A = 2π ρ in two parts: A = A_SR + A_LR, with
+# ∇∧A_SR = 2π(ρ-ρref)
+# ∇∧A_LR = 2πρref
+# where ρref is a gaussian centered at the origin and with the same integral as ρ (=M)
+# This way, ρ-ρref has zero mass, and A_SR is shorter range.
+# The long-range part is computed explicitly
+
+function ρref_real(x::T, y::T, M, σ) where {T <: Real}
+    r = hypot(x, y)
+    # gaussian normalized to have integral M
+    M * exp(-T(1)/2 * (r/σ)^2) / (σ^2 * (2T(π))^(2/2))
+end
+
+function magnetic_field_produced_by_ρref(x::T, y::T, M, σ) where {T <: Real}
+    # The solution of ∇∧A = 2π ρref is ϕ(r) [-y;x] where ϕ satisfies
+    # ∇∧A = 2ϕ + r ϕ' => rϕ' + 2 ϕ = 2π ρref
+    # Wolfram alpha (after a bit of coaxing) says the solution to
+    # rϕ' + 2ϕ = C exp(-α x^2)
+    # is
+    # 1/x^2 (c1 - C exp(-α x^2)/2α), which coupled with smoothness at 0 gives
+    # C/(2α)/x^2*(1 - exp(-α x^2))
+    r = hypot(x, y)
+    r == 0 && return magnetic_field_produced_by_ρref(1e-8, 0.0, M, σ) # hack
+
+    α = 1 / (2*σ^2)
+    C = 2T(π) * M / (σ^2 * (2T(π))^(2/2))
+    r = hypot(x, y)
+    ϕ = 1/2 * C / α / r^2 * (1 - exp(-α*r^2))
+    ϕ * @SVector[-y, x]
+end
+
+function make_div_free(basis::PlaneWaveBasis{T}, A) where {T}
+    out = [zeros(complex(T), basis.fft_size), zeros(complex(T), basis.fft_size)]
+    A_fourier = [from_real(basis, A[α]).fourier for α = 1:2]
+    for (iG, G) in enumerate(G_vectors(basis))
+        vec = [A_fourier[1][iG], A_fourier[2][iG]]
+        G = [G[1], G[2]]
+        # project on divergence free fields, ie in Fourier
+        # project A(G) on the orthogonal of G
+        if iG != 1
+            out[1][iG], out[2][iG] = vec - (G'vec) * G / sum(abs2, G)
+        else
+            out[1][iG], out[2][iG] = vec
+        end
+    end
+    [from_fourier(basis, out[α]).real for α = 1:2]
+end
 
 struct Anyonic
+    hbar
     β
 end
 function (A::Anyonic)(basis)
@@ -22,17 +71,38 @@ function (A::Anyonic)(basis)
     @assert basis.model.lattice[2, 1] == basis.model.lattice[1, 2] == 0
     @assert basis.model.lattice[1, 1] == basis.model.lattice[2, 2]
 
-    TermAnyonic(basis, eltype(basis)(A.β))
+    TermAnyonic(basis, eltype(basis)(A.hbar), eltype(basis)(A.β))
 end
 
-struct TermAnyonic{T <: Real} <: Term
+struct TermAnyonic{T <: Real, Tρ <: RealFourierArray, TA} <: Term
     basis::PlaneWaveBasis{T}
+    hbar::T
     β::T
+    ρref::Tρ
+    Aref::TA
+end
+function TermAnyonic(basis::PlaneWaveBasis{T}, hbar, β) where {T}
+    # compute correction magnetic field
+    ρref = zeros(T, basis.fft_size)
+    Aref = [zeros(T, basis.fft_size), zeros(T, basis.fft_size)]
+    M = basis.model.n_electrons  # put M to zero to disable
+    σ = 2
+    for (ir, r) in enumerate(r_vectors(basis))
+        rcart = basis.model.lattice * (r - @SVector[.5, .5, .0])
+        ρref[ir] = ρref_real(rcart[1], rcart[2], M, σ)
+        Aref[1][ir], Aref[2][ir] = magnetic_field_produced_by_ρref(rcart[1], rcart[2], M, σ)
+    end
+    # Aref is not divergence-free in the finite basis, so we explicitly project it
+    # This is because we assume elsewhere (eg to ensure self-adjointness of the Hamiltonian)
+    Aref = make_div_free(basis, Aref)
+    ρref = from_real(basis, ρref)
+    TermAnyonic(basis, hbar, β, ρref, Aref)
 end
 
 function ene_ops(term::TermAnyonic, ψ, occ; ρ, kwargs...)
     basis = term.basis
     T = eltype(basis)
+    hbar = term.hbar
     β = term.β
     @assert ψ !== nothing # the hamiltonian depends explicitly on ψ
 
@@ -42,26 +112,28 @@ function ene_ops(term::TermAnyonic, ψ, occ; ρ, kwargs...)
     # => A(G1, G2) = -2π i ρ(G1, G2) * [-G2;G1;0]/(G1^2 + G2^2)
     A1 = zeros(complex(T), basis.fft_size)
     A2 = zeros(complex(T), basis.fft_size)
+    ρ_fourier = ρ.fourier  # unpack before hot loop for perf
+    ρref_fourier = term.ρref.fourier
     for (iG, G) in enumerate(G_vectors_cart(basis))
         G2 = sum(abs2, G)
         if G2 != 0
-            A1[iG] = +2T(π) * G[2] / G2 * ρ.fourier[iG] * im
-            A2[iG] = -2T(π) * G[1] / G2 * ρ.fourier[iG] * im
+            A1[iG] = +2T(π) * G[2] / G2 * (ρ_fourier[iG] - ρref_fourier[iG]) * im
+            A2[iG] = -2T(π) * G[1] / G2 * (ρ_fourier[iG] - ρref_fourier[iG]) * im
         end
     end
-    Areal = [from_fourier(basis, A1).real,
-             from_fourier(basis, A2).real,
+    Areal = [from_fourier(basis, A1).real + term.Aref[1],
+             from_fourier(basis, A2).real + term.Aref[2],
              zeros(T, basis.fft_size)]
 
-    # 2β (-i∇)⋅A + β^2 |A|^2
+    # 2 hbar β (-i∇)⋅A + β^2 |A|^2
     ops_energy = [MagneticFieldOperator(basis, basis.kpoints[1],
-                                        2β .* Areal),
+                                        2 .* hbar .* β .* Areal),
                   RealSpaceMultiplication(basis, basis.kpoints[1],
                                           β^2 .* (abs2.(Areal[1]) .+ abs2.(Areal[2])))]
 
     # Now compute effective local potential - 2β x^⟂/|x|² ∗ (βAρ + J)
     J = compute_current(basis, ψ, occ)
-    eff_current = [J[α].real .+
+    eff_current = [hbar .* J[α].real .+
                    β .* ρ.real .* Areal[α] for α = 1:2]
     eff_current_fourier = [from_real(basis, eff_current[α]).fourier for α = 1:2]
     # eff_pot = - 2β x^⟂/|x|² ∗ eff_current

test/anyons.jl

@@ -0,0 +1,20 @@
+using DFTK
+using LinearAlgebra
+using Test
+
+@testset "Anyons" begin
+    # Test that the magnetic field satisfies ∇∧A = 2π ρref, ∇⋅A = 0
+    x = 1.23
+    y = -1.8
+    ε = 1e-8
+    M = 2.31
+    σ = 1.81
+    dAdx = (  DFTK.magnetic_field_produced_by_ρref(x + ε, y, M, σ)
+            - DFTK.magnetic_field_produced_by_ρref(x,     y, M, σ)) / ε
+    dAdy = (  DFTK.magnetic_field_produced_by_ρref(x, y + ε, M, σ)
+            - DFTK.magnetic_field_produced_by_ρref(x, y,     M, σ)) / ε
+    curlA = dAdx[2] - dAdy[1]
+    divA = dAdx[1] + dAdy[2]
+    @test norm(curlA - 2π*DFTK.ρref_real(x, y, M, σ)) < 1e-4
+    @test abs(divA) < 1e-6
+end

test/hamiltonian_consistency.jl

@@ -69,6 +69,8 @@ let
     pot(x, y, z) = (x - a/2)^2 + (y - a/2)^2
     Apot(x, y, z) = .2 * [y - a/2, -(x - a/2), 0]
     Apot(X) = Apot(X...)
-    test_consistency_term(Magnetic(Apot);  kgrid=[1, 1, 1], lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
-    test_consistency_term(DFTK.Anyonic(1); kgrid=[1, 1, 1], lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
+    test_consistency_term(Magnetic(Apot);
+                          kgrid=[1, 1, 1], lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
+    test_consistency_term(DFTK.Anyonic(2, 3.2);
+                          kgrid=[1, 1, 1], lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
 end

test/runtests.jl

@@ -86,6 +86,7 @@ Random.seed!(0)
 
     if "all" in TAGS
         include("ewald.jl")
+        include("anyons.jl")
         include("energy_nuclear.jl")
         include("occupation.jl")
         include("energies_guess_density.jl")