view_crystal updated

ext/PlottingExt/ViewCrystal.jl

@@ -36,7 +36,7 @@ function reference_bonds_upto(cryst, nbonds, ndims)
     refbonds = filter(reference_bonds(cryst, max_dist)) do b
         return !(b.i == b.j && iszero(b.n))
     end
-    
+
     # Verify max_dist heuristic
     if length(refbonds) > 10nbonds
         @warn "Found $(length(refbonds)) bonds using max_dist of $max_dist"
@@ -100,7 +100,7 @@ function exchange_decomposition(J)
 
     # Now vecs is a pure rotation
     @assert vecs'*vecs ≈ I && det(vecs) ≈ 1
-    
+
     # Quaternion that rotates Cartesian coordinates into principle axes of J.
     axis, angle = Sunny.matrix_to_axis_angle(Mat3(vecs))
     q = iszero(axis) ? Makie.Quaternionf(0,0,0,1) : Makie.qrotation(axis, angle)
@@ -287,7 +287,12 @@ function label_atoms(cryst; ismagnetic)
             push!(ret, isempty(typ) ? "Position $rstr" : "'$typ' at $rstr")
         end
         if ismagnetic
-            # TODO: Show onsite couplings?
+            # See similar logic in print_site()
+            refatoms = [b.i for b in Sunny.reference_bonds(cryst, 0.0)]
+            i_ref = Sunny.findfirstval(i_ref -> Sunny.is_related_by_symmetry(cryst, i, i_ref), refatoms)
+            R_site = Sunny.rotation_between_sites(cryst, i, i_ref)
+            push!(ret, Sunny.allowed_g_tensor_string(cryst, i_ref; R_site, digits=8, atol=1e-12))
+            push!(ret, "See print_site($i; i_ref=$i_ref)")
         end
         join(ret, "\n")
     end
@@ -339,7 +344,7 @@ function draw_atoms_or_dipoles(; ax, full_crystal_toggle, dipole_menu, cryst, sy
                 labels = label_atoms(xtal; ismagnetic)[idxs]
                 inspector_label = (_plot, index, _position) -> labels[index]
             end
-            
+
             # Show dipoles. Mostly consistent with code in plot_spins.
             if !isnothing(sys) && ismagnetic
                 sites = Sunny.position_to_site.(Ref(sys), rs)

src/Symmetry/Crystal.jl

@@ -78,7 +78,7 @@ struct Crystal
     recipvecs :: Mat3                    # Reciprocal lattice vectors (conventional)
     positions :: Vector{Vec3}            # Positions in fractional coords
     types     :: Vector{String}          # Types
-    classes   :: Vector{Int}             # Class indices
+    classes   :: Vector{Int}             # Symmetry-equivalent class indices
     sg        :: Spacegroup              # Spacegroup symmetries and setting
     symprec   :: Float64                 # Tolerance to imperfections in symmetry
 end

src/Symmetry/Printing.jl

@@ -54,7 +54,6 @@ function fractional_mat3_to_string(m; digits=4, atol=1e-12)
     return "["*join(rowstrs, "; ")*"]"
 end
 
-
 # Like number_to_math_string(), but outputs a string that can be prefixed to a
 # variable name.
 function coefficient_to_math_string(x::T; digits=4, atol=1e-12) where T <: Real
@@ -110,6 +109,19 @@ function formatted_matrix(elemstrs::AbstractMatrix{String}; prefix)
     return "$prefix["*join(join.(eachrow(padded_elems), " "), spacing)*"]"
 end
 
+function int_to_underscore_string(x::Int)
+    subscripts = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉']
+    chars = collect(repr(x))
+    if chars[begin] == '-'
+        popfirst!(chars)
+        sign = "₋"
+    else
+        sign = ""
+    end
+    digits = map(c -> parse(Int, c), chars)
+    return sign * prod(subscripts[digits.+1])
+end
+
 
 """
     print_bond(cryst::Crystal, b::Bond; b_ref=b)
@@ -123,7 +135,7 @@ function print_bond(cryst::Crystal, b::Bond; b_ref=b, io=stdout)
     digits = 14
     # Tolerance below which coefficients are dropped
     atol = 1e-12
-    
+
     if b.i == b.j && iszero(b.n)
         print_site(cryst, b.i; i_ref=b.i, io)
     else
@@ -133,7 +145,7 @@ function print_bond(cryst::Crystal, b::Bond; b_ref=b, io=stdout)
         printstyled(io, "Bond($(b.i), $(b.j), $(b.n))"; bold=true, color=:underline)
         println(io)
         (m_i, m_j) = (coordination_number(cryst, b.i, b), coordination_number(cryst, b.j, b))
-        dist_str = number_to_simple_string(global_distance(cryst, b); digits, atol=1e-12)
+        dist_str = number_to_simple_string(global_distance(cryst, b); digits=10, atol=1e-12)
         if m_i == m_j
             println(io, "Distance $dist_str, coordination $m_i")
         else
@@ -192,9 +204,7 @@ anisotropies.
 """
 function print_suggested_frame(cryst::Crystal, i::Int)
     R = suggest_frame_for_atom(cryst, i)
-
     R_strs = [number_to_math_string(x; digits=14, atol=1e-12) for x in R]
-
     println(formatted_matrix(R_strs; prefix="R = "))
 end
 
@@ -221,47 +231,37 @@ function print_site(cryst, i; i_ref=i, R=Mat3(I), ks=[2,4,6], io=stdout)
         println(io, "Type '$(cryst.types[i])', position $(fractional_vec3_to_string(r)), multiplicity $m")
     end
 
-    # Tolerance below which coefficients are dropped
+    digits = 14
     atol = 1e-12
-    # How many digits to use in printing coefficients
-    digits = 10
 
+    R_site = rotation_between_sites(cryst, i, i_ref)
+    println(io, allowed_g_tensor_string(cryst, i_ref; R_global, R_site, digits, atol))
+    println(io, allowed_anisotropy_string(cryst, i_ref; R_global, R_site, digits, atol, ks))
+end
+
+
+function rotation_between_sites(cryst, i, i_ref)
     # Rotation that maps from i_ref to i
     if i == i_ref
-        R_site = Mat3(I)
+        return Mat3(I)
     else
         syms = symmetries_between_atoms(cryst, i, i_ref)
         isempty(syms) && error("Atoms $i and $i_ref are not symmetry equivalent.")
-        R_site = cryst.latvecs * first(syms).R * inv(cryst.latvecs)
-        R_site *= det(R_site) # Remove possible inversion (appropriate for spin pseudo-vector)
+        R = cryst.latvecs * first(syms).R * inv(cryst.latvecs)
+        return R * det(R) # Remove possible inversion (appropriate for spin pseudo-vector)
     end
+end
 
+function allowed_g_tensor_string(cryst, i_ref; R_global=Mat3(I), R_site, digits, atol)
     basis = basis_for_symmetry_allowed_couplings(cryst, Bond(i_ref, i_ref, [0, 0, 0]); R_global)
     basis = map(basis) do b
         R_site * b * R_site' # == transform_coupling_by_symmetry(b, R_site, true)
     end
-
     basis_strs = coupling_basis_strings(zip('A':'Z', basis); digits, atol)
-    println(io, formatted_matrix(basis_strs; prefix="Allowed g-tensor: "))
-
-    print_allowed_anisotropy(cryst, i_ref; R_global, R_site, atol, digits, ks, io)
+    return formatted_matrix(basis_strs; prefix="Allowed g-tensor: ")
 end
 
-function int_to_underscore_string(x::Int)
-    subscripts = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉']
-    chars = collect(repr(x))
-    if chars[begin] == '-'
-        popfirst!(chars)
-        sign = "-"
-    else
-        sign = ""
-    end
-    digits = map(c -> parse(Int, c), chars)
-    return sign * prod(subscripts[digits.+1])
-end
-
-
-function print_allowed_anisotropy(cryst::Crystal, i_ref::Int; R_global::Mat3, R_site::Mat3, atol, digits, ks, io=stdout)
+function allowed_anisotropy_string(cryst::Crystal, i_ref::Int; R_global::Mat3, R_site::Mat3, digits, atol, ks)
     prefix="    "
 
     lines = String[]
@@ -302,10 +302,10 @@ function print_allowed_anisotropy(cryst::Crystal, i_ref::Int; R_global::Mat3, R_
             push!(lines, prefix * join(terms, " + "))
         end
     end
-    println(io, "Allowed anisotropy in Stevens operators:")
-    println(io, join(lines, " +\n"))
 
+    ret = "Allowed anisotropy in Stevens operators:\n" * join(lines, " +\n")
     if R_global != I
-        println(io, "Modified reference frame! Use R*g*R' or rotate_operator(op, R).")
+        ret *= "\nModified reference frame! Use R*g*R' or rotate_operator(op, R)."
     end
+    return ret
 end

src/Symmetry/SymmetryAnalysis.jl

@@ -109,6 +109,11 @@ function symmetries_between_bonds(cryst::Crystal, b1::BondPos, b2::BondPos)
     return ret
 end
 
+# Are i1 and i2 symmetry-equivalent sites?
+function is_related_by_symmetry(cryst::Crystal, i1::Int, i2::Int)
+    return cryst.classes[i1] == cryst.classes[i2]
+end
+
 # Is there a symmetry operation that transforms `b1` into either `b2` or its
 # reverse?
 function is_related_by_symmetry(cryst::Crystal, b1::Bond, b2::Bond)

test/test_symmetry.jl

@@ -350,37 +350,37 @@ end
             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
+        Distance 0.3535533906, coordination 6
         Connects [0, 0, 0] to [1/4, 1/4, 0]
         Allowed exchange matrix: [A C -D
                                   C A -D
                                   D D  B]
         Allowed DM vector: [-D D 0]
 
         Bond(3, 5, [0, 0, 0])
-        Distance 0.61237243569579, coordination 12
+        Distance 0.6123724357, coordination 12
         Connects [1/2, 1/2, 0] to [1/4, 0, 1/4]
         Allowed exchange matrix: [  A  C-E  D-F
                                   C+E    B -C+E
                                   D+F -C-E    A]
         Allowed DM vector: [E F -E]
 
         Bond(1, 3, [-1, 0, 0])
-        Distance 0.70710678118655, coordination 6
+        Distance 0.7071067812, coordination 6
         Connects [0, 0, 0] to [-1/2, 1/2, 0]
         Allowed exchange matrix: [A D C
                                   D A C
                                   C C B]
 
         Bond(1, 3, [0, 0, 0])
-        Distance 0.70710678118655, coordination 6
+        Distance 0.7071067812, coordination 6
         Connects [0, 0, 0] to [1/2, 1/2, 0]
         Allowed exchange matrix: [A D C
                                   D A C
                                   C C B]
 
         Bond(1, 2, [-1, 0, 0])
-        Distance 0.79056941504209, coordination 12
+        Distance 0.790569415, coordination 12
         Connects [0, 0, 0] to [-3/4, 1/4, 0]
         Allowed exchange matrix: [A  D -F
                                   D  B  E
@@ -453,6 +453,15 @@ end
     positions = Sunny.fcc_crystal().positions
     cryst = Crystal(latvecs, positions; types = ["A", "B", "B", "B"])
 
+    capt = IOCapture.capture() do
+        print_suggested_frame(cryst, 1)
+    end
+    @test capt.output == """
+        R = [0.70803177573023 -0.70618057503467                 0
+             0.40878233631266  0.40985392929053 -0.81542428107328
+             0.57583678770556  0.57734630170186  0.57886375066688]
+        """
+
     R = [0.70803177573023 -0.70618057503467                 0
          0.40878233631266  0.40985392929053 -0.81542428107328
          0.57583678770556  0.57734630170186  0.57886375066688]