CommonSolve Interface

Project.toml

@@ -1,32 +1,29 @@
 name = "AutoBZCore"
 uuid = "66bd3e16-1600-45cf-8f55-0b550710682b"
 authors = ["lxvm <[email protected]>"]
-version = "0.3.8"
+version = "0.4.0"
 
 [deps]
 AutoSymPTR = "78a0c066-08f1-49a8-82f0-b29cd485e1d3"
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 FourierSeriesEvaluators = "2a892dea-6eef-4bb5-9d1c-de966c9f6db5"
 FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 IteratedIntegration = "3ecdc4d6-ee34-4049-885a-a4e3631db98b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
-HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 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]
 AtomsBaseExt = "AtomsBase"
-HDF5Ext = "HDF5"
 SymmetryReduceBZExt = ["SymmetryReduceBZ", "Polyhedra"]
 UnitfulExt = "Unitful"
 WannierIOExt = "WannierIO"
@@ -36,17 +33,15 @@ Aqua = "0.7"
 AtomsBase = "0.3"
 AutoSymPTR = "0.4"
 ChunkSplitters = "2"
+CommonSolve = "0.2"
 Elliptic = "1"
 FourierSeriesEvaluators = "1"
 FunctionWrappers = "1"
 GeneralizedGaussianQuadrature = "0.1"
 HCubature = "1.4"
-HDF5 = "0.16.15"
 IteratedIntegration = "0.5"
 LinearAlgebra = "1.9"
-Printf = "1.9"
 QuadGK = "2.6"
-Reexport = "1"
 StaticArrays = "1"
 SymmetryReduceBZ = "0.2"
 WannierIO = "0.1,0.2"
@@ -56,9 +51,8 @@ julia = "1.9"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 AtomsBase = "a963bdd2-2df7-4f54-a1ee-49d51e6be12a"
 Elliptic = "b305315f-e792-5b7a-8f41-49f472929428"
-HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 GeneralizedGaussianQuadrature = "958e0c08-f14d-42e8-a0ab-84193b3783f2"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -68,4 +62,4 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 WannierIO = "cb1bc77f-5443-4951-af9f-05b616a3e422"
 
 [targets]
-test = ["Aqua", "Elliptic", "LinearAlgebra", "GeneralizedGaussianQuadrature", "Test", "HDF5", "StaticArrays", "OffsetArrays", "SymmetryReduceBZ", "Polyhedra"]
+test = ["Aqua", "Elliptic", "LinearAlgebra", "GeneralizedGaussianQuadrature", "Test", "StaticArrays", "OffsetArrays", "SymmetryReduceBZ", "Polyhedra"]

aps_example/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 AutoBZCore = "66bd3e16-1600-45cf-8f55-0b550710682b"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+FourierSeriesEvaluators = "2a892dea-6eef-4bb5-9d1c-de966c9f6db5"
 HChebInterp = "78faba9b-a54b-441f-8118-62407cbe4e59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"

aps_example/README.md

@@ -1,4 +1,6 @@
-[Scripts updated for AutoBZCore v0.3]
+[Scripts in development for AutoBZCore v0.4]
+
+Warning: this example may not work well using Julia 1.10. 
 
 Hello, to reproduce the code example in these
 [slides](https://web.mit.edu/lxvm/www/slides/autobz_aps.pdf) follow these steps:

aps_example/aps_example.jl

@@ -20,36 +20,73 @@ for (i, h, n) in zip(hrdat.Rvectors, hrdat.H, hrdat.Rdegens)
     H_R[CartesianIndex(Tuple(i))] = h / n
 end
 
-using AutoBZCore, LinearAlgebra
+using FourierSeriesEvaluators, LinearAlgebra
 
-bz = load_bz(CubicSymIBZ(), "svo.wout")
-# bz = load_bz(IBZ(), "svo.wout") # works with SymmetryReduceBZ.jl installed
 h = FourierSeries(H_R, period=1.0)
 
 η = 1e-2                    # 10 meV (scattering amplitude)
-dos_integrand(h_k::FourierValue, η, ω) = -imag(tr(inv((ω+im*η)*I - h_k.s)))/pi
-integrand = FourierIntegrand(dos_integrand, h, η)
+ω_min = 10
+ω_max = 15
+p0 = (; η, ω=(ω_min + ω_max)/2) # initial parameters
+# BUG cannot redefine this function without breaking functionwrappers in v1.10
+# https://github.com/JuliaLang/julia/issues/52635#issuecomment-2150808569
+greens_function(k, h_k, (; η, ω)) = tr(inv((ω+im*η)*I - h_k))
+prototype = let k = FourierSeriesEvaluators.period(h)
+    greens_function(k, h(k), p0)
+end
+
+using AutoBZCore
+bz = load_bz(CubicSymIBZ(), "svo.wout")
+# bz = load_bz(IBZ(), "svo.wout") # works with SymmetryReduceBZ.jl installed
 
-dos_solver_iai = IntegralSolver(integrand, bz, IAI(); abstol=1e-3)
-dos_solver_ptr = IntegralSolver(integrand, bz, PTR(npt=100); abstol=1e-3)
+integrand = FourierIntegralFunction(greens_function, h, prototype)
+prob_dos = AutoBZProblem(integrand, bz, p0; abstol=1e-3)
 
 using HChebInterp
 
-dos_iai = hchebinterp(dos_solver_iai, 10, 15; atol=1e-2)
-dos_ptr = hchebinterp(dos_solver_ptr, 10, 15; atol=1e-2)
+cheb_order = 15
+
+function dos_solver(prob, alg)
+    solver = init(prob, alg)
+    ω -> begin
+        solver.p = (; solver.p..., ω)
+        solve!(solver).value
+    end
+end
+function threaded_dos_solver(prob, alg; nthreads=min(cheb_order, Threads.nthreads()))
+    solvers = [init(prob, alg) for _ in 1:nthreads]
+    BatchFunction() do ωs
+        out = Vector{typeof(prototype)}(undef, length(ωs))
+        Threads.@threads for i in 1:nthreads
+            solver = solvers[i]
+            for j in i:nthreads:length(ωs)
+                ω = ωs[j]
+                solver.p = (; solver.p..., ω)
+                out[j] = solve!(solver).value
+            end
+        end
+        return out
+    end
+end
+
+dos_solver_iai = dos_solver(prob_dos, IAI(QuadGKJL()))
+@time greens_iai = hchebinterp(dos_solver_iai, ω_min, ω_max; atol=1e-2, order=cheb_order)
+
+dos_solver_ptr = dos_solver(prob_dos, PTR(; npt=100))
+@time greens_ptr = hchebinterp(dos_solver_ptr, ω_min, ω_max; atol=1e-2, order=cheb_order)
 
 using CairoMakie
 
 set_theme!(fontsize=24, linewidth=4)
 
 fig1 = Figure()
-ax1 = Axis(fig1[1,1], limits=((10,15), (0,det(bz.B)*6)), xlabel="ω (eV)", ylabel="SVO DOS (eV⁻¹ Å⁻³)")
-p1 = lines!(ax1, 10:η/100:15, dos_iai; label="IAI(), η=$η")
+ax1 = Axis(fig1[1,1], limits=((10,15), (0,6)), xlabel="ω (eV)", ylabel="SVO DOS (eV⁻¹)")
+p1 = lines!(ax1, 10:η/100:15, ω -> -imag(greens_iai(ω))/pi/det(bz.B); label="IAI, η=$η")
 axislegend(ax1)
 save("iai_svo_dos.pdf", fig1)
 
 fig2 = Figure()
-ax2 = Axis(fig2[1,1], limits=((10,15), (0,det(bz.B)*6)), xlabel="ω (eV)", ylabel="SVO DOS (eV⁻¹ Å⁻³)")
-p2 = lines!(ax2, 10:η/100:15, dos_ptr; label="PTR(), η=$η")
+ax2 = Axis(fig2[1,1], limits=((10,15), (0,6)), xlabel="ω (eV)", ylabel="SVO DOS (eV⁻¹)")
+p2 = lines!(ax2, 10:η/100:15, ω -> -imag(greens_ptr(ω))/pi/det(bz.B); label="PTR, η=$η")
 axislegend(ax2)
 save("ptr_svo_dos.pdf", fig2)

docs/Project.toml

@@ -1,2 +1,5 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FourierSeriesEvaluators = "2a892dea-6eef-4bb5-9d1c-de966c9f6db5"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

docs/make.jl

@@ -23,7 +23,6 @@ makedocs(
         "Algorithms" => "algorithms.md",
         "Reference" => "reference.md",
         "Extensions" => "extensions.md",
-        "Density of States" => "dos.md",
     ],
 )
 

docs/src/algorithms.md

@@ -1,10 +1,12 @@
-# Integral algorithms
+# Algorithms
+
+## `IntegralProblem` algorithms
 
 ```@docs
 AutoBZCore.IntegralAlgorithm
 ```
 
-## Quadrature
+### Quadrature
 
 ```@docs
 AutoBZCore.QuadratureFunction
@@ -14,35 +16,42 @@ AutoBZCore.ContQuadGKJL
 AutoBZCore.MeroQuadGKJL
 ```
 
-## Cubature
+### Cubature
 
 ```@docs
 AutoBZCore.HCubatureJL
 AutoBZCore.MonkhorstPack
 AutoBZCore.AutoSymPTRJL
 ```
 
-## Meta-algorithms
+### Meta-algorithms
 
 ```@docs
 AutoBZCore.NestedQuad
-AutoBZCore.EvalCounter
-AutoBZCore.AbsoluteEstimate
 ```
 
-# BZ-specific integral algorithms
+## `AutoBZProblem` algorithms
 
 In order to make algorithms domain-agnostic, the BZ loaded from
-[`load_bz`](@ref) can be called with the algorithms below, which are wrappers
-for algorithms above with the additional capability of mapping integrals over
-the IBZ to the FBZ.
+[`load_bz`](@ref) can be called with the algorithms below, which are aliases
+for algorithms above
 
 ```@docs
 AutoBZCore.AutoBZAlgorithm
 AutoBZCore.IAI
 AutoBZCore.TAI
 AutoBZCore.PTR
 AutoBZCore.AutoPTR
-AutoBZCore.PTR_IAI
-AutoBZCore.AutoPTR_IAI
-```
+```
+
+## `DOSProblem` algorithms
+
+Currently the available algorithms are an initial release and we would like to include
+the following reference algorithms that are also common in the literature in a future release:
+- (Linear) Tetrahedron Method
+- Adaptive Gaussian broadening
+
+```@docs
+AutoBZCore.DOSAlgorithm
+AutoBZCore.GGR
+```