Avoid crashes when reporting symmetry-allowed interactions (SunnySuit…

…e#317)

Fixes SunnySuite#260.

Project.toml

@@ -1,7 +1,7 @@
 name = "Sunny"
 uuid = "2b4a2ac8-8f8b-43e8-abf4-3cb0c45e8736"
 authors = ["The Sunny team"]
-version = "0.7.2"
+version = "0.7.3"
 
 [deps]
 CrystalInfoFramework = "6007d9b0-c6b2-11e8-0510-1d10e825f3f1"

docs/src/versions.md

@@ -1,5 +1,11 @@
 # Version History
 
+## v0.7.3
+(In development)
+
+* Fix error in `print_symmetry_table` for slightly-distorted crystal cells ([PR
+  #317](https://github.com/SunnySuite/Sunny.jl/pull/317)).
+
 ## v0.7.2
 (Sep 11, 2024)
 

src/Operators/Rotation.jl

@@ -137,7 +137,6 @@ function unitary_tensor_for_rotation(R::Mat3; Ns)
     end
 end
 
-# TODO: Replace this with a function that takes generators.
 """
     rotate_operator(A, R)
 

src/Symmetry/AllowedAnisotropy.jl

@@ -61,38 +61,23 @@ function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int,
     D = unitary_irrep_for_rotation(R; N=2k+1)
     C = D * C
 
-    # Find an orthonormal basis for the columns of A, discarding linearly
-    # dependent columns.
-    C = colspace(C; atol=1e-12)
-
-    # It is tempting to sparsify here to make the ouput look nicer. Don't do
-    # this because (empirically) it is observed to significantly degrade
-    # accuracy in stevens_basis_for_symmetry_allowed_anisotropies().
-
-    # C = sparsify_columns(C; atol=1e-12)
-
-    return C
-end
-
-function stevens_basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int, R::Mat3)
-    # Each column of C represents a coefficient vector c that can be contracted
-    # with spherical tensors T to realize an allowed anisotropy, Λ = cᵀ T.
-    C = basis_for_symmetry_allowed_anisotropies(cryst, i; k, R)
+    # Transform columns c of C to columns b of B such that 𝒜 = cᵀ T = bᵀ 𝒪.
+    # Effectively, this reexpresses the symmetry-allowed operators 𝒜 in the
+    # basis of Stevens operators 𝒪.
+    B = mapslices(C; dims=1) do c
+        transform_spherical_to_stevens_coefficients(k, c)
+    end
 
-    # Transform each column c to coefficients b that satisfy bᵀ 𝒪 = cᵀ T
-    B = [transform_spherical_to_stevens_coefficients(k, c) for c in eachcol(C)]
+    # If 𝒜 is symmetry-allowed, then its Hermitian and anti-Hermitian parts are
+    # independently symmetry-allowed. These are associated with the real and
+    # imaginary parts of B. Use the real and imaginary parts of B to construct
+    # an over-complete set of symmetry-allowed operators that are guaranteed to
+    # be Hermitian. Create a widened real matrix from these two parts, and
+    # eliminate linearly-dependent vectors from the column space.
+    B = colspace(hcat(real(B), imag(B)); atol=1e-12)
 
-    # Concatenate columns into single matrix
-    B = reduce(hcat, B; init=zeros(ComplexF64, 2k+1,0))
-
     # Find linear combination of columns that sparsifies B
-    B = sparsify_columns(B; atol=1e-12)
-
-    # All coefficients must now be real
-    @assert norm(imag(B)) < 1e-12
-    B = real(B)
-
-    return B
+    return sparsify_columns(B; atol=1e-12)
 end
 
 

src/Symmetry/Printing.jl

@@ -190,7 +190,7 @@ end
 
 
 """
-    print_suggested_frame(cryst, i; digits=4)
+    print_suggested_frame(cryst, i)
 
 Print a suggested reference frame, as a rotation matrix `R`, that can be used as
 input to `print_site()`. The purpose is to simplify the description of allowed
@@ -227,12 +227,11 @@ function print_site(cryst, i; R=Mat3(I), ks=[2,4,6], io=stdout)
     # Tolerance below which coefficients are dropped
     atol = 1e-12
     # How many digits to use in printing coefficients
-    digits = 14
+    digits = 10
 
     R = convert(Mat3, R) # Rotate to frame of R
     basis = basis_for_symmetry_allowed_couplings(cryst, Bond(i, i, [0,0,0]))
-    # TODO: `R` should be passed to `basis_for_symmetry_allowed_couplings` to
-    # get a nicer basis.
+    # TODO: `basis_for_symmetry_allowed_couplings` should accept R instead
     basis = [R * b * R' for b in basis]
     basis_strs = coupling_basis_strings(zip('A':'Z', basis); digits, atol)
     println(io, formatted_matrix(basis_strs; prefix="Allowed g-tensor: "))
@@ -260,23 +259,28 @@ function print_allowed_anisotropy(cryst::Crystal, i::Int; R::Mat3, atol, digits,
     lines = String[]
     cnt = 1
     for k in ks
-        B = stevens_basis_for_symmetry_allowed_anisotropies(cryst, i; k, R)
+        B = basis_for_symmetry_allowed_anisotropies(cryst, i; k, R)
 
         if size(B, 2) > 0
             terms = String[]
             for b in reverse(collect(eachcol(B)))
+                # rescale column so that the largest component is 1
+                b /= argmax(abs, b)
 
-                if any(x -> 1e-12 < abs(x) < 1e-6, b)
+                if any(x -> atol < abs(x) < sqrt(atol), b)
                     @info """Found a very small but nonzero expansion coefficient.
                              This may indicate a slightly misaligned reference frame."""
                 end
 
-                # rescale column by its minimum nonzero value
-                _, min = findmin(b) do x
-                    abs(x) < 1e-12 ? Inf : abs(x)
+                # rescale by up to 60× if it makes all coefficients integer
+                denoms = denominator.(rationalize.(b; tol=atol))
+                if all(<=(60), denoms)
+                    factor = lcm(denominator.(rationalize.(b; tol=atol)))
+                    if factor <= 60
+                        b .*= factor
+                    end
                 end
-                b /= b[min]
-
+
                 # reverse b elements to print q-components in ascending order, q=-k...k
                 ops = String[]
                 for (b_q, q) in zip(reverse(b), -k:k)
@@ -299,8 +303,6 @@ function print_allowed_anisotropy(cryst::Crystal, i::Int; R::Mat3, atol, digits,
     println(io, join(lines, " +\n"))
 
     if R != I
-        println(io)
-        println(io, "Modified reference frame! Transform using `rotate_operator(op, R)` where")
-        println(io, formatted_matrix(number_to_math_string.(R); prefix="R = "))
+        println(io, "Modified reference frame! Transform using `rotate_operator(op, R)`.")
     end
 end

test/test_symmetry.jl

@@ -254,8 +254,8 @@ end
                            B B A]
         Allowed anisotropy in Stevens operators:
             c₁*(𝒪[2,-2]+2𝒪[2,-1]+2𝒪[2,1]) +
-            c₂*(-7𝒪[4,-3]-2𝒪[4,-2]+𝒪[4,-1]+𝒪[4,1]+7𝒪[4,3]) + c₃*(𝒪[4,0]+5𝒪[4,4]) +
-            c₄*(-11𝒪[6,-6]-8𝒪[6,-3]+𝒪[6,-2]-8𝒪[6,-1]-8𝒪[6,1]+8𝒪[6,3]) + c₅*(𝒪[6,0]-21𝒪[6,4]) + c₆*((9/5)𝒪[6,-6]+(24/5)𝒪[6,-5]+𝒪[6,-2]+(8/5)𝒪[6,-1]+(8/5)𝒪[6,1]+(24/5)𝒪[6,5])
+            c₂*(7𝒪[4,-3]+2𝒪[4,-2]-𝒪[4,-1]-𝒪[4,1]-7𝒪[4,3]) + c₃*(𝒪[4,0]+5𝒪[4,4]) +
+            c₄*(11𝒪[6,-6]+8𝒪[6,-3]-𝒪[6,-2]+8𝒪[6,-1]+8𝒪[6,1]-8𝒪[6,3]) + c₅*(-𝒪[6,0]+21𝒪[6,4]) + c₆*(9𝒪[6,-6]+24𝒪[6,-5]+5𝒪[6,-2]+8𝒪[6,-1]+8𝒪[6,1]+24𝒪[6,5])
         
         Bond(1, 2, [0, 0, 0])
         Distance 0.35355339059327, coordination 6
@@ -320,11 +320,7 @@ end
             c₁*𝒪[2,0] +
             c₂*𝒪[4,-3] + c₃*𝒪[4,0] +
             c₄*𝒪[6,-3] + c₅*𝒪[6,0] + c₆*𝒪[6,6]
-        
-        Modified reference frame! Transform using `rotate_operator(op, R)` where
-        R = [1/√2      0  1/√2
-             1/√6 -√2/√3 -1/√6
-             1/√3   1/√3 -1/√3]
+        Modified reference frame! Transform using `rotate_operator(op, R)`.
         """
 
     capt = IOCapture.capture() do
@@ -337,6 +333,51 @@ end
                 0     0 -1
              1/√2  1/√2  0]
         """
+
+    # Test for https://github.com/SunnySuite/Sunny.jl/issues/260
+
+    a = 6.22
+    distortion = 0.15
+    latvecs = lattice_vectors(a, a, a, 90+distortion, 90+distortion, 90+distortion)
+    positions = Sunny.fcc_crystal().positions
+    crystal = Crystal(latvecs, positions; types = ["A", "B", "B", "B"])
+
+    R = [0.70803177573023 -0.70618057503467                 0
+         0.40878233631266  0.40985392929053 -0.81542428107328
+         0.57583678770556  0.57734630170186  0.57886375066688]
+    capt = IOCapture.capture() do
+        print_site(crystal, 1; R)
+        print_site(crystal, 2; R)
+    end
+    @test capt.output == """
+        Atom 1
+        Type 'A', position [0, 0, 0], multiplicity 1
+        Allowed g-tensor: [A-0.9921699533B               0             0
+                                         0 A-0.9921699533B             0
+                                         0               0 A+2.00000689B]
+        Allowed anisotropy in Stevens operators:
+            c₁*𝒪[2,0] +
+            c₂*𝒪[4,-3] + c₃*𝒪[4,0] +
+            c₄*𝒪[6,-3] + c₅*𝒪[6,0] + c₆*𝒪[6,6]
+        Modified reference frame! Transform using `rotate_operator(op, R)`.
+        Atom 2
+        Type 'B', position [1/2, 1/2, 0], multiplicity 3
+        Allowed g-tensor: [A-0.9973854271B                                                                   0                                                                   0
+                                         0                0.3350832418A+0.335961638B+0.6649167582C-1.33332874D 0.4720195577A+0.4732569225B-0.4720195577C-0.4658456409D-1.41236599E
+                                         0 0.4720195577A+0.4732569225B-0.4720195577C-0.4658456409D+1.41236599E               0.6649167582A+0.6666597888B+0.3350832418C+1.33332874D]
+        Allowed anisotropy in Stevens operators:
+            c₁*𝒪[2,-1] + c₂*𝒪[2,0] + c₃*𝒪[2,2] +
+            c₄*𝒪[4,-3] + c₅*𝒪[4,-1] + c₆*𝒪[4,0] + c₇*𝒪[4,2] + c₈*𝒪[4,4] +
+            c₉*𝒪[6,-5] + c₁₀*𝒪[6,-3] + c₁₁*𝒪[6,-1] + c₁₂*𝒪[6,0] + c₁₃*𝒪[6,2] + c₁₄*𝒪[6,4] + c₁₅*𝒪[6,6]
+        Modified reference frame! Transform using `rotate_operator(op, R)`.
+        """
+
+    # These operators should be symmetry allowed
+    s = 4
+    sys = System(crystal, [1 => Moment(; s, g=2), 2 => Moment(; s, g=2)], :dipole)
+    O = stevens_matrices(s)
+    set_onsite_coupling!(sys, O[6,-1]+0.997385420O[6,1], 2)
+    set_onsite_coupling!(sys, rotate_operator(O[6,2], R), 2)
 end