Suggest standardize when a nonconventional cell is inferred (#461)

Also accept nonconventional hexagonal and tetragonal cells for lower-symmetry spacegroups (triclinic, monoclinic, and orthorhombic).

Author: kbarros

docs/src/versions.md

@@ -31,6 +31,8 @@
 * Introduce [`set_spin_s_at!`](@ref) to set the local quantum spin-``s``
   ([#454](@ref)).
 * Add missing 3D support for [`q_space_grid`](@ref) ([#457](@ref)).
+* If user-provided lattice vectors do not match crystallographic conventions,
+  suggest the use of [`standardize`](@ref) ([#461](@ref)).
 * Fixes to [`SpinWaveTheoryKPM`](@ref) ([#462](@ref)).
 
 ## v0.7.8

src/Symmetry/AllowedCouplings.jl

@@ -136,11 +136,11 @@ end
 # Check whether a coupling matrix J is consistent with symmetries of a bond
 function is_coupling_valid(cryst::Crystal, b::BondPos, J)
     J isa Number && return true
-    
+
     for (symop, parity) in symmetries_between_bonds(cryst, b, b)
         R = cryst.latvecs * symop.R * inv(cryst.latvecs)
         J′ = transform_coupling_by_symmetry(J, R*det(R), parity)
-        # For non-conventional unit cells, the rotation matrices R produced by
+        # For nonconventional unit cells, the rotation matrices R produced by
         # spglib might only be accurate to about 11 digits. If an ambiguous
         # situation is detected, throw an informative error.
         if !isapprox(J, J′; rtol=1e-10)

src/Symmetry/Crystal.jl

@@ -90,7 +90,12 @@ end
 # representing fractions (between 0 and 1) of the lattice vectors `latvecs`.
 # All symmetry information is automatically inferred.
 function Crystal(latvecs, positions; types::Union{Nothing, Vector{String}}=nothing, symprec=1e-5)
-    print_crystal_warnings(latvecs, positions)
+    if length(positions) >= 100
+        @warn """Very large crystal cells are not recommended. To model chemical inhomogeneity:
+                 1. Create the `Crystal` with an idealized chemical cell.
+                 2. Create the `System` with large volume dimensions, `dims`.
+                 3. Use `to_inhomogeneous` and related functions to set inhomogeneities."""
+    end
     latvecs = convert(Mat3, latvecs)
     positions = [convert(Vec3, p) for p in positions]
     if isnothing(types)
@@ -106,21 +111,23 @@ function Crystal(latvecs, positions, symbol::Union{Int, String}; types::Union{No
                  setting=nothing, choice::Union{Nothing, String}=nothing, symprec=1e-5)
     if !isnothing(setting)
         if symbol isa Integer
-            @warn """`setting` argument is deprecated! Omit to get the ITA standard setting. 
-                     Alternatively, use `choice` instead."""
+            @warn "`setting` argument is deprecated! Omit to get the ITA standard setting or use `choice`."
         else
             @warn "`setting` argument is deprecated! Use `choice` instead."
         end
         choice = setting
     end
 
-    print_crystal_warnings(latvecs, positions)
     latvecs = Mat3(latvecs)
     positions = Vec3.(positions)
     if isnothing(types)
         types = fill("", length(positions))
     end
 
+    # Enforce right-handedness convention for aesthetics
+    det(latvecs) < 0 && error("Lattice vectors are not right-handed.")
+
+    # Will report an informative error if latvecs not consistent with symbol
     sgt = unique_spacegroup_type(symbol, latvecs; choice)
     sg = Spacegroup(Int(sgt.hall_number))
     return crystal_from_spacegroup(latvecs, positions, types, sg; symprec)
@@ -148,17 +155,6 @@ function spg_get_dataset_scaled(cell, symprec)
 end
 
 
-function print_crystal_warnings(latvecs, positions)
-    det(latvecs) < 0 && @warn "Lattice vectors are not right-handed."
-    if length(positions) >= 100
-        @info """This a very large crystallographic cell, which Sunny does not handle well. If
-                 the intention is to model chemical inhomogeneity, the recommended steps are:
-                 (1) Create a `Crystal` with idealized chemical cell. (2) Create a `System`
-                 with many cells, as set by the `dims` parameter. (3) Use `to_inhomogeneous`
-                 and related functions to introduce model inhomogeneities."""
-    end
-end
-
 """
     natoms(cryst::Crystal)
 
@@ -302,29 +298,47 @@ function conventionalize_setting(latvecs::Mat3, setting::SymOp, sgnum::Int)
 end
 
 function crystal_from_inferred_symmetry(latvecs::Mat3, positions::Vector{Vec3}, types::Vector{String}; symprec, suppress_warnings=false)
-    # Print a warning if non-conventional lattice vectors are detected.
-    try cell_type(latvecs) catch e @warn e.msg end
-
     validate_positions(positions; symprec)
 
-    cell = Spglib.Cell(latvecs, positions, types)
-    d = spg_get_dataset_scaled(cell, symprec)
+    d = spg_get_dataset_scaled(Spglib.Cell(latvecs, positions, types), symprec)
     classes = d.crystallographic_orbits
     # classes = d.equivalent_atoms
     symops = SymOp.(d.rotations, d.translations)
     label = spacegroup_label(Int(d.hall_number))
     sgnum = Int(d.spacegroup_number)
+    cell_std = cell_type(standard_setting[sgnum])
     setting = mapping_to_standard_setting_from_spglib_dataset(d)
     setting = conventionalize_setting(latvecs, setting, sgnum)
 
-    # Spglib-inferred symops are unreliable for large cells. Empty the list to
-    # prevent faulty symmetry analysis.
-    if length(positions) > d.n_std_atoms
-        if !suppress_warnings
-            types_str = allunique(types) ? "" : " Alternatively, break site symmetry with `types=[...]`."
-            @warn "This cell is non-standard and missing symmetry information. Consider `standardize`.$types_str"
+    # Spglib-inferred symops are unreliable for large cells. Empty the list
+    # to prevent faulty symmetry analysis.
+    volume_ratio = length(positions) // d.n_std_atoms
+    volume_ratio > 1 && empty!(symops)
+    volume_ratio_str = number_to_math_string(volume_ratio)
+
+    if !suppress_warnings
+        if volume_ratio > 1
+            volume_ratio_str = number_to_math_string(volume_ratio)
+            types_str = allunique(types) ? "" : " or distinct `types=[...]`"
+            @error "Symmetry analysis disabled! Cell is $volume_ratio_str times too large. Fix with `standardize`$types_str."
+        elseif volume_ratio < 1
+            @info "Cell is $volume_ratio_str the standard size for spacegroup $sgnum. Consider `standardize`."
+        elseif cell_std in (tetragonal, hexagonal, cubic) && setting.R ≉ I
+            # For a tetragonal, hexagonal, trigonal or cubic space group, the
+            # provided cell should already be Spglib-idealized.
+            @info "Nonstandard $cell_std cell for spacegroup $sgnum. Consider `standardize`."
+        elseif !(cell_type(latvecs) in all_compatible_cells(cell_std))
+            # For triclinic, monoclinic, and orthorhombic spacegroups, check
+            # that the provided cell is symmetry-consistent (not necessarily
+            # Spglib-idealized). For orthorhombic, arbitrary axis ordering is
+            # acceptable. For monoclinic, arbitrary β magnitude is acceptable,
+            # whereas Spglib finds a cell with 90 < β < 120. For triclinic,
+            # Spglib's idealization algorithm is even more complex, involving
+            # Niggli reduction.
+            @info "Nonstandard $cell_std cell for spacegroup $sgnum. Consider `standardize`."
+        elseif det(latvecs) < 0
+            @info "Lattice vectors are not right-handed. Consider `standardize`."
         end
-        symops = SymOp[]
     end
 
     sg = Spacegroup(symops, label, sgnum, setting)
@@ -511,7 +525,7 @@ end
 
 # Indices of atoms that wrap into unique positions within the primitive cell.
 # The choice is ambiguous because the un-wrapped positions are not restricted to
-# `latvecs * primitive_cell()`.
+# the parallelpiped spanned by `latvecs * primitive_cell()`.
 function primitive_atoms_arbitrary(cryst::Crystal)
     P = primitive_cell(cryst)
     P ≈ I && return collect(1:natoms(cryst))

src/Symmetry/LatticeUtils.jl

@@ -67,19 +67,21 @@ end
 # spacegroups are characterized by a 3-fold rotational symmetry. All trigonal
 # spacegroups (143-167) admit a hexagonal setting. Some of these (146, 148, 155,
 # 160, 161, 166, 167) additionally admit a rhombohedral setting.
+#
+# The "alt" suffix indicates deviation from the ITA lattice vector conventions.
 @enum CellType begin
     triclinic
     monoclinic
     orthorhombic
     tetragonal
+    tetragonal_alt
     rhombohedral
     hexagonal
+    hexagonal_alt
     cubic
 end
 
-# Infer the CellType (lattice system) from lattice vectors. Report an error if
-# the unit cell is not in conventional form, which would invalidate the table of
-# symops for a given Hall number.
+# Infer the CellType (lattice system) from lattice vectors.
 function cell_type(latvecs)
     a, b, c, α, β, γ = lattice_params(latvecs)
 
@@ -95,7 +97,7 @@ function cell_type(latvecs)
         if a ≈ b
             return tetragonal
         elseif b ≈ c || c ≈ a
-            error("Found a nonconventional tetragonal unit cell. Consider using `lattice_vectors(a, a, c, 90, 90, 90)`.")
+            return tetragonal_alt # nonconventional
         else
             return orthorhombic
         end
@@ -107,7 +109,7 @@ function cell_type(latvecs)
         if γ ≈ 120
             return hexagonal
         else
-            error("Found a nonconventional hexagonal unit cell. Consider using `lattice_vectors(a, a, c, 90, 90, 120)`.")
+            return hexagonal_alt # nonconventional
         end
     end
 
@@ -166,11 +168,11 @@ end
 
 function all_compatible_cells(cell::CellType)
     if cell == triclinic
-        [triclinic, monoclinic, orthorhombic, tetragonal, rhombohedral, hexagonal, cubic]
+        [triclinic, monoclinic, orthorhombic, tetragonal, tetragonal_alt, rhombohedral, hexagonal, hexagonal_alt, cubic]
     elseif cell == monoclinic
-        [monoclinic, orthorhombic, tetragonal, hexagonal, cubic]
+        [monoclinic, orthorhombic, tetragonal, tetragonal_alt, hexagonal, hexagonal_alt, cubic]
     elseif cell == orthorhombic
-        [orthorhombic, tetragonal, cubic]
+        [orthorhombic, tetragonal, tetragonal_alt, cubic]
     elseif cell == tetragonal
         [tetragonal, cubic]
     elseif cell == rhombohedral

src/Symmetry/SpacegroupData.jl

@@ -40,7 +40,7 @@ function suggest_alternative_symbols(sgts)
 end
 
 # Get the single SpacegroupType associated with symbol that is valid for
-# latvecs. Throw an informative error if the setting is ambiguous.
+# latvecs. Throw an informative error if the setting is invalid or ambiguous.
 function unique_spacegroup_type(symbol, latvecs; choice=nothing)
     # All settings for symbol
     sgts = all_spacegroup_types_for_symbol(symbol)
@@ -62,9 +62,15 @@ function unique_spacegroup_type(symbol, latvecs; choice=nothing)
     sgts = filter(sgts) do sgt
         return cell in all_compatible_cells(cell_type(sgt))
     end
-    if isempty(sgts)
-        expected = join(string.(allowed_cells), " or ")
-        error("Expected $expected lattice system but got $cell")
+
+    if cell == tetragonal_alt && cell_type(hall_std) == tetragonal
+        error("Use a conventional tetragonal cell: `lattice_vectors(a, a, c, 90, 90, 90)`")
+    elseif cell == hexagonal_alt && cell_type(hall_std) == hexagonal
+        error("Use a conventional hexagonal cell: `lattice_vectors(a, a, c, 90, 90, 120)`")
+    elseif isempty(sgts)
+        allowed_str = join(string.(allowed_cells), " or ")
+        received = only(replace([cell], hexagonal_alt => hexagonal, tetragonal_alt => tetragonal))
+        error("Expected $allowed_str cell but got $received")
     end
 
     # For monoclinic lattice systems, filter by axis setting

test/test_correlation_sampling.jl

@@ -1,18 +1,17 @@
 @testitem "Correlation sampling" begin
     using LinearAlgebra
 
+    # FCC with nonstandard, primitive lattice vectors
+    latvecs = [[1, 1, 0] [0, 1, 1] [1, 0, 1]] / 2
+    positions = [[0, 0, 0]]
+    msg = "Cell is 1/4 the standard size for spacegroup 225. Consider `standardize`."
+    cryst = @test_logs (:info, msg) Crystal(latvecs, positions)
+
     function simple_model_fcc(; mode, seed=111)
-        dims = (4, 4, 4)
         J = 1.0
-
-        # FCC with nonstandard, primitive lattice vectors
-        latvecs = [[1, 1, 0] [0, 1, 1] [1, 0, 1]] / 2
-        positions = [[0, 0, 0]]
-        cryst = Crystal(latvecs, positions)
-
         s = mode==:SUN ? 1/2 : 1
         κ = mode==:SUN ? 2 : 1
-        sys = System(cryst, [1 => Moment(; s, g=2)], mode; dims, seed)
+        sys = System(cryst, [1 => Moment(; s, g=2)], mode; dims=(4, 4, 4), seed)
         sys.κs .= κ
         set_exchange!(sys, J, Bond(1, 1, [1, 0, 0]))
         return sys

test/test_lswt.jl

@@ -1,11 +1,12 @@
-@testitem "Kitchen Sink" begin
+@testitem "Kitchen sink" begin
     using LinearAlgebra
 
     # Pyrochlore with nonstandard, primitive lattice vectors
     latvecs = [[1, 1, 0] [1, 0, 1] [0, 1, 1]] / 2
     positions = [[5, 5, 1], [5, 1, 5], [1, 5, 5], [5, 5, 5]] / 8
 
-    cryst = @test_logs (:warn, "Lattice vectors are not right-handed.") Crystal(latvecs, positions)
+    msg = "Cell is 1/4 the standard size for spacegroup 227. Consider `standardize`."
+    cryst = @test_logs (:info, msg) Crystal(latvecs, positions)
     natoms = Sunny.natoms(cryst)
 
     moments = [1 => Moment(s=5/2, g=7.2)]
@@ -192,7 +193,8 @@ end
     using LinearAlgebra
 
     latvecs = lattice_vectors(1, 1, 1, 90, 90, 90)
-    cryst = Crystal(latvecs, [[0,0,0], [0.4,0,0]]; types=["A", "B"])
+    msg = "Nonstandard tetragonal cell for spacegroup 99. Consider `standardize`."
+    cryst = @test_logs (:info, msg) Crystal(latvecs, [[0,0,0], [0.4,0,0]]; types=["A", "B"])
 
     sys = System(cryst, [1 => Moment(s=1, g=2), 2 => Moment(s=2, g=2)], :dipole)
     set_pair_coupling!(sys, (S1, S2) -> +(S1'*diagm([2,-1,-1])*S1)*(S2'*diagm([2,-1,-1])*S2), Bond(1, 2, [0,0,0]))

test/test_rescaling.jl

@@ -110,8 +110,8 @@ end
 
 @testitem "Anisotropy SU(N) equivalence" begin
     latvecs = lattice_vectors(1.0, 1.1, 1.0, 90, 90, 90)
-    msg = "Found a nonconventional tetragonal unit cell. Consider using `lattice_vectors(a, a, c, 90, 90, 90)`."
-    cryst = @test_logs (:warn, msg) Crystal(latvecs, [[0, 0, 0]])
+    msg = "Nonstandard tetragonal cell for spacegroup 123. Consider `standardize`."
+    cryst = @test_logs (:info, msg) Crystal(latvecs, [[0, 0, 0]])
 
     # Dipole system with renormalized anisotropy
     sys0 = System(cryst, [1 => Moment(s=3, g=2)], :dipole)
@@ -169,7 +169,8 @@ end
 @testitem "Biquadratic renormalization 2" begin
     # Simple dimer model
     latvecs = lattice_vectors(1, 1, 1, 90, 90, 90)
-    cryst = Crystal(latvecs, [[0, 0, 0], [0.3, 0, 0]]; types=["A", "B"])
+    msg = "Nonstandard tetragonal cell for spacegroup 99. Consider `standardize`."
+    cryst = @test_logs (:info, msg) Crystal(latvecs, [[0, 0, 0], [0.3, 0, 0]]; types=["A", "B"])
     s1 = 3/2
     s2 = 2
     v1 = randn(3)