diff --git a/ext/PlottingExt/ViewCrystal.jl b/ext/PlottingExt/ViewCrystal.jl index 983439127..29a374b7a 100644 --- a/ext/PlottingExt/ViewCrystal.jl +++ b/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) diff --git a/src/Symmetry/Crystal.jl b/src/Symmetry/Crystal.jl index 8a18a2eed..7a13e5a05 100644 --- a/src/Symmetry/Crystal.jl +++ b/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 diff --git a/src/Symmetry/Printing.jl b/src/Symmetry/Printing.jl index e4ecc083d..89205325b 100644 --- a/src/Symmetry/Printing.jl +++ b/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 diff --git a/src/Symmetry/SymmetryAnalysis.jl b/src/Symmetry/SymmetryAnalysis.jl index 6802176d2..d00f89e75 100644 --- a/src/Symmetry/SymmetryAnalysis.jl +++ b/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) diff --git a/test/test_symmetry.jl b/test/test_symmetry.jl index 7f6b41591..b4ae62566 100644 --- a/test/test_symmetry.jl +++ b/test/test_symmetry.jl @@ -350,7 +350,7 @@ 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 @@ -358,7 +358,7 @@ end 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 @@ -366,21 +366,21 @@ end 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]