Features/spin bivector

CHANGELOG.md

@@ -9,6 +9,7 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ### Added
   - `eachvertex` iterator for the vertices of the parallelepiped representation of a unit cell.
+  - `SpinBivector` type representing spin direction in a dimension-agnostic way.
 
 ### Changed
   - The Types section of the documentation has been split up into separate sections for lattice

docs/src/api/data.md

@@ -37,6 +37,7 @@ Electrum.weight
 ## Spin states
 ```@docs
 Electrum.Multiplicity
+Electrum.SpinBivector
 ```
 
 ## Energies and occupancies

src/Electrum.jl

@@ -136,7 +136,7 @@ export energy, occupancy, energies, occupancies, min_energy, max_energy, min_occ
     max_occupancy, fermi
 # Spin data: multiplicity and direction
 include("data/spin.jl")
-export Multiplicity
+export Multiplicity, SpinBivector
 # Real and reciprocal space data grids
 include("data/grids.jl")
 export DataGrid, RealDataGrid, ReciprocalDataGrid

src/data/spin.jl

@@ -39,3 +39,113 @@
 
 Base.UnitRange(s::Multiplicity) = first(s):last(s)
 Base.show(io::IO, s::Multiplicity) = print(io, typeof(s), "()")
+
+#---Spin bivectors---------------------------------------------------------------------------------#
+"""
+    Electrum._is_skew_symmetric(m::AbstractMatrix)
+
+Returns `true` if a matrix is skew-symmetric (equivalently, antisymmetric).
+"""
+function _is_skew_symmetric(m::AbstractMatrix)
+    size(m, 1) == size(m, 2) || return false
+    for a in axes(m, 1)
+        for b in (a+1):lastindex(m, 1)
+            m[a,b] == -m[b,a] || return false
+        end
+    end
+    return true
+end
+
+"""
+    SpinBivector{D,T}
+
+A skew-symmetric matrix representing a spin [bivector](https://en.wikipedia.org/wiki/Bivector) in 
+`D` dimensions. No automatic normalization has been implemented for this type.
+
+# Why the spin axis is a bivector
+
+The identification of a spin direction with a vector is only possible in 3D. This is because the
+axis of rotation is defined by the [Hodge dual](https://en.wikipedia.org/wiki/Hodge_star_operator)
+of the rotation plane. The Hodge dual is a linear map relating `k`-dimensional objects in 
+`D`-dimensional space to `D-k`-dimensional objects in the same space. 
+
+A bivector is a 2-dimensional object, and in three dimensions, its Hodge dual is a vector. This is
+not the case in any other number of dimensions. To allow for generality in describing spin, this
+representation is used instead.
+"""
+struct SpinBivector{D,T} <: StaticMatrix{D,D,T}
+    data::SVector{D,SVector{D,T}}
+    function SpinBivector{D,T}(m::StaticMatrix) where {D,T}
+        @assert _is_skew_symmetric(m) "Input matrix is not skew-symmetric."
+        return new(SVector{D,T}.(eachcol(m)))
+    end
+end
+
+# Needed to resolve method ambiguities
+SpinBivector{D,T}(::StaticArray) where {D,T} = error("Input must be a skew-symmetric matrix.")
+SpinBivector{D}(::StaticArray) where D = error("Input must be a skew-symmetric matrix.")
+SpinBivector(::StaticArray) = error("Input must be a skew-symmetric matrix.")
+
+SpinBivector{D}(m::StaticMatrix) where D = SpinBivector{D,eltype(m)}(m)
+SpinBivector(m::StaticMatrix{D,D}) where D = SpinBivector{D,eltype(m)}(m)
+
+SpinBivector{D,T}(m::AbstractMatrix) where {D,T} = SpinBivector{D,T}(SMatrix{D,D}(m))
+SpinBivector{D}(m::AbstractMatrix) where D = SpinBivector{D,eltype(m)}(SMatrix{D,D}(m))
+
+function Base.getproperty(b::SpinBivector, s::Symbol)
+    s === :matrix && return hcat(getfield(b, :data)...)
+    return getfield(b, s)
+end
+
+Base.propertynames(::SpinBivector; private = false) = private ? (:data, :matrix) : (:matrix,)
+
+# Really only for resolving method ambiguities
+Base.getindex(b::SpinBivector, i::Int) = getindex(b.matrix, i)
+Base.getindex(b::SpinBivector, i::Int...) = getindex(b.matrix, i...)
+
+Base.:(==)(a::SpinBivector, b::SpinBivector) = (a.matrix == b.matrix)
+DataSpace(::Type{<:SpinBivector{D}}) where D = ByRealSpace{D}()
+# Required for StaticArray subtypes
+Tuple(b::SpinBivector) = Tuple(b.matrix)
+
+_vector_length_mismatch() = throw(DimensionMismatch("Vectors must be the same length."))
+
+"""
+    Electrum._wedge_matrix(u::AbstractVector, v::AbstractVector) -> AbstractMatrix
+
+Returns a bivector, the wedge product of vectors `u` and `v`, represented as a skew-symmetric
+matrix. This function will throw a `DimensionMismatch` if `u` and `v` have different lengths.
+"""
+function _wedge_matrix(u::AbstractVector, v::AbstractVector)
+    length(u) === length(v) || _vector_length_mismatch()
+    return u*v' - v*u'
+end
+
+# Try to produce a static matrix whenever possible
+_wedge_matrix(u::StaticVector{D}, v::StaticVector{D}) where D = u*v' - v*u'
+_wedge_matrix(u::StaticVector{D}, v::AbstractVector) where D = _wedge_matrix(u, SVector{D}(v))
+_wedge_matrix(u::AbstractVector, v::StaticVector{D}) where D = _wedge_matrix(SVector{D}(u), v)
+# Needed to resolve method ambiguities
+_wedge_matrix(::StaticVector, ::StaticVector) = _vector_length_mismatch()
+
+# Constructors for taking wedge products implicitly
+"""
+    SpinBivector(u::StaticVector{D}, v::StaticVector{D})
+    SpinBivector{D}(u::AbstractVector, v::AbstractVector)
+    SpinBivector{D,T}(u::AbstractVector, v::AbstractVector)
+
+Constructs a spin bivector from the wedge products of vectors `u` and `v`, which define a bivector.
+
+The first constructor may be used if only one argument is a `StaticVector{D}` (the other will 
+automatically be converted to the correct size).
+"""
+function SpinBivector(u::AbstractVector, v::AbstractVector)
+    w = _wedge_matrix(u,v)
+    w isa StaticArray && return SpinBivector(w)
+    error("Vector lengths cannot be inferred from input types. Use SpinBivector{D} instead.")
+end
+
+(T::Type{<:SpinBivector{D}})(u::AbstractVector, v::AbstractVector) where D = T(_wedge_matrix(u,v))
+
+# TODO: mathematical operations should return reasonable types
+Base.:-(s::SpinBivector) = typeof(s)(-s.matrix)

test/spin.jl

@@ -1,4 +1,4 @@
-@testset "Spin" begin
+@testset "Multiplicity" begin
     @test all(size(Multiplicity(n)) == (n,) for n in 1:100)
     @test all(last(Multiplicity(n)) == (n - 1)//2 for n in 1:100)
     @test all(Multiplicity(n)[2] == -(n - 1)//2 + 1 for n in 2:100)
@@ -21,3 +21,43 @@
     # Text representation
     @test eval(Meta.parse(repr(Multiplicity(3)))) === Multiplicity(3)
 end
+
+@testset "Spin bivectors" begin
+    x = [1, 0, 0]
+    y = [0, 1, 0]
+    (sx, sy) = SVector{3}.((x, y))
+    z_matrix = @SMatrix [0 1 0; -1 0 0; 0 0 0]
+    s_z = SpinBivector(z_matrix)
+    @test Electrum._is_skew_symmetric(z_matrix)
+    @test !Electrum._is_skew_symmetric([1 2 3; 4 5 6; 7 8 9])
+    @test_throws AssertionError SpinBivector{3}([1 2 3; 4 5 6; 7 8 9])
+    @test SpinBivector(sx, y) === s_z
+    @test SpinBivector(y, sx) === -s_z
+    @test SpinBivector{3}([0 1 0; -1 0 0; 0 0 0]) === s_z
+    @test SpinBivector{3,Int}([0 1 0; -1 0 0; 0 0 0]) === s_z
+    # Constructor with wedge product
+    @test SpinBivector(sx, sy) === s_z
+    @test SpinBivector{3}(x, y) === s_z
+    @test SpinBivector{3,Float64}(x, y) == s_z
+    # Wedge product kernel for AbstractVector
+    @test Electrum._wedge_matrix(sx, sy) == Electrum._wedge_matrix(x, y)
+    # Properties
+    @test :data in propertynames(s_z, private = true)
+    @test only(propertynames(s_z)) == :matrix
+    # Indexing
+    @test s_z[2,1] == -1
+    @test s_z[1,2] == 1
+    @test s_z[2] == -1
+    @test s_z[4] == 1
+    # Bad constructor inputs
+    @test_throws ErrorException SpinBivector(x, y)
+    @test_throws ErrorException SpinBivector(@SVector [1, 2, 3])
+    @test_throws ErrorException SpinBivector{3}(@SVector [1, 2, 3])
+    @test_throws ErrorException SpinBivector{3,Int}(@SVector [1, 2, 3])
+    @test_throws ErrorException SpinBivector(SArray{Tuple{2,2,2}}(1, 2, 3, 4, 5, 6, 7, 8))
+    @test_throws DimensionMismatch SpinBivector(SVector{2}(1, 0), SVector{3}(0, 1, 0))
+    @test_throws DimensionMismatch SpinBivector{3}([1, 0, 0], [0, 1])
+    # Traits
+    @test Electrum.DataSpace(SpinBivector{3}) === Electrum.ByRealSpace{3}()
+    @test Electrum.DataSpace(s_z) === Electrum.ByRealSpace{3}()
+end