diff --git a/Project.toml b/Project.toml index 80d8b3edf..5b67ccde1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Sunny" uuid = "2b4a2ac8-8f8b-43e8-abf4-3cb0c45e8736" authors = ["The Sunny team"] -version = "0.7.4" +version = "0.7.5" [deps] Brillouin = "23470ee3-d0df-4052-8b1a-8cbd6363e7f0" diff --git a/docs/src/versions.md b/docs/src/versions.md index 438ef5124..30ef33085 100644 --- a/docs/src/versions.md +++ b/docs/src/versions.md @@ -1,5 +1,13 @@ # Version History +## v0.7.5 +(In development) + +* [`view_crystal`](@ref) shows allowed quadratic anisotropy. + [`print_site`](@ref) accepts an optional reference atom `i_ref`, with default + of `i`. The optional reference bond `b_ref` of [`print_bond`](@ref) now + defaults to `b`. + ## v0.7.4 (Dec 6, 2024) diff --git a/ext/PlottingExt/ViewCrystal.jl b/ext/PlottingExt/ViewCrystal.jl index 7908e221e..617a758a0 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" @@ -60,14 +60,64 @@ function propagate_reference_bond_for_cell(cryst, b_ref) return reduce(vcat, found) end +function anisotropy_on_site(sys, i) + interactions = isnothing(sys) ? nothing : Sunny.interactions_homog(something(sys.origin, sys)) + onsite = interactions[i].onsite + if onsite isa Sunny.HermitianC64 + onsite = Sunny.StevensExpansion(Sunny.matrix_to_stevens_coefficients(onsite)) + end + return onsite :: Sunny.StevensExpansion +end + +# Get the quadratic anisotropy as a 3×3 exchange matrix for atom `i` in the +# chemical cell. +function quadratic_anisotropy(sys, i) + # Get certain Stevens expansion coefficients + (; c0, c2) = anisotropy_on_site(sys, i) -# Get the 3×3 exchange matrix for bond `b` -function exchange_on_bond(interactions, b) - isnothing(interactions) && return zero(Sunny.Mat3) + # Undo RCS renormalization for quadrupolar anisotropy for spin-s + if sys.mode == :dipole + s = (sys.Ns[i] - 1) / 2 + c2 = c2 / Sunny.rcs_factors(s)[2] # Don't mutate c2 in-place! + end + + # Stevens quadrupole operators expressed as 3×3 spin bilinears + quadrupole_basis = [ + [1 0 0; 0 -1 0; 0 0 0], # 𝒪₂₂ = SˣSˣ - SʸSʸ + [0 0 1; 0 0 0; 1 0 0] / 2, # 𝒪₂₁ = (SˣSᶻ + SᶻSˣ)/2 + [-1 0 0; 0 -1 0; 0 0 2], # 𝒪₂₀ = 2SᶻSᶻ - SˣSˣ - SʸSʸ + [0 0 0; 0 0 1; 0 1 0] / 2, # 𝒪₂₋₁ = (SʸSᶻ + SᶻSʸ)/2 + [0 1 0; 1 0 0; 0 0 0], # 𝒪₂₋₂ = SˣSʸ + SʸSˣ + ] + + # The c0 coefficient incorporates a factor of S². For quantum spin + # operators, S² = s(s+1) I. For the large-s classical limit, S² = s² is a + # scalar. + S² = if sys.mode == :dipole_uncorrected + # Undoes extraction in `operator_to_stevens_coefficients`. Note that + # spin magnitude s² is set to κ², as originates from `onsite_coupling` + # for p::AbstractPolynomialLike. + sys.κs[i]^2 + else + # Undoes extraction in `matrix_to_stevens_coefficients` where 𝒪₀₀ = I. + s = (sys.Ns[i]-1) / 2 + s * (s+1) + end + + return c2' * quadrupole_basis + only(c0) * I / S² +end + +function coupling_on_bond(interactions, b) + isnothing(interactions) && return zero(Mat3) pairs = interactions[b.i].pair indices = findall(pc -> pc.bond == b, pairs) - isempty(indices) && return zero(Sunny.Mat3) - return pairs[only(indices)].bilin * Mat3(I) + return isempty(indices) ? nothing : pairs[only(indices)] +end + +# Get the 3×3 exchange matrix for bond `b` +function exchange_on_bond(interactions, b) + coupling = coupling_on_bond(interactions, b) + return isnothing(coupling) ? zero(Mat3) : coupling.bilin * Mat3(I) end # Get largest exchange interaction scale. For symmetric part, this is the @@ -100,7 +150,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) @@ -178,7 +228,7 @@ function draw_exchange_geometries(; ax, obs, ionradius, pts, scaled_exchanges) end function draw_bonds(; ax, obs, ionradius, exchange_mag, cryst, interactions, bonds, refbonds, color) - + # Map each bond to line segments in global coordinates segments = map(bonds) do b (; ri, rj) = Sunny.BondPos(cryst, b) @@ -226,6 +276,10 @@ function draw_bonds(; ax, obs, ionradius, exchange_mag, cryst, interactions, bon dmvecstr = join(Sunny.number_to_simple_string.(dmvec; digits=3), ", ") J_matrix_str *= "\nDM: [$dmvecstr]" end + c = coupling_on_bond(interactions, b) + if !isnothing(c) && (!iszero(c.biquad) || !isempty(c.general.data)) + J_matrix_str *= "\n + higher order terms" + end end return """ @@ -274,20 +328,33 @@ function is_type_degenerate(cryst, i) end # Construct atom labels for use in DataInspector -function label_atoms(cryst; ismagnetic) +function label_atoms(cryst; ismagnetic, sys) return map(1:natoms(cryst)) do i typ = cryst.types[i] rstr = Sunny.fractional_vec3_to_string(cryst.positions[i]) ret = [] - if ismagnetic && is_type_degenerate(cryst, i) - c = cryst.classes[i] - push!(ret, isempty(typ) ? "Class $c at $rstr" : "'$typ' (class $c) at $rstr") - else - push!(ret, isempty(typ) ? "Position $rstr" : "'$typ' at $rstr") - end + (; multiplicity, letter) = Sunny.get_wyckoff(cryst, i) + wyckstr = "$multiplicity$letter" + typstr = isempty(typ) ? "" : "'$typ', " + push!(ret, typstr * "Wyckoff $wyckstr, $rstr") + if ismagnetic - # TODO: Show onsite couplings? + if isnothing(sys) + # 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, prefix="Aniso: ", digits=8, atol=1e-12)) + else + aniso = quadratic_anisotropy(sys, i) + basis_strs = Sunny.number_to_simple_string.(aniso; digits=3) + push!(ret, Sunny.formatted_matrix(basis_strs; prefix="Aniso: ")) + (; c4, c6) = anisotropy_on_site(sys, i) + if !iszero(c4) || !iszero(c6) + push!(ret, " + higher order terms") + end + end end join(ret, "\n") end @@ -336,10 +403,10 @@ function draw_atoms_or_dipoles(; ax, full_crystal_toggle, dipole_menu, cryst, sy # Labels for non-ghost atoms inspector_label = nothing if !isghost - labels = label_atoms(xtal; ismagnetic)[idxs] + labels = label_atoms(xtal; ismagnetic, sys)[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/Operators/Rotation.jl b/src/Operators/Rotation.jl index d54669c9c..f7a259685 100644 --- a/src/Operators/Rotation.jl +++ b/src/Operators/Rotation.jl @@ -143,8 +143,10 @@ end Rotates the local quantum operator `A` according to the ``3×3`` rotation matrix `R`. """ -function rotate_operator(A::HermitianC64, R) +function rotate_operator(A::AbstractMatrix{ComplexF64}, R) isempty(A) && return A + A ≈ A' || error("Non-Hermitian operator") + A = hermitianpart(A) R = convert(Mat3, R) N = size(A, 1) U = unitary_irrep_for_rotation(R; N) diff --git a/src/Operators/Stevens.jl b/src/Operators/Stevens.jl index fb0bb0c6b..01ed65c4e 100644 --- a/src/Operators/Stevens.jl +++ b/src/Operators/Stevens.jl @@ -72,7 +72,8 @@ function stevens_abstract_polynomials(; J, k::Int) end -# Construct Stevens operators as polynomials in the spin operators. +# Construct Stevens operators as polynomials in the spin operators. Listed in +# descending order q = k,..-k. function stevens_matrices_of_dim(k::Int; N::Int) if k >= N return fill(Hermitian(zeros(ComplexF64, N, N)), 2k+1) @@ -134,6 +135,9 @@ end const stevens_αinv = map(inv, stevens_α) +# Expands matrix A in Stevens operators. The coefficients are returned as an +# OffsetArray c[k] with indices k = 0..6. Elements of c[k][:] are the Stevens +# coefficients in descending order q = k..-k. function matrix_to_stevens_coefficients(A::HermitianC64) N = size(A,1) @assert N == size(A,2) @@ -162,15 +166,15 @@ function operator_for_stevens_rotation(k, R) return real(V) end -# Let c denote coefficients of an operator expansion 𝒜 = c† 𝒪. Under the +# Let c denote coefficients of an operator expansion 𝒜 = cᵀ 𝒪. Under the # rotation R, Stevens operators transform as 𝒪 → V 𝒪. Alternatively, we can # treat the Stevens operators as fixed, provided the coefficients transform as -# c† → c† V, or c → V† c. -function rotate_stevens_coefficients(c, R::Mat3) +# cᵀ → cᵀ V, or c → Vᵀ c. +function rotate_stevens_coefficients(c::AbstractVector{Float64}, R::Mat3) N = length(c) k = Int((N-1)/2) V = operator_for_stevens_rotation(k, R) - return V' * c + return transpose(V) * c end diff --git a/src/Operators/Symbolic.jl b/src/Operators/Symbolic.jl index def550803..26389be2c 100644 --- a/src/Operators/Symbolic.jl +++ b/src/Operators/Symbolic.jl @@ -91,9 +91,9 @@ end # Extract Stevens operator coefficients from spin polynomial -function operator_to_stevens_coefficients(p::DP.AbstractPolynomialLike, S) +function operator_to_stevens_coefficients(p::DP.AbstractPolynomialLike, S²) p = expand_in_stevens_operators(p) - p = DP.subs(p, spin_squared_symbol => S^2) + p = DP.subs(p, spin_squared_symbol => S²) return map(stevens_symbols) do 𝒪ₖ map(𝒪ₖ) do 𝒪kq j = findfirst(==(𝒪kq), DP.monomials(p)) @@ -102,7 +102,7 @@ function operator_to_stevens_coefficients(p::DP.AbstractPolynomialLike, S) end end -function rotate_operator(p, R) +function rotate_operator(p::DP.AbstractPolynomialLike, R) 𝒮′ = R * [𝒮ˣ, 𝒮ʸ, 𝒮ᶻ] DP.subs(p, 𝒮ˣ => 𝒮′[1], 𝒮ʸ => 𝒮′[2], 𝒮ᶻ => 𝒮′[3]) end diff --git a/src/Symmetry/AllowedAnisotropy.jl b/src/Symmetry/AllowedAnisotropy.jl index 82f0cf61b..79dc80c0f 100644 --- a/src/Symmetry/AllowedAnisotropy.jl +++ b/src/Symmetry/AllowedAnisotropy.jl @@ -16,30 +16,38 @@ function transform_spherical_to_stevens_coefficients(k, c) return transpose(stevens_αinv[k]) * c end +# The bare Stevens expansion coefficients b::Vector{Float64} are not understood +# by rotate_operator(). Embed them in a large StevensExpansion object that can +# be passed to is_anisotropy_valid(). +function is_stevens_expansion_valid(cryst, i, b) + N = length(b) + c0 = (N == 1) ? b : zeros(1) + c2 = (N == 5) ? b : zeros(5) + c4 = (N == 9) ? b : zeros(9) + c6 = (N == 13) ? b : zeros(13) + stvexp = StevensExpansion(6, c0, c2, c4, c6) + return is_anisotropy_valid(cryst, i, stvexp) +end -function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int, R::Mat3) - # The symmetry operations for the point group at atom i. Each one encodes a - # rotation/reflection. +function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int, R_global::Mat3, atol::Float64) + # The symmetry operations for the point group at atom i. symops = symmetries_for_pointgroup_of_atom(cryst, i) # The Wigner D matrices for each symop Ds = map(symops) do s # R is an orthogonal matrix that transforms positions, x → x′ = R x. It - # might or might not include a reflection, i.e., det R = ±1. - sR = cryst.latvecs * s.R * inv(cryst.latvecs) - - # Unlike position x, spin S = [Sx, Sy, Sz] is a _pseudo_ vector, which - # means that, under reflection, the output gains an additional minus - # sign. That is, the orthogonal transformation R applied to spin has the - # action, S → S′ = ± R S, where the minus sign corresponds to the case - # det(R) = -1. More simply, we may write S′ = Q S, where Q = det(R) R. - Q = det(sR) * sR - - # The Wigner D matrix, whose action on a spherical tensor corresponds to - # the 3x3 rotation Q (see more below). + # might or might not include an inversion, i.e., det R = ±1. + R = cryst.latvecs * s.R * inv(cryst.latvecs) + + # Unlike position, spin angular momentum is a pseudo-vector, which means + # it is invariant under global inversion. The transformation of spin is + # always a pure rotation, S → S′ = Q S, where Q = R det R. + Q = det(R) * R + + # The Wigner D matrix corresponding to Q (see more below). return unitary_irrep_for_rotation(Q; N=2k+1) end - + # A general operator in the spin-k representation can be decomposed in the # basis of spherical tensors, 𝒜 = ∑_q c_q T_kq, for some coefficients c_q. # Spherical tensors transform as T_kq → D^{*}_qq′ T_kq′. Alternatively, we @@ -53,14 +61,6 @@ function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int, # operator 𝒜. C = sum(D' for D in Ds) - # Transform coefficients c to c′ in rotated Stevens operators, T′ = D* T, - # where the Wigner D matrix is associated with the rotation R. That is, find - # c′ satisfying c′ᵀ T′ = c T. Recall c′ᵀ T′ = (c′ᵀ D*) T = (D† c′)ᵀ T. The - # constraint becomes D† c′ = c. Since D is unitary, we have c′ = D c. We - # apply this transformation to each column c of C. - D = unitary_irrep_for_rotation(R; N=2k+1) - C = D * C - # 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 𝒪. @@ -68,16 +68,49 @@ function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int, transform_spherical_to_stevens_coefficients(k, c) end + # Apply a global rotation to the Cartesian coordinate system. Stevens + # operators rotate as 𝒪′ = V 𝒪. Coefficients satisfying b′ᵀ 𝒪′ = bᵀ 𝒪 + # must therefore transform as b′ = V⁻ᵀ b. + V = operator_for_stevens_rotation(k, R_global) + B = transpose(V) \ B + # 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) + B = colspace(hcat(real(B), imag(B)); atol) # Find linear combination of columns that sparsifies B - return sparsify_columns(B; atol=1e-12) + B = sparsify_columns(B; atol) + + # Scale each column to make the final expression prettier + return map(eachcol(B)) do b + b /= argmax(abs, 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 + + # Magnify 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 + + # Check that the operator bᵀ 𝒪 satisifies point group symmetries in the + # original coordinate system. Multiply by Vᵀ to effectively undo the + # global rotation of Stevens operators 𝒪. + @assert is_stevens_expansion_valid(cryst, i, transpose(V) * b) + + b + end + + return B end diff --git a/src/Symmetry/AllowedCouplings.jl b/src/Symmetry/AllowedCouplings.jl index 73ef839f6..8385be6ad 100644 --- a/src/Symmetry/AllowedCouplings.jl +++ b/src/Symmetry/AllowedCouplings.jl @@ -222,15 +222,28 @@ end # Returns a list of ``3×3`` matrices that form a linear basis for the # symmetry-allowed coupling matrices associated with bond `b`. -function basis_for_symmetry_allowed_couplings_aux(cryst::Crystal, b::BondPos) +function basis_for_symmetry_allowed_couplings_aux(cryst::Crystal, b::BondPos; R_global::Mat3) # Expected floating point precision for 9x9 matrix operations atol = 1e-12 + # Output basis vectors x should be reported in a rotated Cartesian + # coordinate system prior to sparsification. The transformation x → K x is + # equivalent to J → R J Rᵀ, where J = reshape(x, 3, 3). Note that K† = K⁻¹. + K = kron(R_global, R_global) + + # P is a product of projection operators, each of which impose one + # constraint. An eigenvector with eigenvalue of 1 satisfies all constraints, + # and is therefore an allowed coupling. P = symmetry_allowed_couplings_operator(cryst, b) - # Any solution to the original symmetry constraints `R J Rᵀ = J` or `R J Rᵀ - # = Jᵀ` decomposes into purely symmetric/antisymmetric solutions. Therefore - # we can pick a basis that separates into symmetric and antisymmetric parts. - # We will do so by decomposing P. By construction, P = P_sym+P_asym. + + # Global rotation of the coordinate system. Solutions to K P K⁻¹ x̃ = x̃ + # yield solutions to P x = x where x̃ = K x. + P = K * P * K' + + # Each constraint has the form `R J Rᵀ = J` or `R J Rᵀ = Jᵀ`. If J is a + # solution, then its symmetric and antisymmetric parts are individually + # solutions. Knowing this, decompose P to solve for the symmetric and + # antisymmetric parts separately. By construction, P = P_sym+P_asym. P_sym = P * sym_basis * sym_basis' P_asym = P * asym_basis * asym_basis' @@ -238,21 +251,21 @@ function basis_for_symmetry_allowed_couplings_aux(cryst::Crystal, b::BondPos) acc_asym = Vector{Float64}[] # If any "reference" basis vectors are eigenvalues of P_sym with eigenvalue - # 1, use them as outputs, and remove them from P_sym + # 1, use them as outputs, and remove them from P_sym. for x in eachcol(sym_basis) if isapprox(P_sym*x, x; atol) push!(acc_sym, x) P_sym = P_sym * (I - x*x') end end - # Same for P_asym + # Same for P_asym. for x in eachcol(asym_basis) if isapprox(P_asym*x, x; atol) push!(acc_asym, x) P_asym = P_asym * (I - x*x') end end - + # Search for eigenvectors of P_sym with eigenvalue 1. These provide an # orthonormal basis for symmetric couplings. v = nullspace(P_sym-I; atol) @@ -272,15 +285,15 @@ function basis_for_symmetry_allowed_couplings_aux(cryst::Crystal, b::BondPos) return map(acc) do x # Normalize each basis vector so that its maximum component is 1. The # shift by atol avoids unnecessary sign change in the case where the - # maximum magnitude values of x appear symmetrically as ±c for some c. - _, i = findmax(abs.(x.+atol)) - x = x / x[i] + # maximum magnitude values of x appear symmetrically as ±c. + x /= argmax(c -> abs(c + atol), x) - # Reinterpret as 3x3 matrix + # Reinterpret as 3x3 matrix. x = Mat3(reshape(x, 3, 3)) - - # Double check that x indeed satifies the necessary symmetries - @assert is_coupling_valid(cryst, b, x) + + # Check that the coupling J satisifies the point group symmetries in the + # original coordinate system (after undoing the R_global rotation). + @assert is_coupling_valid(cryst, b, R_global' * x * R_global) return x end @@ -297,8 +310,8 @@ function transform_coupling_for_bonds(cryst, b, b_ref, J_ref) return transform_coupling_by_symmetry(J_ref, R*det(R), parity) end -function basis_for_symmetry_allowed_couplings(cryst::Crystal, b::Bond; b_ref=b) - basis = basis_for_symmetry_allowed_couplings_aux(cryst, BondPos(cryst, b_ref)) +function basis_for_symmetry_allowed_couplings(cryst::Crystal, b::Bond; b_ref=b, R_global=Mat3(I)) + basis = basis_for_symmetry_allowed_couplings_aux(cryst, BondPos(cryst, b_ref); R_global) # Transform coupling basis from `b_ref` to `b` if b == b_ref 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 20cbdb608..5ad3e0115 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,22 +109,35 @@ 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, bond::Bond; b_ref::Bond) + print_bond(cryst::Crystal, b::Bond; b_ref=b) -Prints symmetry information for bond `bond`. A symmetry-equivalent reference -bond `b_ref` can optionally be provided to fix the meaning of the coefficients -`A`, `B`, ... +Prints symmetry information for bond `b`. An optional symmetry-equivalent +reference bond `b_ref` can be provided to keep a consistent meaning of the free +parameters `A`, `B`, etc. """ -function print_bond(cryst::Crystal, b::Bond; b_ref=nothing, io=stdout) +function print_bond(cryst::Crystal, b::Bond; b_ref=b, io=stdout) # How many digits to use in printing coefficients digits = 14 # Tolerance below which coefficients are dropped atol = 1e-12 - + if b.i == b.j && iszero(b.n) - print_site(cryst, b.i; io) + print_site(cryst, b.i; i_ref=b.i, io) else ri = cryst.positions[b.i] rj = cryst.positions[b.j] + b.n @@ -133,7 +145,7 @@ function print_bond(cryst::Crystal, b::Bond; b_ref=nothing, 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 @@ -145,12 +157,6 @@ function print_bond(cryst::Crystal, b::Bond; b_ref=nothing, io=stdout) println(io, "Connects '$(cryst.types[b.i])' at $(fractional_vec3_to_string(ri)) to '$(cryst.types[b.j])' at $(fractional_vec3_to_string(rj))") end - # If `b_ref` is nothing, select it from `reference_bonds` - b_ref = @something b_ref begin - d = global_distance(cryst, b) - ref_bonds = reference_bonds(cryst, d; min_dist=d) - only(filter(b′ -> is_related_by_symmetry(cryst, b, b′), ref_bonds)) - end basis = basis_for_symmetry_allowed_couplings(cryst, b; b_ref) basis_strs = coupling_basis_strings(zip('A':'Z', basis); digits, atol) println(io, formatted_matrix(basis_strs; prefix="Allowed exchange matrix: ")) @@ -184,7 +190,7 @@ b)` for every bond `b` in `reference_bonds(cryst, max_dist)`, where function print_symmetry_table(cryst::Crystal, max_dist; io=stdout) validate_crystal(cryst) for b in reference_bonds(cryst, max_dist) - print_bond(cryst, b; b_ref=b, io) + print_bond(cryst, b; io) end end @@ -198,89 +204,84 @@ 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 """ - print_site(cryst, i; R=I) + print_site(cryst, i; i_ref=i, R=I) Print symmetry information for the site `i`, including allowed g-tensor and -allowed anisotropy operator. An optional rotation matrix `R` can be provided to -define the reference frame for expression of the anisotropy. +allowed anisotropy operator. An optional symmetry-equivalent reference atom +`i_ref` can be provided to keep a consistent meaning of the free parameters. An +optional rotation matrix `R` can map to a new Cartesian reference frame for +expression of the allowed anisotropy. """ -function print_site(cryst, i; R=Mat3(I), ks=[2,4,6], io=stdout) +function print_site(cryst, i; i_ref=i, R=Mat3(I), ks=[2,4,6], io=stdout) + R_global = convert(Mat3, R) r = cryst.positions[i] class_i = cryst.classes[i] - m = count(==(class_i), cryst.classes) printstyled(io, "Atom $i\n"; bold=true, color=:underline) + (; multiplicity, letter) = get_wyckoff(cryst, i) + wyckstr = "$multiplicity$letter" + if isempty(cryst.types[i]) - println(io, "Position $(fractional_vec3_to_string(r)), multiplicity $m") + println(io, "Position $(fractional_vec3_to_string(r)), Wyckoff $wyckstr") else - println(io, "Type '$(cryst.types[i])', position $(fractional_vec3_to_string(r)), multiplicity $m") + println(io, "Type '$(cryst.types[i])', position $(fractional_vec3_to_string(r)), Wyckoff $wyckstr") end - # Tolerance below which coefficients are dropped + digits = 14 atol = 1e-12 - # How many digits to use in printing coefficients - 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: `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: ")) - - print_allowed_anisotropy(cryst, i; R, atol, digits, ks, io) + 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 int_to_underscore_string(x::Int) - subscripts = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'] - chars = collect(repr(x)) - if chars[begin] == '-' - popfirst!(chars) - sign = "-" + +function rotation_between_sites(cryst, i, i_ref) + # Rotation that maps from i_ref to i + if i == i_ref + return Mat3(I) else - sign = "" + syms = symmetries_between_atoms(cryst, i, i_ref) + isempty(syms) && error("Atoms $i and $i_ref are not symmetry equivalent.") + R = cryst.latvecs * first(syms).R * inv(cryst.latvecs) + return R * det(R) # Remove possible inversion (appropriate for spin pseudo-vector) end - digits = map(c -> parse(Int, c), chars) - return sign * prod(subscripts[digits.+1]) end +function allowed_g_tensor_string(cryst, i_ref; R_global=Mat3(I), R_site, prefix="Allowed g-tensor: ", 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) + return formatted_matrix(basis_strs; prefix) +end -function print_allowed_anisotropy(cryst::Crystal, i::Int; R::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[] cnt = 1 for k in ks - B = basis_for_symmetry_allowed_anisotropies(cryst, i; k, R) + B = basis_for_symmetry_allowed_anisotropies(cryst, i_ref; k, R_global, atol) + + # B is an allowed basis for i_ref, but we want to print the allowed + # basis for i. These sites are symmetry equivalent under the rotation + # R_site. V is the corresponding linear operator that acts on Stevens + # operators, 𝒪′ = V 𝒪. Coefficients satisfying b′ᵀ 𝒪′ = bᵀ 𝒪 then + # transform as b′ = V⁻ᵀ b. + V = operator_for_stevens_rotation(k, R_site) + B = [transpose(V) \ b for b in B] 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 -> 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 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 - + for b in reverse(B) # 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) @@ -291,7 +292,11 @@ function print_allowed_anisotropy(cryst::Crystal, i::Int; R::Mat3, atol, digits, end # clean up printing of term - ops = length(ops) == 1 ? ops[1] : "("*join(ops, "+")*")" + ops = if length(ops) > 1 || startswith(only(ops), '-') + "("*join(ops, "+")*")" + else + only(ops) + end ops = replace(ops, "+-" => "-") push!(terms, "c" * int_to_underscore_string(cnt) * "*" * ops) cnt += 1 @@ -299,10 +304,10 @@ function print_allowed_anisotropy(cryst::Crystal, i::Int; R::Mat3, atol, digits, push!(lines, prefix * join(terms, " + ")) end end - println(io, "Allowed anisotropy in Stevens operators:") - println(io, join(lines, " +\n")) - if R != I - println(io, "Modified reference frame! Transform using `rotate_operator(op, R)`.") + ret = "Allowed anisotropy in Stevens operators:\n" * join(lines, " +\n") + if R_global != I + 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/src/System/OnsiteCoupling.jl b/src/System/OnsiteCoupling.jl index 6327e637e..bf42fee85 100644 --- a/src/System/OnsiteCoupling.jl +++ b/src/System/OnsiteCoupling.jl @@ -24,8 +24,8 @@ function onsite_coupling(sys, site, p::DP.AbstractPolynomialLike) error("Symbolic operator only valid for system with mode `:dipole_uncorrected`.") end - s = sys.κs[site] - c = operator_to_stevens_coefficients(p, s) + S² = sys.κs[site]^2 + c = operator_to_stevens_coefficients(p, S²) # No renormalization here because symbolic polynomials `p` are associated # with the large-s limit. @@ -148,7 +148,7 @@ function set_onsite_coupling!(sys::System, op, i::Int) if !is_anisotropy_valid(sys.crystal, i, onsite) error("""Symmetry-violating anisotropy: $op. - Use `print_site(crystal, $i)` for more information.""") + Use `print_site(cryst, $i)` for more information.""") end cryst = sys.crystal diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl index 8bc27a7fa..125392cf2 100644 --- a/src/System/PairExchange.jl +++ b/src/System/PairExchange.jl @@ -220,17 +220,17 @@ function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float # Verify that couplings are symmetry-consistent if !is_coupling_valid(sys.crystal, bond, bilin) error("""Symmetry-violating bilinear exchange $bilin. - Use `print_bond(crystal, $bond)` for more information.""") + Use `print_bond(cryst, $bond)` for more information.""") end if !is_coupling_valid(sys.crystal, bond, biquad) biquad_str = formatted_matrix(number_to_math_string.(biquad); prefix=" ") error("""Symmetry-violating biquadratic exchange (written in Stevens basis) $biquad_str - Use `print_bond(crystal, $bond)` for more information.""") + Use `print_bond(cryst, $bond)` for more information.""") end if !is_coupling_valid(sys.crystal, bond, tensordec) error("""Symmetry-violating coupling. - Use `print_bond(crystal, $bond)` for more information.""") + Use `print_bond(cryst, $bond)` for more information.""") end # Print a warning if an interaction already exists for bond diff --git a/src/System/Types.jl b/src/System/Types.jl index cdace1b23..54af9c4ac 100644 --- a/src/System/Types.jl +++ b/src/System/Types.jl @@ -1,6 +1,8 @@ -# Stevens function expansion, renormalized for dipole projection +# Coefficients for an expansion in Stevens functions. Each component ck lists +# coefficients in descending order q = k,..-k. In :dipole mode, the +# renormalization factors (`rcs_factor`) are included in these coefficients. struct StevensExpansion - kmax ::Int + kmax :: Int c0 :: SVector{1, Float64} c2 :: SVector{5, Float64} c4 :: SVector{9, Float64} diff --git a/test/test_operators.jl b/test/test_operators.jl index d9f161b10..d1b32d03d 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -330,11 +330,11 @@ end # Test Stevens coefficients extraction S = spin_matrices(Inf) O = stevens_matrices(Inf) - S_mag = π + S² = π p = S'*S * O[4, 2] - c = Sunny.operator_to_stevens_coefficients(p, S_mag) + c = Sunny.operator_to_stevens_coefficients(p, S²) @test iszero(c[1]) && iszero(c[2]) && iszero(c[3]) && iszero(c[5]) && iszero(c[6]) - @test c[4] ≈ [0.0, 0.0, S_mag^2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + @test c[4] ≈ [0.0, 0.0, S², 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] # Test round trip Stevens -> spin -> Stevens c_ref = map(OffsetArrays.OffsetArray(0:6, 0:6)) do k diff --git a/test/test_symmetry.jl b/test/test_symmetry.jl index 21d50e4aa..00d262210 100644 --- a/test/test_symmetry.jl +++ b/test/test_symmetry.jl @@ -340,7 +340,7 @@ end end @test capt.output == """ Atom 1 - Position [0, 0, 0], multiplicity 16 + Position [0, 0, 0], Wyckoff 16c Allowed g-tensor: [A B B B A B B B A] @@ -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 @@ -389,6 +389,24 @@ end """ + capt = IOCapture.capture() do + print_site(cryst, 5; i_ref=2) + end + @test capt.output == """ + Atom 5 + Position [1/4, 0, 1/4], Wyckoff 16c + Allowed g-tensor: [ A -B B + -B A -B + 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𝒪[6,-6]+24𝒪[6,-5]+5𝒪[6,-2]+8𝒪[6,-1]-8𝒪[6,1]-24𝒪[6,5]) + """ + + 𝒪 = stevens_matrices(4) + @test Sunny.is_anisotropy_valid(cryst, 5, 7𝒪[4,-3]+2𝒪[4,-2]-𝒪[4,-1]+𝒪[4,1]+7𝒪[4,3]) + capt = IOCapture.capture() do print_suggested_frame(cryst, 2) end @@ -404,15 +422,15 @@ end end @test capt.output == """ Atom 2 - Position [1/4, 1/4, 0], multiplicity 16 - Allowed g-tensor: [A-B 0 0 - 0 A-B 0 - 0 0 A+2B] + Position [1/4, 1/4, 0], Wyckoff 16c + Allowed g-tensor: [A 0 0 + 0 A 0 + 0 0 B] 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)`. + Modified reference frame! Use R*g*R' or rotate_operator(op, R). """ cryst = Sunny.hyperkagome_crystal() @@ -433,44 +451,50 @@ end 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"]) + 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] capt = IOCapture.capture() do - print_site(crystal, 1; R) - print_site(crystal, 2; R) + print_site(cryst, 1; R) + print_site(cryst, 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] + Type 'A', position [0, 0, 0], Wyckoff 3a + Allowed g-tensor: [A 0 0 + 0 A 0 + 0 0 B] 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)`. + Modified reference frame! Use R*g*R' or 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] + Type 'B', position [1/2, 1/2, 0], Wyckoff 9e + Allowed g-tensor: [A 0 0 + 0 B D+E + 0 D-E C] 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)`. + Modified reference frame! Use R*g*R' or 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) + @test Sunny.is_anisotropy_valid(cryst, 2, 𝒪[6,-1]+0.997385420𝒪[6,1]) + @test Sunny.is_anisotropy_valid(cryst, 2, rotate_operator(𝒪[6,2], R)) end