Support loading magnetic structure from MCIF

docs/src/versions.md

@@ -15,6 +15,10 @@
   independent of the calculator ([Issue
   264](https://github.com/SunnySuite/Sunny.jl/issues/264)). Consequently, color
   ranges in plots may need to be rescaled.
+* [`Crystal`](@ref) can now infer a chemical unit cell from an mCIF file.
+  `System` now supports [`set_dipoles_from_mcif!`](@ref). Through spglib, one
+  can now [`standardize`](@ref) any `Crystal`, with an option to idealize site
+  positions.
 
 ## v0.5.11
 (June 2, 2024)

src/MCIF.jl

@@ -0,0 +1,112 @@
+"""
+    set_dipoles_from_mcif!(sys::System, filename::AbstractString)
+
+Load the magnetic supercell according to an mCIF file. System `sys` must already
+be resized to the correct supercell dimensions.
+"""
+function set_dipoles_from_mcif!(sys::System, filename::AbstractString)
+    cif = CIF.Cif(Path(filename))
+    # For now, assumes there is only one data collection per .cif
+    cif = cif[first(keys(cif))]
+
+    a = parse_cif_float(cif["_cell_length_a"][1])
+    b = parse_cif_float(cif["_cell_length_b"][1])
+    c = parse_cif_float(cif["_cell_length_c"][1])
+    α = parse_cif_float(cif["_cell_angle_alpha"][1])
+    β = parse_cif_float(cif["_cell_angle_beta"][1])
+    γ = parse_cif_float(cif["_cell_angle_gamma"][1])
+    supervecs = sys.crystal.latvecs .* sys.latsize
+    supervecs2 = lattice_vectors(a, b, c, α, β, γ)
+
+    # TODO: Tolerance to permutations (with sign flips) of lattice vectors
+    if !isapprox(supervecs, supervecs2; rtol=sys.crystal.symprec)
+        error("""Invalid supercell dimensions,
+                   System: $supervecs
+                   mCIF:   $supervecs2
+                 """)
+    end
+
+    oneof(fields...) = findfirstval(in(keys(cif)), fields)
+    # The first entry is the IUCR standard field name. If missing, search for
+    # alternate field names that appear in legacy files.
+    mcif_fields = (;
+        magn_operation_xyz=oneof("_space_group_symop_magn_operation.xyz", "_space_group_symop.magn_operation_xyz"),
+        magn_centering_xyz=oneof("_space_group_symop_magn_centering.xyz", "_space_group_symop.magn_centering_xyz"),
+        moment_label=oneof("_atom_site_moment.label", "_atom_site_moment_label"),
+        moment_crystalaxis_x=oneof("_atom_site_moment.crystalaxis_x", "_atom_site_moment_crystalaxis_x"),
+        moment_crystalaxis_y=oneof("_atom_site_moment.crystalaxis_y", "_atom_site_moment_crystalaxis_y"),
+        moment_crystalaxis_z=oneof("_atom_site_moment.crystalaxis_z", "_atom_site_moment_crystalaxis_z"),
+    )
+
+    sym_table = CIF.get_loop(cif, mcif_fields.magn_operation_xyz)
+    magn_operations = MSymOp.(sym_table[:, mcif_fields.magn_operation_xyz])
+    sym_table = CIF.get_loop(cif, mcif_fields.magn_centering_xyz)
+    magn_centerings = MSymOp.(sym_table[:, mcif_fields.magn_centering_xyz])
+
+    dip_table = CIF.get_loop(cif, mcif_fields.moment_label)
+    labels = String.(dip_table[:, mcif_fields.moment_label])
+
+    Mxs = parse_cif_float.(dip_table[:, mcif_fields.moment_crystalaxis_x])
+    Mys = parse_cif_float.(dip_table[:, mcif_fields.moment_crystalaxis_y])
+    Mzs = parse_cif_float.(dip_table[:, mcif_fields.moment_crystalaxis_z])
+    moments = Vec3.(zip(Mxs, Mys, Mzs))
+
+    geom_table = CIF.get_loop(cif, "_atom_site_label")
+    all_labels = String.(geom_table[:, "_atom_site_label"])
+    xs = parse_cif_float.(geom_table[:, "_atom_site_fract_x"])
+    ys = parse_cif_float.(geom_table[:, "_atom_site_fract_y"])
+    zs = parse_cif_float.(geom_table[:, "_atom_site_fract_z"])
+    idxs = indexin(labels, all_labels)
+    positions = Vec3.(zip(xs[idxs], ys[idxs], zs[idxs]))
+
+    set_dipoles_from_mcif_aux!(sys; positions, moments, magn_operations, magn_centerings)
+end
+
+# The keyword parameters follow the mCIF format. As a consequence, all (x,y,z)
+# components are provided in the coordinate system defined by the lattice
+# vectors of the supercell. Similarly, the symmetry operations act in this
+# coordinate system.
+function set_dipoles_from_mcif_aux!(sys; positions, moments, magn_operations, magn_centerings)
+    supervecs = sys.crystal.latvecs .* sys.latsize
+
+    # Use the zero vector as a marker for unvisited sites
+    fill!(sys.dipoles, zero(Vec3))
+
+    for (r, μ) in zip(positions, moments)
+        for s1 in magn_operations, s2 in magn_centerings
+            # Symmetry operation that composes centering an then a
+            # crystallographic operation
+            s = s1 * s2
+
+            # Apply symmetry operation to position, and then map to original
+            # fractional coordinates
+            r_new = s.R * r + s.T
+            r_new = orig_crystal(sys).latvecs \ supervecs * r_new
+
+            # Apply symmetry operation to magnetic moment, and then map to
+            # Cartesian coordinates. Because the moment is a pseudo-vector, it
+            # is invariant to an inversion in space. For this reason, remove the
+            # determinant of s.R. If the parity s.p = ±1 is negative, this
+            # implies time reversal, which reverses the magnetic dipole.
+            μ_new = (s.R * det(s.R) * s.p) * μ
+            μ_new = supervecs * μ_new
+
+            site = try 
+                position_to_site(sys, r_new)
+            catch _
+                rethrow(ErrorException("mCIF position $r_new is missing in the chemical cell"))
+            end
+
+            # Get spin dipole by inverting the `magnetic_moment` transformation
+            dipole = inv(- sys.units.μB * sys.gs[site]) * μ_new
+            s_prev = sys.dipoles[site]
+            set_dipole!(sys, dipole, site)
+            s_prev == zero(Vec3) || s_prev ≈ sys.dipoles[site] || error("Conflicting dipoles at site $site")
+        end
+    end
+
+    unvisited = [site for site in eachsite(sys) if iszero(sys.dipoles[site])]
+    if !isempty(unvisited)
+        error("Missing dipoles for sites $(collect(Tuple.(unvisited)))")
+    end
+end

src/MathBasics.jl

@@ -55,3 +55,22 @@ function tracelesspart(A)
     @assert allequal(size(A))
     return A - tr(A) * I / size(A,1)
 end
+
+# https://github.com/JuliaLang/julia/issues/44996
+function findfirstval(f, a)
+    i = findfirst(f, a)
+    return isnothing(i) ? nothing : a[i]
+end
+
+# Let F be the matrix with 1's on the antidiagonal. Then AF (or FA) is the
+# matrix A with columns (or rows) reversed. If (Q, R) = AF is the QR
+# decomposition of AF, then (QF, FRF) is the QL decomposition of A.
+function ql_slow(A)
+    AF = reduce(hcat, reverse(eachcol(A)))
+    Q, R = qr(AF)
+    # TODO: Perform these reversals in-place
+    QF = reduce(hcat, reverse(eachcol(collect(Q))))
+    RF = reduce(hcat, reverse(eachcol(collect(R))))
+    FRF = reduce(vcat, reverse(transpose.(eachrow(collect(RF)))))
+    return QF, FRF
+end

src/Operators/Rotation.jl

@@ -94,9 +94,7 @@ end
 # Return the closest matrix to R that is exactly orthogonal.
 #
 # This might be useful in cases where the crystal lattice vectors have numerical
-# error, and spglib produces only approximately orthogonal matrices (up to
-# symprec).  TODO: In the future, use spglib's feature "spg_standardize_cell()",
-# or write our own.
+# error, leading to inexact symops.  Alternatively one can call `standardize`.
 function to_orthogonal(R)
     # https://math.stackexchange.com/q/2215359
     U, _, V = svd(R)

src/Sunny.jl

@@ -40,7 +40,7 @@ include("Symmetry/AllowedCouplings.jl")
 include("Symmetry/AllowedAnisotropy.jl")
 include("Symmetry/Parsing.jl")
 include("Symmetry/Printing.jl")
-export Crystal, subcrystal, lattice_vectors, lattice_params, primitive_cell_shape, Bond,
+export Crystal, subcrystal, standardize, lattice_vectors, lattice_params, primitive_cell_shape, Bond,
     reference_bonds, print_site, print_bond, print_symmetry_table, print_suggested_frame
 
 include("Units.jl")
@@ -75,6 +75,9 @@ export minimize_energy!
 include("FormFactor.jl")
 export FormFactor
 
+include("MCIF.jl")
+export set_dipoles_from_mcif!
+
 include("SpinWaveTheory/SpinWaveTheory.jl")
 include("SpinWaveTheory/HamiltonianDipole.jl")
 include("SpinWaveTheory/HamiltonianSUN.jl")

src/Symmetry/Crystal.jl

@@ -10,11 +10,14 @@ An object describing a crystallographic unit cell and its space group symmetry.
 Constructors are as follows:
 
 
-    Crystal(filename; symprec=1e-5)
+    Crystal(filename; override_symmetry=false, symprec=nothing)
 
-Reads the crystal from a `.cif` file located at the path `filename`.  The
-optional parameter `symprec` controls the precision tolerance for spacegroup
-symmetries.
+Reads the crystal from a `.cif` file located at the path `filename`. If
+`override_symmetry=true`, the spacegroup will be inferred based on atom
+positions and the returned unit cell may be reduced in size. For an mCIF file,
+the return value is the magnetic supercell, unless `override_symmetry=true`. If
+a precision for spacegroup symmetries cannot be inferred from the cif file, it
+must be specified with `symprec`.
 
     Crystal(latvecs, positions; types=nothing, symprec=1e-5)
 
@@ -175,6 +178,78 @@ function sort_sites!(cryst::Crystal)
 end
 
 
+# Attempt to find a permutation matrix P (with sign flips) such that latvecs*P
+# has a more standard form.
+function permute_to_standardize_lattice_vectors(latvecs)
+    P = Mat3(I) # Iteratively build permutation matrix
+
+    for ij in eachindex(latvecs)
+        if norm(latvecs[ij]) < 1e-12
+            latvecs[ij] = 0
+        end
+    end
+
+    # Conventionally, a1 should be aligned with x direction
+    a1, a2, a3 = eachcol(latvecs)
+    if iszero(a2[2]) && iszero(a2[3]) && norm(a2) ≈ norm(a1)
+        P = P * [0 1 0; 1 0 0; 0 0 1] # Permute (a1, a2)
+    elseif iszero(a3[2]) && iszero(a3[3]) && norm(a3) ≈ norm(a1)
+        P = P * [0 0 1; 0 1 0; 1 0 0] # Permute (a1, a3)
+    end
+
+    # Conventionally, a2 should be in xy plane
+    _, a2, a3 = eachcol(latvecs * P)
+    if iszero(a3[3]) && norm(a3) ≈ norm(a2)
+        P = P * [1 0 0; 0 0 1; 0 1 0] # Permute (a2, a3)
+    end
+
+    # Flip columns so that the diagonal elements are positive
+    signs = [a < 0 ? -1 : 1 for a in diag(latvecs*P)]
+    P = P * Diagonal(signs)
+
+    # To preserve volume-orientation, we may need to flip one lattice
+    # vector. Pick a2 arbitrarily.
+    @assert det(P) ≈ 1 || det(P) ≈ -1
+    if det(P) ≈ -1
+        P = P * Diagonal([1,-1,1])
+    end
+
+    @assert det(P) ≈ 1
+    @assert P' ≈ inv(P)
+    return P
+end
+
+
+"""
+    standardize(cryst::Crystal; idealize=true)
+
+Return the symmetry-inferred standardized crystal unit cell. If `idealize=true`,
+then the lattice vectors and site positions will be adapted. See "definitions
+and conventions" of the [spglib
+documentation](https://spglib.readthedocs.io/en/stable/) for more information.
+"""
+function standardize(cryst::Crystal; idealize=true)
+    !isnothing(cryst.root) && error("Call this function on root crystal instead")
+
+    (; symprec) = cryst
+    cell = Spglib.Cell(cryst.latvecs, cryst.positions, cryst.types)
+    (; lattice, positions, atoms) = Spglib.standardize_cell(cell, symprec; no_idealize=!idealize)
+
+    if !idealize
+        # Here, spglib can return strange permutations of the lattice vectors.
+        # Attempt to permute lattice vectors back to a standard order (with
+        # sign-flips as needed).
+        P = permute_to_standardize_lattice_vectors(lattice)
+        # These transformations preserve global positions, `lattice * r`
+        lattice = lattice * P
+        positions = [P' * r for r in positions]
+    end
+
+    ret = crystal_from_inferred_symmetry(Mat3(lattice), Vec3.(positions), atoms; symprec)
+    sort_sites!(ret)
+    return ret
+end
+
 function crystal_from_inferred_symmetry(latvecs::Mat3, positions::Vector{Vec3}, types::Vector{String}; symprec=1e-5, check_cell=true)
     # Print a warning if non-conventional lattice vectors are detected.
     try cell_type(latvecs) catch e @warn e.msg end
@@ -199,21 +274,9 @@ function crystal_from_inferred_symmetry(latvecs::Mat3, positions::Vector{Vec3},
         ratio = length(positions) / d.n_std_atoms
         if ratio > 1
             ratio_str = number_to_simple_string(ratio; digits=4)
-            (a, b, c, α, β, γ) = number_to_math_string.(lattice_params(d.std_lattice); digits=8)
-            latstr = "lattice_vectors($a, $b, $c, $α, $β, $γ)"
-            position_str = "["*join(fractional_vec3_to_string.(d.std_positions), ", ")*"]"
-            crystal_build_str = if allequal(d.std_types)
-                "Crystal(latvecs, positions)"
-            else
-                types_str = "["*join(repr.(string.(d.std_types)), ", ")*"]"
-                "types = $types_str\n    Crystal(latvecs, positions; types)"
-            end
-            @warn """There exists another unit cell that is $ratio_str times smaller. Symmetry
-                     information may be incomplete. Consider using the smaller unit cell:
-                         latvecs = $latstr
-                         positions = $position_str
-                         $crystal_build_str
-                     Alternatively, site symmetry can be broken with a `types` argument."""
+            types_str = allunique(types) ? "" : "\nAlternatively, select `types` argument to break site symmetry."
+            @warn """The symmetry-inferred, conventional unit cell is $ratio_str times smaller. Obtain
+                     it with `standardize`.$types_str"""
         end
     end
 
@@ -349,6 +412,14 @@ function crystal_from_symbol(latvecs::Mat3, positions::Vector{Vec3}, types::Vect
     end
 end
 
+function symops_subset(symops1, symops2; symprec)
+    return all(symops1) do s
+        any(symops2) do s′
+            isapprox(s, s′; atol=symprec)
+        end
+    end
+end
+
 # Builds a crystal from an explicit set of symmetry operations and a minimal set of positions
 function crystal_from_symops(latvecs::Mat3, positions::Vector{Vec3}, types::Vector{String}, symops::Vector{SymOp}, spacegroup::String; symprec=1e-5)
     all_positions = Vec3[]
@@ -385,16 +456,8 @@ function crystal_from_symops(latvecs::Mat3, positions::Vector{Vec3}, types::Vect
     inferred = crystal_from_inferred_symmetry(latvecs, all_positions, all_types; symprec, check_cell=false)
 
     # Compare the inferred symops to the provided ones
-    is_subgroup = all(symops) do s
-        any(inferred.symops) do s′
-            isapprox(s, s′; atol=symprec)
-        end
-    end
-    is_supergroup = all(inferred.symops) do s
-        any(symops) do s′
-            isapprox(s, s′; atol=symprec)
-        end
-    end
+    is_subgroup = symops_subset(symops, inferred.symops; symprec)
+    is_supergroup = symops_subset(inferred.symops, symops; symprec)
 
     if !is_subgroup
         @warn """User provided symmetry operation could not be inferred by Spglib,
@@ -715,7 +778,7 @@ end
 
 function pyrochlore_crystal(; a=1.0)
     latvecs = lattice_vectors(a, a, a, 90, 90, 90)
-    return Crystal(latvecs, [[0, 0, 0]], 227, setting="2")
+    return Crystal(latvecs, [[5/8, 1/8, 1/8]], 227, setting="1")
 end
 
 function hyperkagome_crystal(; a=1.0)

src/Symmetry/MSymOp.jl

@@ -29,55 +29,3 @@ end
 function Base.:*(s1::MSymOp, s2::MSymOp)
     MSymOp(s1.R * s2.R, s1.T + s1.R * s2.T, s1.p*s2.p)
 end
-
-
-
-# This is a low-level function that will eventually be used to support the
-# reading of MCIF files. All keyword parameters follow the MCIF format: their
-# components should be interpreted in the coordinate system defined by the
-# lattice vectors of the supercell.
-function propagate_mcif_dipoles(sys; positions, dipoles, magn_operations, magn_centerings)
-    orig_crystal(sys) == sys.crystal || error("System employs a reshaped chemical cell")
-
-    # Defined such that r = D r′ where r is a "usual" position, while r′ is a
-    # position given in multiples of the system's (super) lattice vectors
-    D = Diagonal(Vec3(sys.latsize))
-
-    # Convert positions, dipoles, and symops to multiples of lattice vectors for
-    # the chemical cell
-    positions = [D * r for r in positions]
-    dipoles = [D * d for d in dipoles]
-    magn_operations = [MSymOp(D * s.R * inv(D), D * s.T, s.p) for s in magn_operations]
-    magn_centerings = [MSymOp(D * s.R * inv(D), D * s.T, s.p) for s in magn_centerings]
-
-    # Use the zero vector as a marker for unvisited sites
-    fill!(sys.dipoles, zero(Vec3))
-
-    for (r, d) in zip(positions, dipoles)
-        for s1 in magn_operations, s2 in magn_centerings
-            s = s1 * s2
-            # Transformed dipole. Because spin is a pseudo-vector, it is
-            # invariant to an inversion in space. For this reason, remove the
-            # determinant of R. If the parity p = ±1 is negative, this imposes
-            # time reversal, which effectively flips the spin.
-            d_new = (s.R * det(s.R) * s.p) * d
-            r_new = s.R * r + s.T
-            site = try 
-                position_to_site(sys, r_new)
-            catch _
-                rethrow(ErrorException("MCIF position $r_new is missing in the chemical cell"))
-            end
-
-            state = dipolar_state(sys, site, d_new)
-            if !iszero(sys.dipoles[site])
-                sys.dipoles[site] ≈ state.s || error("Conflicting dipoles at site $site")
-            end
-            setspin!(sys, state, site)
-        end
-    end
-
-    unvisited = [site for site in eachsite(sys) if iszero(sys.dipoles[site])]
-    if !isempty(unvisited)
-        error("Missing dipoles for sites $(collect(Tuple.(unvisited)))")
-    end
-end