add ImplicitIntegration extension (#14)

* add ImplicitIntegration extension

* add docstring to docs

* bump ImplicitIntegration v0.1.1

Project.toml

@@ -17,12 +17,14 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
+ImplicitIntegration = "bc256489-3a69-4a66-afc4-127cc87e6182"
 Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 SymmetryReduceBZ = "49a35663-c880-4242-bebb-1ec8c0fa8046"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 WannierIO = "cb1bc77f-5443-4951-af9f-05b616a3e422"
 
 [extensions]
+AutoBZCoreImplicitIntegrationExt = "ImplicitIntegration"
 AtomsBaseExt = "AtomsBase"
 SymmetryReduceBZExt = ["SymmetryReduceBZ", "Polyhedra"]
 UnitfulExt = "Unitful"
@@ -39,6 +41,7 @@ FourierSeriesEvaluators = "1"
 FunctionWrappers = "1"
 GeneralizedGaussianQuadrature = "0.1"
 HCubature = "1.4"
+ImplicitIntegration = "0.1.1"
 IteratedIntegration = "0.5"
 LinearAlgebra = "1.9"
 OffsetArrays = "1"
@@ -68,4 +71,4 @@ UnitfulAtomic = "a7773ee8-282e-5fa2-be4e-bd808c38a91a"
 WannierIO = "cb1bc77f-5443-4951-af9f-05b616a3e422"
 
 [targets]
-test = ["Aqua", "Elliptic", "LinearAlgebra", "GeneralizedGaussianQuadrature", "Test", "StaticArrays", "OffsetArrays", "SymmetryReduceBZ", "Polyhedra", "AtomsBase", "Unitful", "UnitfulAtomic", "WannierIO"]
+test = ["Aqua", "Elliptic", "LinearAlgebra", "GeneralizedGaussianQuadrature", "Test", "StaticArrays", "OffsetArrays", "SymmetryReduceBZ", "Polyhedra", "AtomsBase", "Unitful", "UnitfulAtomic", "WannierIO", "ImplicitIntegration"]

docs/src/algorithms.md

@@ -54,4 +54,5 @@ the following reference algorithms that are also common in the literature in a f
 ```@docs
 AutoBZCore.DOSAlgorithm
 AutoBZCore.GGR
+AutoBZCore.ImplicitIntegrationJL
 ```

ext/AutoBZCoreImplicitIntegrationExt.jl

@@ -0,0 +1,41 @@
+module AutoBZCoreImplicitIntegrationExt
+
+using AutoBZCore
+using FourierSeriesEvaluators
+using ImplicitIntegration
+using LinearAlgebra
+using StaticArrays
+
+function AutoBZCore.init_cacheval(h, domain, p, alg::ImplicitIntegrationJL)
+    h isa FourierSeries || throw(ArgumentError("GGR currently supports Fourier series Hamiltonians"))
+    p isa SymmetricBZ || throw(ArgumentError("GGR supports BZ parameters from load_bz"))
+    j = JacobianSeries(h)
+    wj = FourierSeriesEvaluators.workspace_allocate(j, FourierSeriesEvaluators.period(j))
+    return wj
+end
+
+
+function AutoBZCore.dos_solve(h, domain, p, alg::ImplicitIntegrationJL, cacheval;
+    abstol=1e-8, reltol=nothing, maxiters=nothing)
+    domain isa Number || throw(ArgumentError("GGR supports domains of individual eigenvalues"))
+    p isa SymmetricBZ || throw(ArgumentError("GGR supports BZ parameters from load_bz"))
+    E = domain
+    bz = p
+    j = cacheval
+    ϕ = k -> begin
+        H = h(k)-E*I
+        real(StaticArrays._det(StaticArrays.Size(H), H))
+    end
+    f = k -> begin
+        hk, vk = j(k)
+        values, U = eigen(Hermitian(hk))
+        resid, n = findmin(x -> abs(x-E), values)
+        u = U[:,n]
+        return inv(sqrt(mapreduce(v -> abs2(dot(u,v,u)/u'u), +, vk))) # compute 1/|∇ϵₙ|
+    end
+    result = integrate(f, ϕ, map(zero, FourierSeriesEvaluators.period(h)), FourierSeriesEvaluators.period(h); surface=true, tol=abstol)
+    return AutoBZCore.DOSSolution(result.val, AutoBZCore.Success, (; loginfo=result.logger))
+end
+
+
+end

src/AutoBZCore.jl

@@ -93,7 +93,7 @@ include("fourier.jl")
 export DOSProblem
 include("dos_interfaces.jl")
 
-export GGR
+export GGR, ImplicitIntegrationJL
 include("dos_algorithms.jl")
 include("dos_ggr.jl")
 

src/dos_algorithms.jl

@@ -24,3 +24,14 @@ struct GGR <: DOSAlgorithm
     npt::Int
 end
 GGR(; npt=50) = GGR(npt)
+
+
+"""
+    ImplicitIntegrationJL(; kws...)
+
+This algorithm is implemented in an extension. Try it with `using ImplicitIntegration`.
+"""
+struct ImplicitIntegrationJL{K} <: DOSAlgorithm
+    kws::K
+end
+ImplicitIntegrationJL(; kws...) = ImplicitIntegrationJL(kws)

test/dos.jl

@@ -1,6 +1,7 @@
 using Test, AutoBZCore, LinearAlgebra, StaticArrays, OffsetArrays, Elliptic
 using GeneralizedGaussianQuadrature: generalizedquadrature
 using FourierSeriesEvaluators, QuadGK
+using ImplicitIntegration
 
 # test set of known DOS examples
 
@@ -9,8 +10,8 @@ using FourierSeriesEvaluators, QuadGK
 tb_graphene = let t=1.0
     ax = CartesianIndices((-2:2, -2:2))
     hm = OffsetMatrix([zero(MMatrix{2,2,typeof(t),4}) for i in ax], ax)
-    hm[1,1][1,2]   = hm[1,-2][1,2] = hm[-2,1][1,2] = t
-    hm[-1,-1][2,1] = hm[-1,2][2,1] = hm[2,-1][2,1] = t
+    hm[0,0][1,2] = hm[0,-1][1,2] = hm[-1,0][1,2] = t
+    hm[0,0][2,1] = hm[0,1][2,1]  = hm[1,0][2,1]  = t
     FourierSeries(SMatrix{2,2,typeof(t),4}.(hm), period=1.0)
 end
 
@@ -102,9 +103,10 @@ for (model, solution, bandwidth, bzkind) in (
     bz = load_bz(bzkind, I(ndims(model)))
     prob = DOSProblem(model, float(zero(B)), bz)
     E = Float64[-B - 1, -0.8B, -0.6B, -0.2B, 0.1B, 0.3B, 0.5B, 0.7B, 0.9B, B + 2]
-    for alg in (GGR(; npt=200),)
-        cache = AutoBZCore.init(prob, alg)
-        for e in E
+    for alg in (GGR(; npt=200), ImplicitIntegrationJL())
+        cache = AutoBZCore.init(prob, alg; abstol=1e-2)
+        alg isa ImplicitIntegrationJL && !(bzkind isa FBZ) && continue
+        @testset "model=$(nameof(solution)) alg=$alg" for e in E
             cache.domain = e
             @test AutoBZCore.solve!(cache).value ≈ solution(e) atol=1e-2
         end