Significant overhaul to fix #5:

- make `Cell`, `KPath`, and `KPathInterpolant` consistently return points in the lattice basis, and facilitate the conversion for plotting methods automatically
- drop `Rs::nothing` constructor for `irrfbz_path`
- bump major version number (v0.4.0)
- update tests and docs
- add copied ("vendored") pieces from Crystalline to enable conversion to primitive reciprocal basis

Project.toml

@@ -1,7 +1,7 @@
 name = "Brillouin"
 uuid = "23470ee3-d0df-4052-8b1a-8cbd6363e7f0"
 authors = ["Thomas Christensen <tchr@mit.edu> and contributors"]
-version = "0.3.2"
+version = "0.4.0"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

README.md

@@ -18,36 +18,38 @@ julia> Rs   = ([1.0, 0.0, 0.0], [-0.5, sqrt(3)/2, 0.0],   [0, 0, 1.25]) # or `Cr
 julia> Gs   = reciprocalbasis(Rs)
 julia> cell = wignerseitz(Gs)
 Cell{3} (8 faces, 12 vertices):
- verts: [4.18879, 0.0, -2.513274]
-        [2.094395, -3.627599, -2.513274]
-        [4.18879, 0.0, 2.513274]
-        [2.094395, 3.627599, 2.513274]
-        [2.094395, 3.627599, -2.513274]
-        [2.094395, -3.627599, 2.513274]
-        [-2.094395, 3.627599, -2.513274]
-        [-4.18879, 0.0, -2.513274]
-        [-2.094395, -3.627599, -2.513274]
-        [-2.094395, -3.627599, 2.513274]
-        [-4.18879, 0.0, 2.513274]
-        [-2.094395, 3.627599, 2.513274]
- faces: [5, 4, 3, 1]
-        [8, 9, 10, 11]
-        [2, 1, 3, 6]
-        [2, 6, 10, 9]
-        [7, 5, 1, 2, 9, 8]
-        [4, 12, 11, 10, 6, 3]
-        [4, 5, 7, 12]
-        [11, 12, 7, 8]
- basis: [6.283185, 3.627599, -0.0]
-        [0.0, 7.255197, 0.0]
-        [0.0, -0.0, 5.026548]
+  verts: [0.666667, -0.333333, -0.5]
+         [0.333333, -0.666667, -0.5]
+         [0.666667, -0.333333, 0.5]
+         [0.333333, 0.333333, 0.5]
+         [0.333333, 0.333333, -0.5]
+         [0.333333, -0.666667, 0.5]
+         [-0.333333, 0.666667, -0.5]
+         [-0.666667, 0.333333, -0.5]
+         [-0.333333, -0.333333, -0.5]
+         [-0.333333, -0.333333, 0.5]
+         [-0.666667, 0.333333, 0.5]
+         [-0.333333, 0.666667, 0.5]
+  faces: [5, 4, 3, 1]
+         [8, 9, 10, 11]
+         [2, 1, 3, 6]
+         [2, 6, 10, 9]
+         [7, 5, 1, 2, 9, 8]
+         [4, 12, 11, 10, 6, 3]
+         [4, 5, 7, 12]
+         [11, 12, 7, 8]
+  basis: [6.283185, 3.627599, -0.0]
+         [0.0, 7.255197, 0.0]
+         [0.0, -0.0, 5.026548]
 ```
-The resulting Brillouin zone can be plotted using e.g. [PlotlyJS.jl](https://github.com/JuliaPlots/PlotlyJS.jl) (or 3D-capable backends of [AbstractPlotting.jl](https://github.com/JuliaPlots/AbstractPlotting.jl) such as [GLMakie.jl](https://github.com/JuliaPlots/GLMakie.jl)):
+The returned vertices are in the coordiantes of the provided reciprocal basis (to convert, see `cartesianize(!)`); this is the default behavior in Brillouin.
+
+The Brillouin zone can be plotted using e.g. [PlotlyJS.jl](https://github.com/JuliaPlots/PlotlyJS.jl) (or 3D-capable backends of [AbstractPlotting.jl](https://github.com/JuliaPlots/AbstractPlotting.jl) such as [GLMakie.jl](https://github.com/JuliaPlots/GLMakie.jl)):
 ```jl
 julia> using PlotlyJS
 julia> plot(cell)
 ```
-Examples of interactive visualizations are [included in the documentation](https://thchr.github.io/Brillouin.jl/stable/wignerseitz/).
+Examples of interactive visualizations are [included in the documentation](https://thchr.github.io/Brillouin.jl/stable/wignerseitz/). Visualizations will automatically be displayed in Cartesian coordinates.
 
 Irreducible **k**-paths are returned by `irrfbz_path`, and can similarly be visualized (see [examples in documentation](https://thchr.github.io/Brillouin.jl/stable/kpaths/)):
 ```jl
@@ -63,5 +65,8 @@ KPath{3} (7 points, 3 paths, 13 points in paths):
   paths: [:Γ, :M, :K, :Γ, :A, :L, :H, :A]
          [:L, :M]
          [:H, :K, :H₂]
+  basis: [6.283185, 3.627599, -0.0]
+         [0.0, 7.255197, 0.0]
+         [0.0, -0.0, 5.026548]
 ```
 The resulting object can be interpolated, using either `interpolate(kp, N)` or `splice(kp, N)`, which produces an `KPathInterpolant` iterable whose elements interpolate the connected paths (and enable convenient plotting of [band structure diagrams](https://thchr.github.io/Brillouin.jl/stable/kpaths/#Band-structure)).

docs/src/kpaths.md

@@ -9,14 +9,13 @@ sgnum = 227
 Rs = [[1,0,0], [0,1,0], [0,0,1]] # conventional direct basis
 kp = irrfbz_path(sgnum, Rs)
 ```
-The path data is sourced from the [HPKOT paper](http://dx.doi.org/10.1016/j.commatsci.2016.10.015) (or, equivalently, the [SeeK-path](https://github.com/giovannipizzi/seekpath) Python package).
+The path data is sourced from the [HPKOT paper](http://dx.doi.org/10.1016/j.commatsci.2016.10.015) (or, equivalently, the [SeeK-path](https://github.com/giovannipizzi/seekpath) Python package). 
 
-The resulting `KPath` structure initially gives the **k**-point coordinates in the basis of the *primitive* reciprocal basis. To convert to a Cartesian basis, we use [`cartesianize!`](@ref):
+The coordinates of the path are given with respect to the *primitive* reciprocal basis (here, `[[-2π,2π,2π], [2π,-2π,2π], [2π,2π,-2π]]`). To convert to a Cartesian basis, we can use [`cartesianize`](@ref) or [`cartesianize!`](@ref) (in-place):
 ```@example kpath
-pGs = 2π.*[[-1,1,1], [1,-1,1], [1,1,-1]] # primitive reciprocal basis
-cartesianize!(kp, pGs)
+cartesianize(kp)
 ```
-We can visualize the **k**-path using [PlotlyJS.jl](https://github.com/JuliaPlots/PlotlyJS.jl):
+We can visualize the **k**-path using [PlotlyJS.jl](https://github.com/JuliaPlots/PlotlyJS.jl) (conversion to a Cartesian basis for plotting is automatic):
 ```@example kpath
 using PlotlyJS
 Pᵏ = plot(kp)
@@ -32,17 +31,17 @@ Main.HTMLPlot(Pᶜ⁺ᵏ) # hide
 
 ## Interpolation
 Interpolation of a `KPath` structure can be achieved using either [`interpolate(::KPath, ::Integer)`](@ref) or [`splice(::KPath, ::Integer)`](@ref), returning a `KPathInterpolant`.
-As an example, `interpolate(kp, N)` returns an interpolation with a *target* of `N` interpolation points, distributed as equidistantly as possible:
+As an example, `interpolate(kp, N)` returns an interpolation with a *target* of `N` interpolation points, distributed as equidistantly as possible (with the distance metric evaluated in Cartesian space):
 ```@example kpath
 kpi = interpolate(kp, 100)
 ```
 The returned `KPathInterpolant` implements the `AbstractVector` interface, with iterants returning `SVector{D, Float64}` elements.
 To get a conventional "flat" vector, we can simply call `collect(kpi)`.
 
-Internally, `KPathInterpolant`includes additional structure and information: namely, the high-symmetry points and associated labels along the path as well as a partitioning into connected vs. disconnected path segments.
+Internally, `KPathInterpolant` includes additional structure and information: namely, the high-symmetry points and associated labels along the path as well as a partitioning into connected vs. disconnected path segments.
 
 ## Band structure
-The additional structure in `KPathInterpolation`s enable convenient and clear visualization of band structure diagrams in combination with PlotlyJS.
+The additional structure of `KPathInterpolation` enables convenient and clear visualizations of band structure diagrams in combination with PlotlyJS.
 
 To illustrate this, suppose we are considering a tight-binding problem for an *s*-orbital situated at the 1a Wyckoff position. Such a problem has a single band with dispersion [^1] (assuming a cubic side length ``a = 1``):
 ```math

docs/src/wignerseitz.md

@@ -9,6 +9,7 @@ using Brillouin
 Rs = [[0.0, 0.5, 0.5], [0.5, 0.0, 0.5], [0.5, 0.5, 0.0]]
 cᴿ = wignerseitz(Rs)
 ```
+Note that the coordinates of the Wigner-Seitz cell vertices are referred to the basis `Rs`; to convert to Cartesian space, see [`cartesianize`](@ref) and [`cartesianize!`](@ref) (in-place).
 
 We can plot the generated cells using e.g. [PlotlyJS.jl](https://github.com/JuliaPlots/PlotlyJS.jl) via `plot(cᴿ)` (or, alternatively, via a 3D-capable backend of [AbstractPlotting.jl](https://github.com/JuliaPlots/AbstractPlotting.jl)):
 ```@example wignerseitz-cF

src/Brillouin.jl

@@ -12,6 +12,37 @@ const SHOWDIGITS = 6
 # is just a hardcopy of output from Crystalline's equivalent function
 include("CrystallineBravaisVendor.jl")
 # ---------------------------------------------------------------------------------------- #
+@enum BasisEnum CARTESIAN LATTICE
+# ---------------------------------------------------------------------------------------- #
+function latticize! end
+latticize(v::AVec{<:Real}, basismatrix::AbstractMatrix{<:Real}) = basismatrix\v
+latticize(v::AVec{<:Real}, basis::AVec{<:AVec{<:Real}}) = latticize(v, hcat(basis...))'
+latticize(x) = (x.basisenum[] === CARTESIAN ? latticize!(deepcopy(x)) : deepcopy(x))
+
+function cartesianize! end
+cartesianize(v::AVec{<:Real}, basis::AVec{<:AVec{<:Real}}) = v'basis
+cartesianize(x) = (x.basisenum[] === LATTICE ? cartesianize!(deepcopy(x)) : deepcopy(x))
+
+"""
+    basis(x::Union{KPath, KPathInterpolant, Cell})
+
+Return the (reciprocal or direct) lattice basis associated with `x`, in Cartesian
+coordinates.
+
+Methods in Brillouin will by default return points in the lattice basis, i.e., points are
+referred to `basis(x)`. This corresponds to the setting `x.basisenum[] == LATTICE`.
+Coordinates may instead be referred to a Cartesian basis, corresponding to
+`x.basisenum[] == CARTESIAN` by using [`cartesianize`](@ref). The result of `basis(x)`,
+however, is invariant to this and always refers to the lattice basis in Cartesian
+coordinates.
+"""
+basis(x) = x.basis
+
+export latticize!, latticize,
+    cartesianize!, cartesianize,
+    basis
+# ---------------------------------------------------------------------------------------- #
+
 # MAIN FUNCTIONALITY
 
 include("KPaths.jl")

src/CrystallineBravaisVendor.jl

@@ -1,7 +1,12 @@
 module CrystallineBravaisVendor
+
+using ..Brillouin: BasisLike
+using StaticArrays
+
+# ---------------------------------------------------------------------------------------- #
 # this is just a raw output of `bravaistype` from Crystalline, to avoid depending on all
 # of Crystalline just for `bravaistype`; basically just hardcoded "memoization"; note that
-# we do not normalize "oA" to "oC" (i.e., this is `bravaistype(_, _, normalize=false)`)
+# we do *not* normalize "oA" to "oC" (i.e., this is `bravaistype(_, _, normalize=false)`)
 
 @inline function boundscheck_sgnum(sgnum::Integer, D::Integer)
     if D == 3
@@ -56,4 +61,57 @@ function bravaistype(sgnum::Integer, Dᵛ::Val{D}=Val(3)) where D
 end
 bravaistype(sgnum::Integer, D::Integer) = bravaistype(sgnum, Val(D))
 
+# ---------------------------------------------------------------------------------------- #
+# vendor the tools from Crystalline needed to convert back and forth between
+# conventional/reciprocal lattice bases
+
+# ... relies on `bravaistype` *not* normalizing "oA" to "oC"
+centering(sgnum::Integer, Dᵛ::Val{D}) where D = last(bravaistype(sgnum, Dᵛ))
+
+const PRIMITIVE_BASIS_MATRICES = (
+    # 1D
+    Base.ImmutableDict('p'=>SMatrix{1,1,Float64}(1)),                # primitive
+    # 2D
+    Base.ImmutableDict('p'=>SMatrix{2,2,Float64}([1 0; 0 1]),        # primitive/simple
+                  'c'=>SMatrix{2,2,Float64}([1 1; -1 1]./2)),   # centered      
+    # 3D
+    Base.ImmutableDict(
+        'P'=>SMatrix{3,3,Float64}([1 0 0; 0 1 0; 0 0 1]),       # primitive/simple
+        'F'=>SMatrix{3,3,Float64}([0 1 1; 1 0 1; 1 1 0]./2),    # face-centered
+        'I'=>SMatrix{3,3,Float64}([-1 1 1; 1 -1 1; 1 1 -1]./2), # body-centered
+        'R'=>SMatrix{3,3,Float64}([2 -1 -1; 1 1 -2; 1 1 1]./3), # rhombohedrally-centered
+        'A'=>SMatrix{3,3,Float64}([2 0 0; 0 1 -1; 0 1 1]./2),   # base-centered (along x)
+        'C'=>SMatrix{3,3,Float64}([1 1 0; -1 1 0; 0 0 2]./2))   # base-centered (along z)
+    )
+@inline function primitivebasismatrix(cntr::Char, ::Val{D}=Val(3)) where D
+    D ∉ 1:3 && throw(DomainError(D, "dimension must be 1, 2, or 3"))
+    return PRIMITIVE_BASIS_MATRICES[D][cntr]
+end
+
+function reciprocalbasis(Rs::BasisLike{D}) where D
+    Rm = hcat(Rs...)
+    Gm = 2π.*inv(transpose(Rm))
+    return SVector{D}(ntuple(i->SVector{D,Float64}(Gm[:,i]), Val(D)))
+end
+
+function transform_reciprocal(Gs::BasisLike{D}, P::AbstractMatrix{<:Real}) where D
+    Gm′ = hcat(Gs...)/P'
+    return SVector{D}(ntuple(i->SVector{D}(Gm′[:,i]), Val(D)))
+end
+
+for (f, op_P) in ((:primitivize_reciprocal, :identity), (:conventionalize_reciprocal, :inv))
+    @eval function $f(basis::BasisLike{D}, sgnum::Integer) where D
+        cntr = centering(sgnum, Val(D))
+        return $f(basis, cntr)
+    end
+    @eval function $f(Gs::BasisLike{D}, cntr::Char) where D
+        if cntr == 'P' || cntr == 'p' # the conventional and primitive bases coincide
+            return Gs
+        else         
+            P = primitivebasismatrix(cntr, Val(D))        
+            return transform_reciprocal(Gs, $op_P(P))
+        end
+    end
+end
+
 end # module