Finalize overhaul over group loading and add magnetic space group functionality

- Bump Crystalline to v0.5.0

Among many updates, this includes:
- 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

Author: thchr

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

@@ -1,34 +1,8 @@
-#module ParseMSG
-
 using Crystalline
+using StaticArrays
 # --- definitions ---
 # to be transferred to Crystalline's definitions
-using Crystalline: AbstractGroup
-using StaticArrays
 
-struct MSymOperation{D}
-    op :: SymOperation{D}
-    tr :: Bool
-end
-Base.show(io::IO, mop::MSymOperation) = (show(io, mop.op); print(io, mop.tr ? "′" : ""))
-function Crystalline.compose(
-        mop₁::MSymOperation{D},
-        mop₂::MSymOperation{D},
-        modτ::Bool=true) where D
-    return MSymOperation{D}(compose(mop₁.op, mop₂.op, modτ), xor(mop₁.tr, mop₂.tr))
-end
-Base.:*(mop₁::MSymOperation{D}, mop₂::MSymOperation{D}) where D = compose(mop₁, mop₂)
-function Base.isapprox(mop₁::MSymOperation{D}, mop₂::MSymOperation{D}, vs...; kws...) where D
-    isapprox(mop₁.op, mop₂.op, vs...; kws...) && mop₁.tr == mop₂.tr
-end
-struct MSpaceGroup{D} <: AbstractGroup{D}
-    # Note: we use BNS numbers, not  OG numbers (see Sec. 4, Acta Cryst. A78, 99 (2022))
-    # BNS number: `num[1]`: F space-group associated number
-    #             `num[2]`: crystal-system sequential number
-    num :: Tuple{Int, Int}
-    operations :: Vector{MSymOperation{D}}
-end
-Crystalline.label(msg::MSpaceGroup) = MSG_BNS_LABELs_D[msg.num]::String
 # --- parsing ---
 #cd(@__DIR__)
 io = open((@__DIR__)*"/../data/operations/msgs/iso-magnetic_table_bns.txt")
@@ -69,6 +43,7 @@ function read_group!(io)
     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
@@ -109,6 +84,4 @@ for i in 1:1651
     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
-
-#end # module ParseMSG
+end

src/Crystalline.jl

@@ -38,11 +38,12 @@ export smith, Smith # export, so that loading Crystalline effectively also provi
 
 @reexport using Bravais
 import Bravais: primitivize, conventionalize, cartesianize, transform, centering
-using Bravais: stack, all_centeringtranslations, centeringtranslation
+using Bravais: stack, all_centeringtranslations, centeringtranslation,
+               centering_volume_fraction
 
 # included files and exports
 include("constants.jl")
-export MAX_SGNUM, MAX_SUBGNUM, ENANTIOMORPHIC_PAIRS
+export MAX_SGNUM, MAX_SUBGNUM, MAX_MSGNUM, MAX_MSUBGNUM, ENANTIOMORPHIC_PAIRS
 
 include("utils.jl") # useful utility methods (seldom needs exporting)
 export splice_kvpath, interpolate_kvpath
@@ -71,17 +72,35 @@ export SymOperation,                        # types
        irreplabels, klabels,                # ::BandRep & ::BandRepSet 
        isspinful
 
-include("show.jl") # custom printing for structs defined in src/types.jl
-
 include("notation.jl")
 export schoenflies, iuc, centering, seitz, mulliken
 
+include("subperiodic.jl")
+export SubperiodicGroup
+
+include("magnetic/notation-data.jl")
+include("magnetic/types.jl")
+export MSymOperation, MSpaceGroup
+
+include("tables/rotation_translation.jl")
+include("tables/pointgroup.jl")
+include("tables/spacegroup.jl")
+include("tables/subperiodicgroup.jl")
+include("tables/mspacegroup.jl")
+export pointgroup, spacegroup, subperiodicgroup, mspacegroup
+
+include("group-assembly/assemble_pointgroup.jl")
+include("group-assembly/assemble_spacegroup.jl")
+include("group-assembly/assemble_subperiodicgroup.jl")
+include("group-assembly/assemble_mspacegroup.jl")
+
+include("show.jl") # custom printing for structs defined in src/types.jl
+
 include("orders.jl")
 
 include("symops.jl") # symmetry operations for space, plane, and line groups
-export @S_str, spacegroup, compose,
+export @S_str, compose,
        issymmorph, littlegroup, orbit,
-       pointgroup,
        reduce_ops,
        issubgroup, isnormal,
        generate, generators
@@ -99,8 +118,7 @@ include("symeigs2irrep.jl") # find irrep multiplicities from symmetry eigenvalue
 export find_representation
 
 include("pointgroup.jl") # symmetry operations for crystallographic point groups
-export pointgroup, pgirreps,
-       PG_IUCs, find_isomorphic_parent_pointgroup
+export pgirreps, PG_IUCs, find_isomorphic_parent_pointgroup
 
 include("irreps_reality.jl")
 export calc_reality, realify
@@ -133,9 +151,6 @@ export bandreps, matrix, classification, nontrivial_factors, basisdim
 include("deprecations.jl")
 export get_littlegroups, get_lgirreps, get_pgirreps, WyckPos, kvec, wyck, kstar
 
-include("subperiodic.jl")
-export SubperiodicGroup, subperiodicgroup
-
 include("grouprelations/grouprelations.jl")
 export maximal_subgroups, minimal_supergroups
 

src/bandrep.jl

@@ -253,7 +253,7 @@ end
 function matching_lgirreps(BRS::BandRepSet)
     lgirsd = lgirreps(num(BRS), Val(3))
     # create "physical/real" irreps if `BRS` assumes time-reversal symmetry
-    if BRS.timeinvar 
+    if BRS.timereversal 
         for (klab, lgirs) in lgirsd
             lgirsd[klab] = realify(lgirs)
         end

src/conjugacy.jl

@@ -34,11 +34,7 @@ function classes(
     cntr_ops = if cntr === nothing
         nothing
     else
-        if cntr == 'P' || cntr == 'p'
-            SymOperation{D}[]
-        else
-            SymOperation{D}.(all_centeringtranslations(cntr, Val(D)))
-        end
+        SymOperation{D}.(all_centeringtranslations(cntr, Val(D)))
     end
 
     conj_classes = Vector{Vector{SymOperation{D}}}()

src/constants.jl

@@ -8,7 +8,7 @@ const NULL_ATOL    = 1e-11        # absolute tolerance for nullspace
 Return the number of distinct space group types across dimensions 1, 2, and 3 (indexable by
 dimensionality).
 """
-const MAX_SGNUM    = (2, 17, 230) # number of space groups in dimensions 1, 2, and 3
+const MAX_SGNUM    = (2, 17, 230)
 
 @doc raw"""
     MAX_SUBGNUM :: ImmutableDict
@@ -23,6 +23,21 @@ subperiodic group of dimensionality `D` and periodicity `P`:
 """
 const MAX_SUBGNUM = ImmutableDict((3,2)=>80, (3,1)=>75, (2,1)=>7)
 
+@doc raw"""
+    MAX_MSGNUM :: Tuple{Int,Int,Int}
+
+Analogous to `MAX_SUBGNUM`, but for the number of magnetic space groups.
+"""
+const MAX_MSGNUM  = (7, 80, 1651) # Figure 1.2.1 of Litvin's book
+
+
+@doc raw"""
+    MAX_MSUBGNUM :: Tuple{Int,Int,Int}
+
+Analogous to `MAX_SUBGNUM`, but for the number of magnetic subperiodic groups.
+"""
+const MAX_MSUBGNUM = ImmutableDict((3,2)=>528, (3,1)=>394, (2,1)=>31) # Litvin, Fig. 1.2.1
+
 # --- VECTORS ---
 # arbitrary test vector for e.g. evaluating KVecs lines/planes/volumes; 3D by default
 const TEST_αβγ =  [0.123123123, 0.456456456, 0.789789789]

src/group-assembly/assemble_mspacegroup.jl

@@ -0,0 +1,44 @@
+function mspacegroup(BNS₁::Integer, BNS₂::Integer, Dᵛ::Val{D}=Val(3)) where D
+    D == 3 || throw(DomainError("only 3D magnetic space groups are supported"))
+    msgnum = (BNS₁, BNS₂)
+    codes = MSG_CODES_D[msgnum]
+
+    label = MSG_BNS_LABELs_D[msgnum]
+    cntr = first(label)
+    Ncntr = centering_volume_fraction(cntr, Dᵛ)
+    Nop = (length(codes)+1) # number of operations ÷ by equiv. centering translations
+    operations = Vector{MSymOperation{D}}(undef, Nop * Ncntr)
+
+    operations[1] = MSymOperation{D}(one(SymOperation{D}), false)
+    for (n, code) in enumerate(codes)
+        op = SymOperation{D}(ROTATIONS_3D[code[1]], TRANSLATIONS_3D[code[2]])
+        operations[n+1] = MSymOperation{D}(op, code[3])
+    end
+
+    # add the centering-translation related operations (not included in `codes`)
+    if Ncntr > 1
+        cntr_translations = all_centeringtranslations(cntr, Dᵛ)
+        for (i, t) in enumerate(cntr_translations)
+            for n in 1:Nop
+                mop = operations[n]
+                t′ = reduce_translation_to_unitrange(translation(mop) + t)
+                op′ = SymOperation{D}(mop.op.rotation, t′)
+                mop′ = MSymOperation{D}(op′, mop.tr)
+                operations[n+i*Nop] = mop′
+            end
+        end
+    end
+
+    return MSpaceGroup{D}(msgnum, operations)
+end
+
+function mspacegroup(OG₃::Integer, ::Val{D}=Val(3)) where D
+    D == 3 || throw(DomainError("only 3-dimensional magnetic space groups are supported"))
+    OG₃ ∈ 1:MAX_MSGNUM[D] || _throw_invalid_msgnum(OG₃, D)
+
+    BNS₁, BNS₂ = MSG_OG₃2BNS_NUMs_V[OG₃]
+    return mspacegroup(BNS₁, BNS₂, Val(3))
+end
+@noinline function _throw_invalid_msgnum(OG₃, D)
+    throw(DomainError(OG₃, "the sequential magnetic space group number must be between 1 and $(MAX_MSGNUM[D]) in dimension $D"))
+end

src/group-assembly/assemble_pointgroup.jl

@@ -0,0 +1,58 @@
+"""
+    pointgroup(iuclab::String, ::Union{Val{D}, Integer}=Val(3))  -->  PointGroup{D}
+
+Return the symmetry operations associated with the point group identified with label
+`iuclab` in dimension `D` as a `PointGroup{D}`.
+"""
+function pointgroup(iuclab::AbstractString, Dᵛ::Val{D}=Val(3)) where D
+    @boundscheck _check_valid_pointgroup_label(iuclab, D)
+    pgnum = pointgroup_iuc2num(iuclab, D) # this is not generally a particularly well-established numbering
+    return _pointgroup(iuclab, pgnum, Dᵛ)
+end
+@inline pointgroup(iuclab::String, D::Integer) = pointgroup(iuclab, Val(D))
+
+"""
+    pointgroup(pgnum::Integer, ::Union{Val{D}, Integer}=Val(3), setting::Integer=1)
+                                                                      -->  PointGroup{D}
+
+Return the symmetry operations associated with the point group identfied with canonical
+number `pgnum` in dimension `D` as a `PointGroup{D}`. The connection between a point group's
+numbering and its IUC label is enumerated in `Crystalline.PG_NUM2IUC[D]` and
+`Crystalline.IUC2NUM[D]`.
+
+Certain point groups feature in multiple setting variants: e.g., IUC labels 321 and 312 both
+correspond to `pgnum = 18` and correspond to the same group structure expressed in two
+different settings. The `setting` argument allows choosing between these setting variations.
+"""
+function pointgroup(pgnum::Integer, Dᵛ::Val{D}=Val(3), setting::Integer=1) where D
+    iuclab = pointgroup_num2iuc(pgnum, Dᵛ, setting) # also checks validity of `(pgnum, D)`
+    return _pointgroup(iuclab, pgnum, Dᵛ)
+end
+@inline function pointgroup(pgnum::Integer, D::Integer, setting::Integer=1)
+    return pointgroup(pgnum, Val(D), setting)
+end
+
+function _pointgroup(iuclab::String, pgnum::Integer, Dᵛ::Val{D}) where D
+    codes = D == 3 ? PG_CODES_3D_D[iuclab] : 
+               D == 2 ? PG_CODES_2D_D[iuclab] :
+                        PG_CODES_1D_D[iuclab]
+
+    Nop = (length(codes)+1) # number of operations
+    operations = Vector{SymOperation{D}}(undef, Nop)
+    operations[1] = one(SymOperation{D})
+    for (n, code) in enumerate(codes)
+        op = SymOperation{D}(get_indexed_rotation(code[1], Dᵛ), 
+                             zero(SVector{D,Float64}))
+        operations[n+1] = op
+    end
+
+    return PointGroup{D}(pgnum, iuclab, operations)
+end
+
+function _check_valid_pointgroup_label(iuclab::AbstractString, D)
+    D ∉ (1,2,3) && _throw_invalid_dim(D)
+    if iuclab ∉ PG_IUCs[D]
+        throw(DomainError(iuclab, 
+                "iuc label not found in database (see possible labels in `PG_IUCs[D]`)"))
+    end
+end