Library API

(cell1, cell2, cell3, i) :: Site

Four indices identifying a single site in a System. The first three indices select the unit cell and the last index selects the sublattice, i.e., the $i$th atom within the unit cell.

This object can be used to index dipoles and coherents fields of a System. A Site is also required to specify inhomogeneous interactions via functions such as set_field_at! or set_exchange_at!.

Note that the definition of a cell may change when a system is reshaped. In this case, it is convenient to construct the Site using position_to_site, which always takes a position in fractional coordinates of the original lattice vectors.

BinningParameters(binstart,binend,binwidth;covectors = I(4))
BinningParameters(binstart,binend;numbins,covectors = I(4))

Describes a 4D parallelepided histogram in a format compatible with experimental Inelasitic Neutron Scattering data. See generate_mantid_script_from_binning_parameters to convert BinningParameters to a format understandable by the Mantid software, or load_nxs to load BinningParameters from a Mantid .nxs file.

The coordinates of the histogram axes are specified by multiplication of (q,ω) with each row of the covectors matrix, with q given in [R.L.U.]. Since the default covectors matrix is the identity matrix, the default axes are (qx,qy,qz,ω) in absolute units.

The convention for the binning scheme is that:

• The left edge of the first bin starts at binstart
• The bin width is binwidth
• The last bin contains binend
• There are no "partial bins;" the last bin may contain values greater than binend.

A value can be binned by computing its bin index:

coords = covectors * value
bin_ix = 1 .+ floor.(Int64,(coords .- binstart) ./ binwidth)

Bond(i, j, n)

Represents a bond between atom indices i and j. n is a vector of three integers specifying unit cell displacement in terms of lattice vectors.

An object describing a crystallographic unit cell and its space group symmetry. Constructors are as follows:

Crystal(filename; override_symmetry=false, symprec=nothing)

Reads the crystal from a .cif file located at the path filename. If override_symmetry=true, the spacegroup will be inferred based on atom positions and the returned unit cell may be reduced in size. For an mCIF file, the return value is the magnetic supercell, unless override_symmetry=true. If a precision for spacegroup symmetries cannot be inferred from the CIF file, it must be specified with symprec. The latvecs field of the returned Crystal will be in units of angstrom.

Crystal(latvecs, positions; types=nothing, symprec=1e-5)

Constructs a crystal from the complete list of atom positions positions, with coordinates (between 0 and 1) in units of lattice vectors latvecs. Spacegroup symmetry information is automatically inferred. The optional parameter types is a list of strings, one for each atom, and can be used to break symmetry-equivalence between atoms.

Crystal(latvecs, positions, spacegroup_number; types=nothing, setting=nothing, symprec=1e-5)

Builds a crystal by applying symmetry operators for a given international spacegroup number. For certain spacegroups, there are multiple possible unit cell settings; in this case, a warning message will be printed, and a list of crystals will be returned, one for every possible setting. Alternatively, the optional setting string will disambiguate between unit cell conventions.

Currently, crystals built using only the spacegroup number will be missing some symmetry information. It is generally preferred to build a crystal from a .cif file or from the full specification of the unit cell.

Examples

# Read a Crystal from a .cif file
Crystal("filename.cif")

# Build a BCC crystal in the conventional cubic unit cell by specifying both
# atoms. The spacegroup 229 is inferred.
latvecs = lattice_vectors(1, 1, 1, 90, 90, 90)
positions = [[0, 0, 0], [1/2, 1/2, 1/2]]
Crystal(latvecs, positions)

# Build a CsCl crystal (two simple cubic sublattices). Because of the distinct
# atom types, the spacegroup number 221 is now inferred.
types = ["Na", "Cl"]
cryst = Crystal(latvecs, positions; types)

# Build a diamond cubic crystal from its spacegroup number 227 and a single
# atom position. This spacegroup has two possible settings ("1" or "2"), which
# determine an overall unit cell translation.
positions = [[1/4, 1/4, 1/4]]
cryst = Crystal(latvecs, positions, 227; setting="1")

See also lattice_vectors.

FormFactor(ion::String; g_lande=2)

The magnetic form factor for a given magnetic ion and charge state. When passed to intensities, it rescales structure factor intensities based on the magnitude of the scattering vector, $|𝐪|$.

The parameter ion must be one of the following strings:

Am2, Am3, Am4, Am5, Am6, Am7, Au1, Au2, Au3, Au4, Au5, Ce2, Co0, Co1, Co2, Co3,
Co4, Cr0, Cr1, Cr2, Cr3, Cr4, Cu0, Cu1, Cu2, Cu3, Cu4, Dy2, Dy3, Er2, Er3, Eu2,
Eu3, Fe0, Fe1, Fe2, Fe3, Fe4, Gd2, Gd3, Hf2, Hf3, Ho2, Ho3, Ir0a, Ir0b, Ir0c,
Ir1a, Ir1b, Ir2, Ir3, Ir4, Ir5, Ir6, Mn0, Mn1, Mn2, Mn3, Mn4, Mn5, Mo0, Mo1, Nb0,
Nb1, Nd2, Nd3, Ni0, Ni1, Ni2, Ni3, Ni4, Np3, Np4, Np5, Np6, Os0a, Os0b, Os0c,
Os1a, Os1b, Os2, Os3, Os4, Os5, Os6, Os7, Pd0, Pd1, Pr3, Pt1, Pt2, Pt3, Pt4,
Pt5, Pt6, Pu3, Pu4, Pu5, Pu6, Re0a, Re0b, Re0c, Re1a, Re1b, Re2, Re3, Re4, Re5,
Re6, Rh0, Rh1, Ru0, Ru1, Sc0, Sc1, Sc2, Sm2, Sm3, Ta2, Ta3, Ta4, Tb2, Tb3, Tc0,
Tc1, Ti0, Ti1, Ti2, Ti3, Tm2, Tm3, U3, U4, U5, V0, V1, V2, V3, V4, W0a, W0b,
W0c, W1a, W1b, W2c, W3, W4, W5, Y0, Yb2, Yb3, Zr0, Zr1

The trailing number denotes ionization state. For example, "Fe0" denotes a neutral iron atom, while "Fe2" denotes Fe²⁺. If multiple electronic configurations are possible, they will be distinguished by a trailing letter (a, b, ...). Omitting this letter will print an informative error,

FormFactor("Ir0")

ERROR: Disambiguate form factor according to electronic configuration:
    "Ir0a" -- 6s⁰5d⁹
    "Ir0b" -- 6s¹5d⁸
    "Ir0c" -- 6s²5d⁷

In the dipolar approximation (small $|𝐪|$) the form factor is

$F(s) = ⟨j_0(s)⟩ + [(2-g)/g] ⟨j_2(s)⟩$,

involving $s = |𝐪|/4π$ and the Landé $g$-factor. The $⟨j_l(s)⟩$ are radial averages of the $l$th spherical Bessel function of the magnetic dipole. More details are provided in Ref. [1].

The standard approximation tables involve expansion in Gaussians,

\[⟨j_0(s)⟩ = A e^{-as^2} + B e^{-bs^2} + C e^{-cs^2} + D e^{-ds^2} + E\]

and

\[⟨j_2(s)⟩ = (A e^{-as^2} + B e^{-bs^2} + C e^{-cs^2} + D e^{-ds^2} + E) s^2.\]

For 3d, 4d, rare earth, and actinide ions, Sunny uses the revised tables of P. J. Brown, as documented in the McPhase package [2]. For 5d ions, Sunny uses the tables of Kobayashi, Nagao, Ito [3].

Two special, $𝐪$-independent form factor values are available: one(FormFactor) and zero(FormFactor). The first idealizes the magnetic ion as a perfect point particle, while the second zeros all contributions from the magnetic ion.

References:

1. P. J. Brown, The Neutron Data Booklet, 2nd ed., Sec. 2.5 Magnetic Form Factors (2003)
2. Coefficient tables in McPhase documentation
3. K. Kobayashi, T. Nagao, M. Ito, Acta Cryst. A, 67 pp 473–480 (2011)

ImplicitMidpoint(dt::Float64; atol=1e-12) where N

The implicit midpoint method for integrating the Landau-Lifshitz spin dynamics or its generalization to SU(N) coherent states [1]. One call to the step! function will advance a System by dt units of time. This integration scheme is exactly symplectic and eliminates energy drift over arbitrarily long simulation trajectories.

References:

1. H. Zhang and C. D. Batista, Phys. Rev. B 104, 104409 (2021).

Langevin(dt::Float64; damping::Float64, kT::Float64)

An integrator for Langevin spin dynamics using the explicit Heun method. The damping parameter controls the coupling to an implicit thermal bath. One call to the step! function will advance a System by dt units of time. Can be used to sample from the Boltzmann distribution at temperature kT. An alternative approach to sampling states from thermal equilibrium is LocalSampler, which proposes local Monte Carlo moves. For example, use LocalSampler instead of Langevin to sample Ising-like spins.

Setting damping = 0 disables coupling to the thermal bath, yielding an energy-conserving spin dynamics. The Langevin integrator uses an explicit numerical integrator which cannot prevent energy drift. Alternatively, the ImplicitMidpoint method can be used, which is more expensive but prevents energy drift through exact conservation of the symplectic 2-form.

If the System has mode = :dipole, then the dynamics is the stochastic Landau-Lifshitz equation,

\[ d𝐬/dt = -𝐬 × (ξ - 𝐁 + λ 𝐬 × 𝐁),\]

where $𝐁 = -dE/d𝐬$ is the effective field felt by the expected spin dipole $𝐬$. The components of $ξ$ are Gaussian white noise, with magnitude $√(2 k_B T λ)$ set by a fluctuation-dissipation theorem. The parameter damping sets the phenomenological coupling $λ$ to the thermal bath.

If the System has mode = :SUN, then this dynamics generalizes [1] to a stochastic nonlinear Schrödinger equation for SU(N) coherent states $𝐙$,

\[ d𝐙/dt = -i P [ζ + (1 - i λ̃) ℋ 𝐙].\]

Here, $P$ projects onto the space orthogonal to $𝐙$, and $ζ$ denotes complex Gaussian white noise with magnitude $√(2 k_B T λ̃)$. The local-Hamiltonian $ℋ$ embeds the energy gradient into the 𝔰𝔲(N) Lie algebra, and generates evolution of spin dipoles, quadrupoles, etc. The parameter damping here sets $λ̃$, which is analogous to $λ$ above.

When applied to SU(2) coherent states, the generalized spin dynamics reduces exactly to the stochastic Landau-Lifshitz equation. The mapping is as follows. Normalized coherent states $𝐙$ map to dipole expectation values $𝐬 = 𝐙^{†} Ŝ 𝐙$, where spin operators $Ŝ$ are a spin-$|𝐬|$ representation of SU(2). The local effective Hamiltonian $ℋ = -𝐁 ⋅ Ŝ$ generates rotation of the dipole in analogy to the vector cross product $S × 𝐁$. The coupling to the thermal bath maps as $λ̃ = |𝐬| λ$. Note, therefore, that the scaling of the damping parameter varies subtly between :dipole and :SUN modes.

References:

1. D. Dahlbom et al., Phys. Rev. B 106, 235154 (2022).

LocalSampler(; kT, nsweeps=1.0, propose=propose_uniform)

Monte Carlo simulation involving Metropolis updates to individual spins. One call to the step! function will perform nsweeps of MCMC sampling for a provided System. The default value of 1.0 means that step! performs, on average, one trial update per spin.

Assuming ergodicity, the LocalSampler will sample from thermal equilibrium for the target temperature kT.

The trial spin updates are sampled using the propose function. Options include propose_uniform, propose_flip, and propose_delta. Multiple proposals can be mixed with the macro @mix_proposals.

The returned object stores fields ΔE and Δs, which represent the cumulative change to the net energy and dipole, respectively.

!!! warning "Efficiency considerations

A [`Langevin`](@ref) sampler is frequently much more efficient than a
`LocalSampler` for simulating Heisenberg-like spins that vary continuously. A
`LocalSampler` is appropriate in the special case that the spin states are
effectively discrete. E.g., [`propose_flip`](@ref) is very helpful simulating
Ising-like spins that arise due to a strong easy-axis anisotropy.

Moment(; s, g)

Characterizes a effective spin magnetic moment on an atom. Quantum spin-s is a multiple of 1/2 in units of ħ. The g-factor or tensor defines the magnetic_moment $μ = - g 𝐒$ in units of the Bohr magneton.

Example

Moment(s=3/2, g=2)

SampledCorrelations(sys::System; measure, energies, dt)

An object to accumulate samples of dynamical pair correlations. The measure argument specifies a pair correlation type, e.g. ssf_perp. The energies must be evenly-spaced and starting from 0, e.g. energies = range(0, 3, 100). Select the integration time-step dt according to accuracy and speed considerations. suggest_timestep can help in selecting an appropriate value.

Dynamical correlations will be accumulated through calls to add_sample!, which expects a spin configuration in thermal equilibrium. A classical spin dynamics trajectory will be simulated of sufficient length to achieve the target energy resolution. The resulting data can can then be extracted as pair-correlation intensities with appropriate classical-to-quantum correction factors. See also intensities_static, which integrates over energy.

SampledCorrelationsStatic(sys::System; measure)

An object to accumulate samples of static pair correlations. It is similar to SampledCorrelations, but no time-integration will be performed on calls to add_sample!. The resulting object can be used with intensities_static to calculate statistics from the classical Boltzmann distribution. Dynamical intensities data, however, will be unavailable. Similarly, classical-to-quantum corrections that rely on the excitation spectrum cannot be performed.

SpinWaveTheory(sys::System; measure, regularization=1e-8)

Constructs an object to perform linear spin wave theory. The system must be in an energy minimizing configuration. Enables calculation of dispersion bands. If pair correlations are specified with correspec, one can also calculate intensities_bands and broadened intensities.

The spins in system must be energy-minimized, otherwise the Cholesky step of the Bogoliubov diagonalization procedure will fail. The parameter regularization adds a small positive shift to the diagonal of the dynamical matrix to avoid numerical issues with quasi-particle modes of vanishing energy. Physically, this shift can be interpreted as application of an inhomogeneous field aligned with the magnetic ordering.

SpinWaveTheoryKPM(sys::System; measure, resolution, regularization=1e-8)

Experimental

An alternative to SpinWaveTheory that uses the kernel polynomial method (KPM) to perform intensities calculations. In traditional spin wave theory calculations, one would explicitly diagonalize the dynamical matrix, with a cost that scales like $𝒪(V^3)$ in the volume $V$ of the magnetic cell. KPM instead approximates intensities using polynomial expansion of the dynamical matrix. The computational cost of KPM scales like $𝒪(V P)$ in the polynomial order P, and is favorable to direct diagonalization for sufficiently large magnetic cells.

The polynomial order P scales like the spectral bandwidth of the dynamical matrix divided by the target energy resolution. If the specified resolution is too small (relative to the line broadening kernel), the calculated intensities will exhibit artificial oscillations in energy.

SpinWaveTheorySpiral(sys::System; k, axis, measure, regularization=1e-8)

Analogous to SpinWaveTheory, but interprets the provided system as having a generalized spiral order. This order is described by a single propagation wavevector k, which may be incommensurate. The axis vector defines the polarization plane via its surface normal. Typically the spin configuration in sys and the propagation wavevector k will be optimized using minimize_spiral_energy!. In contrast, axis will typically be determined from symmetry considerations.

The resulting object can be used to calculate the spin wave dispersion, or the structure factor via intensities_bands and intensities.

The algorithm for this calculation was developed in Toth and Lake, J. Phys.: Condens. Matter 27, 166002 (2015) and implemented in the SpinW code.

System(crystal::Crystal, moments, mode; dims=(1, 1, 1), seed=nothing)

A spin system is constructed from the Crystal unit cell, a specification of the spin moments symmetry-distinct sites, and a calculation mode. Interactions can be added to the system using, e.g., set_exchange!. The default supercell dimensions are 1×1×1 chemical cells, but this can be changed with dims.

Spin moments comprise a list of pairs, [i1 => Moment(...), i2 => ...], where i1, i2, ... are a complete set of symmetry-distinct atoms. Each Moment contains spin and $g$-factor information.

The two primary options for mode are :SUN and :dipole. In the former, each spin-$s$ degree of freedom is described as an SU(N) coherent state, i.e. a quantum superposition of $N = 2s + 1$ levels. This formalism can be useful to capture multipolar spin fluctuations or local entanglement effects.

Mode :dipole projects the SU(N) dynamics onto the restricted space of pure dipoles. In practice this means that Sunny will simulate Landau-Lifshitz dynamics, but single-ion anisotropy and biquadratic exchange interactions will be renormalized to improve accuracy. To disable this renormalization, use the mode :dipole_large_s which applies the $s → ∞$ classical limit. For details, see the documentation page: Interaction Renormalization.

An integer seed for the random number generator can optionally be specified to enable reproducible calculations.

All spins are initially polarized in the global $z$-direction.

Units(energy, length)

Physical constants in units of reference energy and length scales. Possible lengths are [:angstrom, :nm]. For atomic scale modeling, it is preferable to work in units of length=:angstrom, which follows the CIF file standard. Possible energy units are [:meV, :K, :THz, :inverse_cm, :T]. Kelvin is converted to energy via the Boltzmann constant $k_B$. Similarly, hertz is converted via the Planck constant $h$, inverse cm via the speed of light $c$, and tesla (field strength) via the Bohr magneton $μ_B$. For a given Units system, one can access any of the length and energy scale symbols listed above.

Examples

# Unit system with [energy] = meV and [length] = Å
units = Units(:meV, :angstrom)

# Use the Boltzmann constant kB to convert 1 kelvin into meV
@assert units.K ≈ 0.0861733326

# Use the Planck constant h to convert 1 THz into meV
@assert units.THz ≈ 4.135667696

# Use the constant h c to convert 1 cm⁻¹ into meV
@assert units.inverse_cm ≈ 0.1239841984

# Use the Bohr magneton μB to convert 1 tesla into meV
@assert units.T ≈ 0.05788381806

# The physical constant μ0 μB² in units of ų meV.
@assert u.vacuum_permeability ≈ 0.6745817653

add_sample!(sc::SampledCorrelations, sys::System)
add_sample!(sc::SampledCorrelationsStatic, sys::System)

Measure pair correlation data for the spin configuration in sys, and accumulate these statistics into sc. For a dynamical SampledCorrelations, this involves time-integration of the provided spin trajectory, recording correlations in both space and time. Conversely, SampledCorrelationsStatic, will record only spatial correlations for the single spin configuration that is provided.

Time-integration will update the spin configuration of sys in-place. To avoid this mutation, consider calling clone_system prior to add_sample!.

clone_correlations(sc::SampledCorrelations)

Create a copy of a SampledCorrelations.

clone_system(sys::System)

Creates a full clone of the system, such that mutable updates to one copy will not affect the other, and thread safety is guaranteed.

dispersion(swt::SpinWaveTheory, qpts)

Given a list of wavevectors qpts in reciprocal lattice units (RLU), returns excitation energies for each band. The return value ret is 2D array, and should be indexed as ret[band_index, q_index].

dmvec(D)

Antisymmetric matrix representation of the Dzyaloshinskii-Moriya pseudo-vector,

  [  0    D[3] -D[2]
   -D[3]   0    D[1]
    D[2] -D[1]   0  ]

By construction, Si'*dmvec(D)*Sj ≈ D⋅(Si×Sj) for any dipoles Si and Sj. This helper function is intended for use with set_exchange!.

domain_average(f, cryst, qpts; rotations, weights)

Calculate an average intensity for the reciprocal-space points qpts under a discrete set of rotations. Rotations, in global coordinates, may be given either as an axis-angle pair or as a 3×3 rotation matrix. Each rotation is weighted according to the elements in weights. The function f should accept a list of rotated q-points and return an intensities calculation.

Example

# 0, 120, and 240 degree rotations about the global z-axis
rotations = [([0,0,1], n*(2π/3)) for n in 0:2]
weights = [1, 1, 1]
res = domain_average(cryst, path; rotations, weights) do path_rotated
    intensities(swt, path_rotated; energies, kernel)
end
plot_intensities(res)

eachsite(sys::System)

An iterator over all Sites in the system.

enable_dipole_dipole!(sys::System, μ0_μB²)

Enables long-range interactions between magnetic dipole moments,

\[ -(μ_0/4π) ∑_{⟨ij⟩} [3 (μ_i⋅𝐫̂_{ij})(μ_j⋅𝐫̂_{ij}) - μ_i⋅μ_j] / r_{ij}^3\]

where the sum is over all pairs of spins (singly counted), including periodic images, regularized using the Ewald summation convention. The magnetic_moment is defined as $μ = -g μ_B 𝐒$, where $𝐒$ is the spin angular momentum dipole. The parameter μ0_μB² specifies the physical constant $μ_0 μ_B^2$, which has dimensions of length³-energy. Obtain this constant for a given system of Units via its vacuum_permeability property.

Example

units = Units(:meV, :angstrom)
enable_dipole_dipole!(sys, units.vacuum_permeability)

See also modify_exchange_with_truncated_dipole_dipole!.

energy(sys::System)

The total system energy. See also energy_per_site.

energy_per_site(sys::System)

The total system energy divided by the number of sites.

excitations!(T, tmp, swt::SpinWaveTheory, q)

Given a wavevector q, solves for the matrix T representing quasi-particle excitations, and returns a list of quasi-particle energies. Both T and tmp must be supplied as $2L×2L$ complex matrices, where $L$ is the number of bands for a single $𝐪$ value.

The columns of T are understood to be contracted with the Holstein-Primakoff bosons $[𝐛_𝐪, 𝐛_{-𝐪}^†]$. The first $L$ columns provide the eigenvectors of the quadratic Hamiltonian for the wavevector $𝐪$. The next $L$ columns of T describe eigenvectors for $-𝐪$. The return value is a vector with similar grouping: the first $L$ values are energies for $𝐪$, and the next $L$ values are the negation of energies for $-𝐪$.

excitations!(T, tmp, swt::SpinWaveTheorySpiral, q; branch)

Calculations on a SpinWaveTheorySpiral additionally require a branch index. The possible branches $(1, 2, 3)$ correspond to scattering processes $𝐪 - 𝐤, 𝐪, 𝐪 + 𝐤$ respectively, where $𝐤$ is the ordering wavevector. Each branch will contribute $L$ excitations, where $L$ is the number of spins in the magnetic cell. This yields a total of $3L$ excitations for a given momentum transfer $𝐪$.

excitations(swt::SpinWaveTheory, q)
excitations(swt::SpinWaveTheorySpiral, q; branch)

Returns a pair (energies, T) providing the excitation energies and eigenvectors. Prefer excitations! for performance, which avoids matrix allocations. See the documentation of excitations! for more details.

gaussian(; {fwhm, σ})

Returns the function exp(-x^2/2σ^2) / √(2π*σ^2). Either fwhm or σ must be specified, where fwhm = (2.355...) * σ is the full width at half maximum.

generate_mantid_script_from_binning_parameters(params::BinningParameters)

Generate a Mantid script which bins data according to the given BinningParameters.

global_position(sys::System, site::Site)

Position of a Site in global coordinates.

To precompute a full list of positions, one can use eachsite as below:

pos = [global_position(sys, site) for site in eachsite(sys)]

intensities(swt::SpinWaveTheory, qpts; energies, kernel, kT=0)
intensities(swt::SampledCorrelations, qpts; energies, kernel=nothing, kT)

Calculates pair correlation intensities for a set of $𝐪$-points in reciprocal space.

Traditional spin wave theory calculations are performed with an instance of SpinWaveTheory. One can alternatively use SpinWaveTheorySpiral to study generalized spiral orders with a single, incommensurate-$𝐤$ ordering wavevector. Another alternative is SpinWaveTheoryKPM, which may be faster than SpinWaveTheory for calculations on large magnetic cells (e.g., to study systems with disorder). In spin wave theory, a nonzero temperature kT will scale intensities by the quantum thermal occupation factor $|1 + n_B(ω)|$ where ``n_B(ω) = 1 / (exp(βω)

• 1)`` is the Bose function.

Intensities can also be calculated for SampledCorrelations associated with classical spin dynamics. In this case, thermal broadening will already be present, and the line-broadening kernel becomes an optional argument. Conversely, the parameter kT becomes required. If positive, it will introduce an intensity correction factor $|βω [1 + n_B(ω)]|$ that undoes the occupation factor for the classical Boltzmann distribution, and applies the quantum thermal occupation factor. The special choice kT = nothing will suppress the classical-to-quantum correction factor, and yield statistics consistent with the classical Boltzmann distribution.

intensities_bands(swt::SpinWaveTheory, qpts; kT=0)

Calculate spin wave excitation bands for a set of q-points in reciprocal space. This calculation is analogous to intensities, but does not perform line broadening of the bands.

intensities_static(sc::SpinWaveTheory, qpts; bounds=(-Inf, Inf), kT=0)
intensities_static(sc::SampledCorrelations, qpts; bounds=(-Inf, Inf), kT)
intensities_static(sc::SampledCorrelationsStatic, qpts)

Like intensities, but integrates the dynamical correlations $\mathcal{S}(𝐪, ω)$ over a range of energies $ω$. By default, the integration bounds are $(-∞, ∞)$, yielding the instantaneous (equal-time) correlations.

In SpinWaveTheory the integral will be realized as a sum over discrete bands. A SampledCorrelations object will have a finite grid of available energies, which will constrain the domain of integration. A SampledCorrelationsStatic object stores no dynamical data; here, the return value represents instantaneous correlations for the classical Boltzmann distribution.

The parameter kT can be used to account for the quantum thermal occupation of excitations at finite temperature. For details, see the documentation in intensities.

lattice_params(latvecs)

Compute the lattice parameters $(a, b, c, α, β, γ)$ for the three lattice vectors provided as columns of latvecs. The inverse mapping is lattice_vectors.

lattice_vectors(a, b, c, α, β, γ)

Return the lattice vectors, as columns of the $3×3$ output matrix, that define the shape of a crystallographic cell in global Cartesian coordinates. Conversely, one can view the output matrix as defining the global Cartesian coordinate system with respect to the lattice system.

The lattice constants $(a, b, c)$ have units of length, and the angles $(α, β, γ)$ are in degrees. The inverse mapping is lattice_params.

Example

latvecs = lattice_vectors(1, 1, 2, 90, 90, 120)
a1, a2, a3 = eachcol(latvecs)
@assert a1 ≈ [1, 0, 0]       # a1 always aligned with global x
@assert a2 ≈ [-1/2, √3/2, 0] # a2 always in global (x,y) plane
@assert a3 ≈ [0, 0, 2]       # a3 may generally be a combination of (x,y,z)

params, signal = load_nxs(filename; field="signal")

Given the name of a Mantid-exported MDHistoWorkspace file, load the BinningParameters and the signal from that file.

To load another field instead of the signal, specify e.g. field="errors_squared". Typical fields include errors_squared, mask, num_events, and signal.

lorentzian(; fwhm)

Returns the function (Γ/2) / (π*(x^2+(Γ/2)^2)) where fwhm = Γ is the full width at half maximum.

magnetic_moment(sys::System, site::Site)

Returns $- g 𝐒$, the local magnetic moment in units of the Bohr magneton. The spin dipole $𝐒$ and $g$-tensor may both be Site dependent.

merge_correlations(scs::Vector{SampledCorrelations)

Accumulate a list of SampledCorrelations into a single, summary SampledCorrelations. Useful for reducing the results of parallel computations.

minimize_energy!(sys::System{N}; maxiters=1000, method=Optim.ConjugateGradient(),
                 g_tol=1e-10, kwargs...) where N

Optimizes the spin configuration in sys to minimize energy. A total of maxiters iterations will be attempted. Convergence is reached when the root mean squared energy gradient goes below g_tol. The remaining kwargs will be forwarded to the optimize method of the Optim.jl package.

minimize_spiral_energy!(sys, axis; maxiters=10_000, k_guess=randn(sys.rng, 3))

Finds a generalized spiral order that minimizes the spiral_energy. This involves optimization of the spin configuration in sys, and the propagation wavevector $𝐤$, which will be returned in reciprocal lattice units (RLU). The axis vector normal to the polarization plane should be provided in global Cartesian coordinates, and will usually be determined by symmetry configurations. The initial k_guess will be random, unless otherwise provided.

See also suggest_magnetic_supercell to find a system shape that is approximately commensurate with the returned propagation wavevector $𝐤$.

modify_exchange_with_truncated_dipole_dipole!(sys::System, cutoff, μ0_μB²)

Like enable_dipole_dipole!, the purpose of this function is to introduce long-range dipole-dipole interactions between magnetic moments. Whereas enable_dipole_dipole! employs Ewald summation, this function instead employs real-space pair couplings with truncation at the specified cutoff distance. If the cutoff is relatively small, then this function may be faster than enable_dipole_dipole!.

polarize_spins!(sys::System, dir)

Polarize all spins in the system along the direction dir.

position_to_site(sys::System, r)

Converts a position r to four indices of a Site. The coordinates of r are given in units of the lattice vectors for the original crystal. This function can be useful for working with systems that have been reshaped using reshape_supercell.

Example

# Find the `site` at the center of a unit cell which is displaced by four
# multiples of the first lattice vector
site = position_to_site(sys, [4.5, 0.5, 0.5])

# Print the dipole at this site
println(sys.dipoles[site])

powder_average(f, cryst, radii, n; seed=0)

Calculate a powder-average over structure factor intensities. The radii, with units of inverse length, define spherical shells in reciprocal space. The Fibonacci lattice yields n points on the sphere, with quasi-uniformity. Sample points on different shells are decorrelated through random rotations. A consistent random number seed will yield reproducible results. The function f should accept a list of q-points and call intensities.

Example

radii = range(0.0, 3.0, 200)
res = powder_average(cryst, radii, 500) do qs
    intensities(swt, qs; energies, kernel)
end
plot_intensities(res)

primitive_cell_shape(cryst::Crystal)

Returns the shape of the primitive cell as a 3×3 matrix, in fractional coordinates of the conventional lattice vectors. May be useful for constructing inputs to reshape_supercell.

Examples

# Valid if `cryst` has not been reshaped
@assert cryst.prim_latvecs ≈ cryst.latvecs * primitive_cell_shape(cryst)

print_bond(cryst::Crystal, bond::Bond; b_ref::Bond)

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, ...

print_site(cryst, 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.

function print_stevens_expansion(op)

Prints a local Hermitian operator as a linear combination of Stevens operators. The operator op may be a finite-dimensional matrix or an abstract spin polynomial in the large-$s$ limit.

Examples

S = spin_matrices(2)
print_stevens_expansion(S[1]^4 + S[2]^4 + S[3]^4)
# Prints: (1/20)𝒪₄₀ + (1/4)𝒪₄₄ + 102/5

S = spin_matrices(Inf)
print_stevens_expansion(S[1]^4 + S[2]^4 + S[3]^4)
# Prints: (1/20)𝒪₄₀ + (1/4)𝒪₄₄ + (3/5)𝒮⁴

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 anisotropies.

print_symmetry_table(cryst::Crystal, max_dist)

Print symmetry information for all equivalence classes of sites and bonds, up to a maximum bond distance of max_dist. Equivalent to calling print_bond(cryst, b) for every bond b in reference_bonds(cryst, max_dist), where Bond(i, i, [0,0,0]) refers to a single site i.

print_wrapped_intensities(sys::System; nmax=10)

For Bravais lattices: Prints up to nmax wavevectors according to their instantaneous (static) structure factor intensities, listed in descending order. For non-Bravais lattices: Performs the same analysis for each spin sublattice independently; the output weights are naïvely averaged over sublattices, without incorporating phase shift information. This procedure therefore wraps all wavevectors into the first Brillouin zone. Each wavevector coordinate is given between $-1/2$ and $1/2$ in reciprocal lattice units (RLU). The output from this function will typically be used as input to suggest_magnetic_supercell.

Because this function does not incorporate phase information in its averaging over sublattices, the printed weights are not directly comparable with experiment. For that purpose, use SampledCorrelationsStatic instead.

propose_delta(magnitude)

Generate a proposal function that adds a Gaussian perturbation to the existing spin state. In :dipole mode, the procedure is to first introduce a random three-vector perturbation $𝐬′ = 𝐬 + |𝐬| ξ$ and then return the properly normalized spin $|𝐬| (𝐬′/|𝐬′|)$. Each component of the random vector $ξ$ is Gaussian distributed with a standard deviation of magnitude; the latter is dimensionless and typically smaller than one.

In :SUN mode, the procedure is analogous, but now involving Gaussian perturbations to each of the $N$ complex components of an SU(N) coherent state.

In the limit of very large magnitude, this function coincides with propose_uniform.

Consider also Langevin sampling, which is rejection free.

propose_flip

Function to propose pure spin flip updates in the context of a LocalSampler. Dipoles are flipped as $𝐬 → -𝐬$. SU(N) coherent states are flipped using the time-reversal operator.

propose_uniform

Function to propose a uniformly random spin update in the context of a LocalSampler. In :dipole mode, the result is a random three-vector with appropriate normalization. In :SUN mode, the result is a random SU(N) coherent state with appropriate normalization.

For low-temperature Monte Carlo simulations, uniform spin proposals can be very inefficient due to a high probability of rejection in the Metropolis accept/reject step. Consider also Langevin sampling, which is rejection free.

q_space_grid(cryst::Crystal, axis1, range1, axis2, range2; offset=[0,0,0], orthogonalize=false)
q_space_grid(cryst::Crystal, axis1, range1, axis2, range2, axis3, range3; orthogonalize=false)

Returns a 2D or 3D grid of q-points with uniform spacing. The volume shape is defined by (axis1, axis2, ...) in reciprocal lattice units (RLU). Elements of (range1, range2, ...) provide coefficients $c_i$ used to define grid positions,

    offset + c1 * axis1 + c2 * axis2 + ...

A nonzero offset is allowed only in the 2D case.

The first range parameter, range1, must be a regularly spaced list of coefficients, e.g., range1 = range(lo1, hi1, n). Subsequent range parameters may be a pair of bounds, without grid spacing information. For example, by selecting range2 = (lo2, hi2), an appropriate step-size will be inferred to provide an approximately uniform sampling density in global Cartesian coordinates.

The axes may be non-orthogonal. To extend to an orthohombic volume in global Cartesian coordinates, set orthogonalize=true.

For a 1D grid, use q_space_path instead.

q_space_path(cryst::Crystal, qs, n; labels=nothing)

Returns a 1D path consisting of n wavevectors sampled piecewise-linearly between the qs. Although the qs are provided in reciprocal lattice units (RLU), consecutive samples are spaced uniformly in the global (inverse-length) coordinate system. Optional labels can be associated with each special q-point, and will be used in plotting functions.

See also q_space_grid.

randomize_spins!(sys::System)

Randomizes all spins under appropriate the uniform distribution.

reference_bonds(cryst::Crystal, max_dist)

Returns a full list of bonds, one for each symmetry equivalence class, up to distance max_dist. The reference bond b for each equivalence class is selected according to a scoring system that prioritizes simplification of the elements in basis_for_symmetry_allowed_couplings(cryst, b).

remove_periodicity!(sys::System, flags)

Remove periodic interactions along each dimension d if flags[d] is true. The system must support inhomogeneous interactions via to_inhomogeneous.

Example

# Remove periodic boundaries along the 1st and 3rd dimensions
remove_periodicity!(sys::System, (true, false, true))

repeat_periodically(sys::System, counts::NTuple{3, Int})

Creates a System identical to sys but repeated a given number of times along each system axis according to the specified counts. This is a special case of the more general reshape_supercell.

See also repeat_periodically_as_spiral, which rotates the spins between periodic copies.

repeat_periodically_as_spiral(sys::System, counts::NTuple{3, Int}; k, axis)

Repeats the magnetic cell of System a number of times along each system axis according to the specified counts. Spins in each system image will be rotated according to the propagation wavevector k (in RLU) and the rotation axis (in global Cartesian coordinates). Coincides with repeat_periodically in the special case of k = [0, 0, 0]

See also minimize_spiral_energy! to find an energy-minimizing wavevector k and spin dipole configuration.

Example

k = minimize_spiral_energy!(sys, axis; k_guess=randn(3))
repeat_periodically_as_spiral(sys, counts; k, axis)

reshape_supercell(sys::System, shape)

Maps an existing System to a new one that has the shape and periodicity of a requested supercell. The columns of the $3×3$ integer matrix shape represent the supercell lattice vectors measured in units of the original crystal lattice vectors. Interactions, spins, and other settings will be inherited from sys.

In the special case that shape is a diagonal matrix, this function coincides with resize_supercell.

See also repeat_periodically.

resize_supercell(sys::System, dims::NTuple{3, Int})

Creates a System with a given number of conventional unit cells in each lattice vector direction. Interactions, spins, and other settings will be inherited from sys.

Equivalent to:

reshape_supercell(sys, [dims[1] 0 0; 0 dims[2] 0; 0 0 dims[3]])

See also reshape_supercell and repeat_periodically.

rotate_operator(A, R)

Rotates the local quantum operator A according to the $3×3$ rotation matrix R.

set_coherent!(sys::System, Z, site::Site)

Set a coherent spin state at a Site using the $N$ complex amplitudes in Z.

For a single quantum spin-$s$, these amplitudes will be interpreted in the eigenbasis of $Ŝ^z$. That is, Z[1] represents the amplitude for the basis state fully polarized along the $ẑ$-direction, and subsequent components represent states with decreasing angular momentum along this axis ($m = s, s-1, …, -s$).

set_dipole!(sys::System, dir, site::Site)

Polarize the spin at a Site along the direction dir.

set_dipoles_from_mcif!(sys::System, filename::AbstractString)

Load the magnetic supercell according to an mCIF file. System sys must already be resized to the correct supercell dimensions.

set_exchange!(sys::System, J, bond::Bond; biquad=0)

Sets an exchange interaction $𝐒_i⋅J 𝐒_j$ along the specified bond. This interaction will be propagated to equivalent bonds in consistency with crystal symmetry. Any previous interactions on these bonds will be overwritten. The parameter bond has the form Bond(i, j, offset), where i and j are atom indices within the unit cell, and offset is a displacement in unit cells.

As a convenience, scalar J can be used to specify a Heisenberg interaction. Also, the function dmvec(D) can be used to construct the antisymmetric part of the exchange, where D is the Dzyaloshinskii-Moriya pseudo-vector. The resulting interaction will be $𝐃⋅(𝐒_i×𝐒_j)$.

The optional numeric parameter biquad multiplies a scalar biquadratic interaction, $(𝐒_i⋅𝐒_j)^2$, with Interaction Renormalization if appropriate. For more general interactions, use set_pair_coupling! instead.

Examples

using LinearAlgebra

# Set a Heisenberg and DM interaction: 2Si⋅Sj + D⋅(Si×Sj)
D = [0, 0, 3]
set_exchange!(sys, 2I + dmvec(D), bond)

# The same interaction as an explicit exchange matrix
J = [2 3 0;
    -3 2 0;
     0 0 2]
set_exchange!(sys, J, bond)

set_exchange_at!(sys::System, J, site1::Site, site2::Site; biquad=0, offset=nothing)

Sets an exchange interaction `𝐒_i⋅J 𝐒_j along the single bond connecting two Sites, ignoring crystal symmetry. Any previous coupling on this bond will be overwritten. The system must support inhomogeneous interactions via to_inhomogeneous.

Use symmetry_equivalent_bonds to find (site1, site2, offset) values that are symmetry equivalent to a given Bond in the original system. For systems that are relatively small, the offset vector (in multiples of unit cells) will resolve ambiguities in the periodic wrapping.

See also set_exchange! for more details on specifying J and biquad. For more general couplings, use set_pair_coupling_at! instead.

set_field!(sys::System, B_μB)

Sets the external magnetic field $𝐁$ scaled by the Bohr magneton $μ_B$. This scaled field has units of energy and couples directly to the dimensionless magnetic_moment. At every site, the Zeeman coupling contributes an energy $+ (𝐁 μ_B) ⋅ (g 𝐒)$, involving the local $g$-tensor and spin angular momentum $𝐒$. Commonly, $g ≈ +2$ such that $𝐒$ is favored to anti-align with the applied field $𝐁$. Note that a given system of Units will implicitly use the Bohr magneton to convert between field and energy dimensions.

Example

# In units of meV, apply a 2 tesla field in the z-direction
units = Units(:meV, :angstrom)
set_field!(sys, [0, 0, 2] * units.T)

set_field_at!(sys::System, B_μB, site::Site)

Sets the external magnetic field $𝐁$ scaled by the Bohr magneton $μ_B$ for a single Site. This scaled field has units of energy and couples directly to the dimensionless magnetic_moment. Note that a given system of Units will implicitly use the Bohr magneton to convert between field and energy dimensions.

See the documentation of set_field! for more information.

set_onsite_coupling!(sys::System, op, i::Int)

Set the single-ion anisotropy for the ith atom of every unit cell, as well as all symmetry-equivalent atoms. The operator op may be provided as an abstract function of the local spin operators, as a polynomial of spin_matrices, or as a linear combination of stevens_matrices.

Examples

# An easy axis anisotropy in the z-direction
set_onsite_coupling!(sys, S -> -D*S[3]^3, i)

# The unique quartic single-ion anisotropy for a site with cubic point group
# symmetry
set_onsite_coupling!(sys, S -> 20*(S[1]^4 + S[2]^4 + S[3]^4), i)

# An equivalent expression of this quartic anisotropy, up to a constant shift
O = stevens_matrices(spin_label(sys, i))
set_onsite_coupling!(sys, O[4,0] + 5*O[4,4], i)

set_onsite_coupling_at!(sys::System, op, site::Site)

Sets the single-ion anisotropy operator op for a single Site, ignoring crystal symmetry. The system must support inhomogeneous interactions via to_inhomogeneous.

See also set_onsite_coupling!.

set_pair_coupling!(sys::System, op, bond)

Sets an arbitrary coupling op along bond. This coupling will be propagated to equivalent bonds in consistency with crystal symmetry. Any previous interactions on these bonds will be overwritten. The parameter bond has the form Bond(i, j, offset), where i and j are atom indices within the unit cell, and offset is a displacement in unit cells. The operator op may be provided as an anonymous function that accepts two spin dipole operators, or as a matrix that acts in the tensor product space of the two sites.

Examples

# Bilinear+biquadratic exchange involving 3×3 matrices J1 and J2
set_pair_coupling!(sys, (Si, Sj) -> Si'*J1*Sj + (Si'*J2*Sj)^2, bond)

# Equivalent expression using an appropriate fixed matrix representation
S = spin_matrices(1/2)
Si, Sj = to_product_space(S, S)
set_pair_coupling!(sys, Si'*J1*Sj + (Si'*J2*Sj)^2, bond)

See also spin_matrices, to_product_space.

set_pair_coupling_at!(sys::System, op, bond)

Sets an arbitrary coupling along the single bond connecting two Sites, ignoring crystal symmetry. Any previous coupling on this bond will be overwritten. The system must support inhomogeneous interactions via to_inhomogeneous.

Use symmetry_equivalent_bonds to find (site1, site2, offset) values that are symmetry equivalent to a given Bond in the original system. For systems that are relatively small, the offset vector (in multiples of unit cells) will resolve ambiguities in the periodic wrapping.

The operator op may be provided as an anonymous function that accepts two spin dipole operators, or as a matrix that acts in the tensor product space of the two sites. The documentation for set_pair_coupling! provides examples constructing op.

set_spin_rescaling!(sys, α)

In dipole mode, rescale all spin magnitudes $S → α S$. In SU(N) mode, rescale all SU(N) coherent states $Z → √α Z$ such that every expectation value rescales like $⟨A⟩ → α ⟨A⟩$.

set_vacancy_at!(sys::System, site::Site)

Make a single site nonmagnetic. Site includes a unit cell and a sublattice index.

spin_label(sys::System, i::Int)

If atom i carries a single spin-$s$ moment, then returns the half-integer label $s$. Otherwise, throws an error.

spin_matrices(s)

Returns a triple of $N×N$ spin matrices, where $N = 2s+1$. These are the generators of SU(2) in the spin-s representation.

If s == Inf, then the return values are abstract symbols denoting infinite-dimensional matrices that commute. These can be useful for repeating historical studies, or modeling micromagnetic systems. A technical discussion appears in the Sunny documentation page: Interaction Renormalization.

Example

S = spin_matrices(3/2)
@assert S'*S ≈ (3/2)*(3/2+1)*I
@assert S[1]*S[2] - S[2]*S[1] ≈ im*S[3]

S = spin_matrices(Inf)
@assert S[1]*S[2] - S[2]*S[1] == 0

See also print_stevens_expansion.

spiral_energy(sys::System; k, axis)

Returns the energy of a generalized spiral phase associated with the propagation wavevector k (in reciprocal lattice units, RLU) and an axis vector that is normal to the polarization plane (in global Cartesian coordinates).

When $𝐤$ is incommensurate, this calculation can be viewed as creating an infinite number of periodic copies of sys. The spins on each periodic copy are rotated about the axis vector, with the angle $θ = 2π 𝐤⋅𝐫$, where 𝐫 denotes the displacement vector between periodic copies of sys in multiples of the lattice vectors of the chemical cell.

The return value is the energy associated with one periodic copy of sys. The special case $𝐤 = 0$ yields result is identical to energy.

See also minimize_spiral_energy! and repeat_periodically_as_spiral.

spiral_energy_per_site(sys::System; k, axis)

The spiral_energy divided by the number of sites in sys. The special case $𝐤 = 0$ yields a result identical to energy_per_site.

ssf_custom(f, sys::System; apply_g=true, formfactors=nothing)

Specify measurement of the spin structure factor with a custom contraction function f. This function accepts a wavevector $𝐪$ in global Cartesian coordinates, and a 3×3 matrix with structure factor intensity components $\mathcal{S}^{αβ}(𝐪,ω)$. Indices $(α, β)$ denote dipole components in global coordinates. The return value of f can be any number or isbits type. With specific choices of f, one can obtain measurements such as defined in ssf_perp and ssf_trace.

By default, the g-factor or tensor is applied at each site, such that the structure factor components are correlations between the magnetic moment operators. Set apply_g = false to measure correlations between the bare spin operators.

The optional formfactors comprise a list of pairs [i1 => FormFactor(...), i2 => ...], where i1, i2, ... are a complete set of symmetry-distinct atoms, and each FormFactor implements $𝐪$-space attenuation for the given atom.

Intended for use with SpinWaveTheory and instances of SampledCorrelations.

Examples

# Measure all 3×3 structure factor components Sᵅᵝ
measure = ssf_custom((q, ssf) -> ssf, sys)

# Measure the structure factor trace Sᵅᵅ
measure = ssf_custom((q, ssf) -> real(sum(ssf)), sys)

See also the Sunny documentation on Structure Factor Conventions.

ssf_custom_bm(f, sys::System; u, v, apply_g=true, formfactors=nothing)

Specify measurement of the spin structure factor with a custom contraction function f. The interface is identical to ssf_custom except that f here receives momentum $𝐪$ and the 3×3 structure factor data $\mathcal{S}^{αβ}(𝐪, ω)$ in the basis of the Blume-Maleev axis system. The wavevectors u and v, provided in reciprocal lattice units, will be used to define the scattering plane. In global Cartesian coordinates, the three orthonormal BM axes (e1, e2, e3) are defined as follows:

e3 = normalize(u × v)  # normal to the scattering plane (u, v)
e1 = normalize(q)      # momentum transfer q within scattering plane
e2 = normalize(e3 × q) # perpendicular to q and in the scattering plane

Example

# Measure imaginary part of S²³ - S³² in the Blume-Maleev axis system for
# the scattering plane [0, K, L].
measure = ssf_custom_bm(sys; u=[0, 1, 0], v=[0, 0, 1]) do q, ssf
    imag(ssf[2,3] - ssf[3,2])
end

ssf_perp(sys::System; apply_g=true, formfactors=nothing)

Specify measurement of the spin structure factor with contraction by $(I-𝐪⊗𝐪/q^2)$. The contracted value provides an estimate of unpolarized scattering intensity. In the singular limit $𝐪 → 0$, the contraction matrix is replaced by its rotational average, $(2/3) I$.

This function is a special case of ssf_custom.

Example

# Select Co²⁺ form factor for atom 1 and its symmetry equivalents
formfactors = [1 => FormFactor("Co2")]
ssf_perp(sys; formfactors)

ssf_trace(sys::System; apply_g=true, formfactors=nothing)

Specify measurement of the spin structure factor, with trace over spin components. This quantity can be useful for checking quantum sum rules.

This function is a special case of ssf_custom.

standardize(cryst::Crystal; idealize=true)

Return the symmetry-inferred standardized crystal unit cell. If idealize=true, then the lattice vectors and site positions will be adapted. See "definitions and conventions" of the spglib documentation for more information.

step!(sys::System, dynamics)

Advance the spin configuration one dynamical time-step. The dynamics object may be a continuous spin dynamics, such as Langevin or ImplicitMidpoint, or it may be a discrete Monte Carlo sampling scheme such as LocalSampler.

stevens_matrices(s)

Returns a generator of Stevens operators in the spin-s representation. The return value O can be indexed as O[k,q], where $0 ≤ k ≤ 6$ labels an irrep of SO(3) and $-k ≤ q ≤ k$. This will produce an $N×N$ matrix where ``N = 2s

• 1``. Linear combinations of Stevens operators can be used as a "physical

basis" for decomposing local observables. To see this decomposition, use print_stevens_expansion.

If s == Inf, then symbolic operators will be returned. In this infinite dimensional limit, the Stevens operators become homogeneous polynomials of commuting spin operators.

Example

O = stevens_matrices(2)
S = spin_matrices(2)

A = (1/20)O[4,0] + (1/4)O[4,4] + (102/5)I
B = S[1]^4 + S[2]^4 + S[3]^4
@assert A ≈ B

See also spin_matrices and Interaction Renormalization.

subcrystal(cryst, types) :: Crystal

Filters sublattices of a Crystal by atom types, keeping the space group unchanged.

subcrystal(cryst, classes) :: Crystal

Filters sublattices of Crystal by equivalence classes, keeping the space group unchanged.

suggest_magnetic_supercell(ks; tol=1e-12, maxsize=100)

Suggests a magnetic supercell, in units of the crystal lattice vectors, that is consistent with periodicity of the wavevectors ks in RLU. If the wavevectors are incommensurate (with respect to the maximum supercell size maxsize), one can select a larger error tolerance tol to find a supercell that is almost commensurate.

Prints a $3×3$ matrix of integers that is suitable for use in reshape_supercell.

Examples

# A magnetic supercell for a single-Q structure. Will print
k1 = [0, -1/4, 1/4]
suggest_magnetic_supercell([k1])       # [1 0 0; 0 2 1; 0 -2 1]

# A larger magnetic supercell for a double-Q structure
k2 = [1/4, 0, 1/4]
suggest_magnetic_supercell([k1, k2])   # [1 2 2; -1 2 -2; -1 2 2]

# If given incommensurate wavevectors, find an approximate supercell that
# is exactly commensurate for nearby wavevectors.
suggest_magnetic_supercell([[0, 0, 1/√5], [0, 0, 1/√7]]; tol=1e-2)

# This prints [1 0 0; 0 1 0; 0 0 16], which becomes commensurate under the
# approximations `1/√5 ≈ 7/16` and `1/√7 ≈ 3/8`.

suggest_timestep(sys, integrator; tol)

Suggests a timestep for the numerical integration of spin dynamics according to a given error tolerance tol. The integrator should be Langevin or ImplicitMidpoint. The suggested $dt$ will be inversely proportional to the magnitude of the effective field $|dE/d𝐬|$ arising from the current spin configuration in sys. The recommended timestep $dt$ scales like √tol, which assumes second-order accuracy of the integrator.

The system sys should be initialized to an equilibrium spin configuration for the target temperature. Alternatively, a reasonably timestep estimate can be obtained from any low-energy spin configuration. For this, one can use randomize_spins! and then minimize_energy!.

Large damping magnitude or target temperature kT will tighten the timestep bound. If damping exceeds 1, it will rescale the suggested timestep by an approximate the factor $1/damping$. If kT is the largest energy scale, then the suggested timestep will scale like 1/(damping*kT). Quantification of numerical error for stochastic dynamics is subtle. The stochastic Heun integration scheme is weakly convergent of order-1, such that errors in the estimates of averaged observables may scale like dt. This implies that the tol argument may actually scale like the square of the true numerical error, and should be selected with this in mind.

symmetry_equivalent_bonds(sys::System, bond::Bond)

Given a Bond for the original (unreshaped) crystal, return all symmetry equivalent bonds in the System. Each returned bond is represented as a pair of Sites, which may be used as input to set_exchange_at! or set_pair_coupling_at!. Reverse bonds are not included in the iterator (no double counting).

Example

for (site1, site2, offset) in symmetry_equivalent_bonds(sys, bond)
    @assert site1 < site2
    set_exchange_at!(sys, J, site1, site2; offset)
end

to_inhomogeneous(sys::System)

Returns a copy of the system that allows for inhomogeneous interactions, which can be set using set_onsite_coupling_at!, set_exchange_at!, and set_vacancy_at!.

Inhomogeneous systems do not support symmetry-propagation of interactions or system reshaping.

to_product_space(A, B, ...)

Given lists of operators acting on local Hilbert spaces individually, return the corresponding operators that act on the tensor product space. In typical usage, the inputs will represent local physical observables and the outputs will be used to define quantum couplings.

@mix_proposals weight1 propose1 weight2 propose2 ...

Macro to generate a proposal function that randomly selects among the provided functions according to the provided probability weights. For use with LocalSampler.

Example

# A proposal function that proposes a spin flip 40% of the time, and a
# Gaussian perturbation 60% of the time.
@mix_proposals 0.4 propose_flip 0.6 propose_delta(0.2)

plot_spins(sys::System; arrowscale=1.0, color=:red, colorfn=nothing,
           colormap=:viridis, colorrange=nothing, show_cell=true, orthographic=false,
           ghost_radius=0, ndims=3, compass=true)

Plot the spin configuration defined by sys. Optional parameters are:

• arrowscale: Scale all arrows by dimensionless factor.
• color: Arrow colors. May be symbolic or numeric. If scalar, will be shared among all sites.
• colorfn: Function that dynamically maps from a site index to a numeric color value. Useful for animations.
• colormap, colorrange: Used to populate colors from numbers following Makie conventions.
• show_cell: Show original crystallographic unit cell.
• orthographic: Use orthographic camera perspective.
• ghost_radius: Show periodic images up to a given distance (length units).
• ndims: Spatial dimensions of system (1, 2, or 3).
• compass: If true, draw Cartesian axes in bottom left.

Calling notify on the return value will animate the figure.

plot_spins!(ax, sys::System; opts...)

Mutating variant of plot_spins that allows drawing into a single panel of a Makie figure.

Example

fig = Figure()
plot_spins!(fig[1, 1], sys1)
plot_spins!(fig[2, 1], sys2)
display(fig)

plot_intensities(res::BandIntensities; colormap=nothing, colorrange=nothing, allpositive=true,
                 saturation=0.9, sensitivity=0.0025, fwhm=nothing, ylims=nothing, units=nothing,
                 into=nothing, title="")
plot_intensities(res::Intensities; colormap=nothing, colorrange=nothing, allpositive=true, 
                 saturation=0.9, units=nothing, into=nothing, title="")
plot_intensities(res::StaticIntensities; colormap=nothing, colorrange=nothing, allpositive=true, 
                 saturation=0.9, units=nothing, into=nothing, title="")
plot_intensities(res::PowderIntensities; colormap=nothing, colorrange=nothing, allpositive=true, 
                 saturation=0.9, units=nothing, into=nothing, title="")

Keyword arguments:

• colormap: Color palette for plotting broadened intensities. See Makie docs for allowed values.
• colorrange: A lower and upper bound on intensities. For heatmaps, these bounds define the intensity values that saturate the colormap.
• allpositive: Should intensities be all positive, apart from numerical error? If true, the default colors will clip below zero intensity. If false, the default colors will be symmetric about zero intensity.
• saturation: If colorrange is not explicitly set, this dimensionless parameter defines the upper saturated intensity value as a quantile of maximum intensities taken over wavevectors.
• sensitivity: When plotting BandIntensities, this defines a lower bound on the visible intensities as a fraction of the upper saturated intensity value.
• fwhm: When plotting BandIntensities, this overrides the full-width at half-maximum value used for Gaussian broadening.
• ylims: Limits of the y-axis.
• units: A Units instance for labeling axes and performing conversions.
• into: A symbol for conversion into a new base energy unit (e.g. :meV, :K, etc.)
• title: An optional title for the plot.

plot_intensities!(panel, res::AbstractIntensities; opts...)

Mutating variant of plot_intensities that allows drawing into a single panel of a Makie figure.

Example

fig = Figure()
plot_intensities!(fig[1, 1], res1; title="Panel 1")
plot_intensities!(fig[2, 1], res2; title="Panel 2")
display(fig)

view_crystal(crystal::Crystal; refbonds=10, orthographic=false, ghost_radius=nothing, ndims=3, compass=true)
view_crystal(sys::System; ...)

Launches a graphical user interface to visualize the Crystal unit cell. If a System is provided, then the 3×3 exchange matrices for each bond will be depicted graphically.

• refbonds: By default, calculate up to 10 reference bonds using the reference_bonds function. An explicit list of reference bonds may also be provided.
• orthographic: Use orthographic camera perspective.
• ghost_radius: Show periodic images up to a given distance. Defaults to the cell size.
• ndims: Spatial dimensions of system (1, 2, or 3).
• compass: If true, draw Cartesian axes in bottom left.

export_vtk(filename,params::BinningParameters,data)

Export a VTK-compatible file to filename (do not include file extension when specifying the file name) which contains the data as VTK Cell Data on a grid parameterized by params.

At least one axis of the BinningParameters must be integrated over, since VTK does not support 4D data. See integrate_axes!.