Add magnetic space group tooling and improve group operation loading (f…

…ix #39) (#54)

Bravais.jl: 
- bumps Bravais version to v0.1.9
- This makes the canonicalization of subperiodic group centering symbols more meaningful; now always going to the associated symbol in `D` dimensions rather than `P`
- this fixes a bug for  `all_centeringtranslations` and simplifies the treatment in `primitivebasismatrix`
- also introduces and exports `centering_volume_fraction` which gives the ratio between the conventional and primitive unit cell volumes

Crystalline:
- Bump Crystalline to v0.5.0
- Add magnetic space group functionality (`mspacegroup` and `MSymOperation`)
- Now, load all group operation data by reference to a fixed set of operators (in `tables/rotation_translation.jl` via a set of codes (stored in `tables/spacegroup.jl` and similar for point/subperiodic/magnetic groups): this is much faster computationally as it avoids having to load anything from disk and solves #39
- Introduce a new `AbstractOperation` abstract type to fold both `SymOperation` and `MSymOperation` (for magnetic operations) under one umbrella
- All assembly of groups is now in `group-assembly/...`
- Various minor fixes/code-improvement to subperiodic groups
- Add an element type parameter to `AbstractGroup`
- Fix various minor things

Bravais/Project.toml

@@ -1,7 +1,7 @@
 name = "Bravais"
 uuid = "ada6cbde-b013-4edf-aa94-f6abe8bd6e6b"
 authors = ["Thomas Christensen <[email protected]>"]
-version = "0.1.8"
+version = "0.1.9"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

Bravais/src/Bravais.jl

@@ -23,6 +23,7 @@ export crystal,
     bravaistype,
     centering,
     primitivebasismatrix,
+    centering_volume_fraction,
     directbasis,
     reciprocalbasis,
     AbstractBasis,

Bravais/src/transform.jl

@@ -69,10 +69,10 @@ function canonicalize_centering(cntr, ::Val{D}, ::Val{P}) where {D,P}
     if D == P # space/plane/line groups
         return cntr
     elseif D == 3 && P == 2 # layer groups
-        return cntr == '𝑝' ? 'p' : 
-               cntr == '𝑐' ? 'c' : error(DomainError(cntr, "invalid layer group centering"))
+        return cntr == '𝑝' ? 'P' : 
+               cntr == '𝑐' ? 'C' : error(DomainError(cntr, "invalid layer group centering"))
     elseif D == 3 && P == 1 # rod groups
-        return cntr == '𝓅' ? 'p' : error(DomainError(cntr, "invalid rod group centering"))
+        return cntr == '𝓅' ? 'P' : error(DomainError(cntr, "invalid rod group centering"))
     elseif D == 2 && P == 1 # frieze groups
         return cntr == '𝓅' ? 'p' : error(DomainError(cntr, "invalid frieze group centering"))
     else
@@ -161,38 +161,18 @@ frequently and more ambiguously, as the crystallographic primitive setting.
            Thus, Bravais.jl and [Spglib.jl](https://github.com/singularitti/Spglib.jl)
            transform to identical primitive settings and are hence mutually compatible.
 """
-@inline function primitivebasismatrix(cntr::Char, ::Val{D}=Val(3), ::Val{D}=Val(D)) where D
-    # space groups
+@inline function primitivebasismatrix(cntr::Char,
+                                      Dᵛ::Val{D}=Val(3), Pᵛ::Val{P}=Val(D)) where {D,P}
     D ∉ 1:3 && _throw_invalid_dim(D)
-    return PRIMITIVE_BASIS_MATRICES[D][cntr]
-end
-@inline function primitivebasismatrix(cntr::Char, Dᵛ::Val{3}, Pᵛ::Val{2})
-    # layer and rod groups
+    P ∉ 1:D && throw(DomainError((D,P), "invalid combination of dimensionality D and periodicity P"))
     cntr = canonicalize_centering(cntr, Dᵛ, Pᵛ)
-    P²ᴰ = PRIMITIVE_BASIS_MATRICES[2][cntr]
-    return @SMatrix [P²ᴰ[1,1] P²ᴰ[2,1] 0.0; P²ᴰ[1,2] P²ᴰ[2,2] 0.0; 0.0 0.0 1.0]
-end
-@inline function primitivebasismatrix(cntr::Char, Dᵛ::Val{3}, Pᵛ::Val{1})
-    # rod groups
-    cntr = canonicalize_centering(cntr, Dᵛ, Pᵛ)
-    P²ᴰ = PRIMITIVE_BASIS_MATRICES[1][cntr]
-    return @SMatrix [P²ᴰ[1,1] 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 1.0]
-end
-@inline function primitivebasismatrix(cntr::Char, Dᵛ::Val{2}, Pᵛ::Val{1})
-    # frieze groups
-    cntr = canonicalize_centering(cntr, Dᵛ, Pᵛ)
-    P¹ᴰ = PRIMITIVE_BASIS_MATRICES[1][cntr]
-    return @SMatrix [P¹ᴰ[1,1] 0.0; 0.0 1.0]
-end
-function primitivebasismatrix(cntr::Char, ::Val{D}, ::Val{P}) where {D,P}
-    # fall-back error
-    throw(DomainError((D,P), "invalid combination of dimensionality D and periodicity P"))
+    return PRIMITIVE_BASIS_MATRICES[D][cntr]
 end
 
 @inline function centeringtranslation(cntr::Char,
                                       Dᵛ::Val{D}=Val(3), Pᵛ::Val{P}=Val(D)) where {D,P}
     D ∉ 1:3 && _throw_invalid_dim(D)
-    P ∉ 1:D && _throw_invalid_dim(P)
+    P ∉ 1:D && throw(DomainError((D,P), "invalid combination of dimensionality D and periodicity P"))
     cntr = canonicalize_centering(cntr, Dᵛ, Pᵛ)
     if D == 3
         if     cntr == 'P'; return zeros(SVector{3})
@@ -212,28 +192,53 @@ end
         if     cntr == 'p'; return zeros(SVector{1})
         else;               _throw_invalid_cntr(cntr, 1)
         end
-    else 
-        _throw_invalid_dim(D)
     end
+    error("unreachable reached")
 end
 
 function all_centeringtranslations(cntr::Char,
                                    Dᵛ::Val{D}=Val(3), Pᵛ::Val{P}=Val(D)) where {D,P}
     D ∉ 1:3 && _throw_invalid_dim(D)
     P ∉ 1:D && _throw_invalid_dim(P)
     cntr = canonicalize_centering(cntr, Dᵛ, Pᵛ)
-    if D == 3 && cntr == 'F'
+    if cntr == 'P' || cntr == 'p'
+        # primitive cell is equal to conventional cell: 0 extra centers
+        return SVector{D,Float64}[]
+    elseif D == 3 && cntr == 'F'
         # primitive cell has 1/4th the volume of conventional cell: 3 extra centers
-        return [SVector((1,0,1)./2), SVector((0,1,1)./2), SVector((1,1,0)./2)]
+        return [SVector((0,1,1)./2), SVector((1,0,1)./2), SVector((1,1,0)./2)]
     elseif D == 3 && cntr == 'R'
         # primitive cell has 1/3rd the volume of conventional cell: 2 extra centers
         return [SVector((2,1,1)./3), SVector((1,2,2)./3)]
-    else
+    else # 'I', 'C', 'c', 'A'
         # primitive cell has half the volume of conventional cell: 1 extra center
         return [centeringtranslation(cntr, Dᵛ)]
     end
 end
 
+"""
+    centering_volume_fraction(cntr, Dᵛ, Pᵛ) --> Int
+
+Return the (integer-) ratio between the volumes of the conventional and primitive unit cells
+for a space or subperiodic group with centering `cntr`, embedding dimension `D`, and
+periodicity dimension `P`.
+"""
+@inline function centering_volume_fraction(cntr::Char,
+                                   Dᵛ::Val{D}=Val(3), Pᵛ::Val{P}=Val(D)) where {D,P}
+    D ∉ 1:3 && _throw_invalid_dim(D)
+    P ∉ 1:D && _throw_invalid_dim(P)
+    cntr = canonicalize_centering(cntr, Dᵛ, Pᵛ)
+    if cntr == 'P' || cntr == 'p'
+        return 1
+    elseif D == 3 && cntr == 'F'
+        return 4
+    elseif D == 3 && cntr == 'R'
+        return 3
+    else # 'I', 'C', 'c', 'A'
+        return 2
+    end
+end
+
 # ---------------------------------------------------------------------------------------- #
 
 """

Project.toml

@@ -1,7 +1,7 @@
 name = "Crystalline"
 uuid = "ae5e2be0-a263-11e9-351e-f94dad1eb351"
 authors = ["Thomas Christensen <[email protected]>"]
-version = "0.4.22"
+version = "0.5.0"
 
 [deps]
 Bravais = "ada6cbde-b013-4edf-aa94-f6abe8bd6e6b"
@@ -29,7 +29,7 @@ CrystallineGraphMakieExt = "GraphMakie"
 CrystallinePyPlotExt = "PyPlot"
 
 [compat]
-Bravais = "0.1.8"
+Bravais = "0.1.9"
 Combinatorics = "1.0"
 DelimitedFiles = "1"
 DocStringExtensions = "0.8, 0.9"

build/crawl_and_write_pgs.jl

@@ -180,8 +180,8 @@ function __crawl_and_write_3d_pgirreps()
 
             # ==== save irreps ====
             # we do not save the point group operations anew; they are already stored in 
-            # "data/operations/pgs/..."; note that we explicitly checked the sorting and
-            # equivalence of operations when pgirs was crawled above (cf. flag
+            # "test/data/xyzt-operations/pgs/..."; note that we explicitly checked the 
+            # sorting and # equivalence of operations when pgirs was crawled above (cf. flag
             # `consistency_checks=true`)
             unmangled_pgiuc = Crystalline.unmangle_pgiuclab(pgiuc) # replace '/'s by '_slash_'s
             irreps_file[unmangled_pgiuc*"/matrices"]  = matrices

build/crawl_and_write_subperiodic.jl

@@ -53,16 +53,16 @@ end
 
 # ---------------------------------------------------------------------------------------- #
 ## Use crawling functions & write information to:
-#      `/data/generators/subperiodic/{layer|rod|frieze}/`
-#      `/data/operations/subperiodic/{layer|rod|frieze}/`
+#      `/test/data/xyzt-operations/generators/subperiodic/{layer|rod|frieze}/`
+#      `/test/data/xyzt-operations/subperiodic/{layer|rod|frieze}/`
 for (D, P, sub, maxnum) in [(3,2,"layer",80), (3,1,"rod",75), (2,1,"frieze",7)]
     println(sub)
     for kind in ["operations", "generators"]
         println("  ", kind)
         for num in 1:maxnum
             println("    ", num)
             gens_str = crawl_subperiodic_xyzt(num, D, P, kind)
-            filename = (@__DIR__)*"/../data/$kind/subperiodic/$sub/"*string(num)*".csv"
+            filename = (@__DIR__)*"/../test/data/xyzt-$kind/subperiodic/$sub/"*string(num)*".csv"
             open(filename; write=true, create=true, truncate=true) do io
                 first = true
                 for str in gens_str

build/parse_and_write_magnetic_space_groups.jl

@@ -0,0 +1,87 @@
+using Crystalline
+using StaticArrays
+# --- definitions ---
+# to be transferred to Crystalline's definitions
+
+# --- parsing ---
+#cd(@__DIR__)
+io = open((@__DIR__)*"/../data/operations/msgs/iso-magnetic_table_bns.txt")
+
+# --- read dictionary of operations ---
+function read_ops!(io, ops_d, stop_string)
+    while (str=readline(io); str ≠ stop_string)
+        parts = split(str, isspace; keepempty=false)
+        code, xyz = parts[1], parts[2]
+        op = SymOperation(xyz)
+        ops_d[code] = op
+    end
+    return ops_d
+end
+
+readuntil(io, "Non-hexagonal groups:"); readline(io)
+nonhex_ops_d = Dict{String, SymOperation{3}}()
+read_ops!(io, nonhex_ops_d, "Hexagonal groups:")
+hex_ops_d    = Dict{String, SymOperation{3}}()
+read_ops!(io, hex_ops_d, "----------")
+
+# --- read individual groups ---
+function parse_opstr(str::AbstractString, ops_d)
+    str′ = strip(str,  ['(', ')','\''])
+    code, translation_str = split(str′, '|')::Vector{SubString{String}}
+    translation_parts = split(translation_str, ',')
+    translation_tup = ntuple(Val(3)) do i 
+        Crystalline.parsefraction(translation_parts[i])
+    end
+    translation = SVector{3, Float64}(translation_tup)
+    op = SymOperation{3}(ops_d[code].rotation, translation)
+    tr = last(str) == '\''
+    return MSymOperation(op, tr)
+end
+
+function read_group!(io)
+    # BNS: number & label
+    readuntil(io, "BNS: ")
+    num_str = readuntil(io, " ") # numbering in format `sgnum.cnum`
+    sgnum, cnum = parse.(Int, split(num_str, "."))::Vector{Int}
+    # `sgnum`: F space-group associated number, `cnum`:  crystal-system sequential number
+    bns_label = replace(readuntil(io, " "), '\''=>'′') # BNS label in IUC-style
+
+    # OG: number and label
+    readuntil(io, "OG: ")
+    N₁N₂N₃_str = readuntil(io, " ") # numbering in format `N₁.N₂.N₃`
+    N₁, N₂, N₃ = parse.(Int, split(N₁N₂N₃_str, "."))::Vector{Int}
+    og_label   = replace(readuntil(io, '\n'), '\''=>'′') # OG label in IUC-style
+
+    # operator strings
+    is_hex = crystalsystem(sgnum) ∈ ("hexagonal" , "trigonal")
+    readuntil(io, "Operators")
+    c = read(io, Char)
+    if c == ' '
+        readuntil(io, "(BNS):")
+    elseif c ≠ ':'
+        error("unexpected parsing failure")
+    end
+    ops_str = readuntil(io, "Wyckoff")
+    op_strs = split(ops_str, isspace; keepempty=false)
+
+    # parse operator strings to operations
+    g = Vector{MSymOperation{3}}(undef, length(op_strs))
+    for (i, str) in enumerate(op_strs)
+        g[i] = parse_opstr(str, is_hex ? hex_ops_d : nonhex_ops_d)
+        isdone = str[end] == '\n'
+        isdone && break
+    end
+
+    return MSpaceGroup{3}((sgnum, cnum), g), bns_label, og_label, (N₁, N₂, N₃)
+end
+msgs = MSpaceGroup{3}[]
+MSG_BNS_LABELs_D = Dict{Tuple{Int, Int}, String}()
+MSG_OG_LABELs_D = Dict{Tuple{Int, Int, Int}, String}()
+MSG_BNS2OG_NUMs_D = Dict{Tuple{Int, Int}, Tuple{Int, Int, Int}}()
+for i in 1:1651
+    msg, bns_label, og_label, (N₁, N₂, N₃) = read_group!(io)
+    push!(msgs, msg)
+    MSG_BNS_LABELs_D[msg.num] = bns_label
+    MSG_OG_LABELs_D[(N₁, N₂, N₃)] = og_label
+    MSG_BNS2OG_NUMs_D[msg.num] = (N₁, N₂, N₃)
+end