Library API

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

Four indices identifying a single site in a System. The first three indices select the lattice cell and the last selects the sublattice (i.e., the 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. C.f. count_bins.

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.

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 an intensity_formula, 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)⟩ + \frac{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 $⟨j_l(s)⟩$ can be approximated as a sum of Gaussians,

\[⟨j_0(s)⟩ = A e^{-as^2} + B e^{-bs^2} + C e^{-cs^2} + D e^{-ds^2} + E \\ ⟨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. Built-in 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.

An alternative approach to sampling is Langevin, which may be preferred for simulating continuous spins, especially in the presence of long-range dipole-dipole interactions (cf. enable_dipole_dipole!).

SampledCorrelations

Basic data type for storing sampled correlation data. A SampleCorrelations is initialized by calling either dynamical_correlations or instant_correlations.

SpinInfo(atom::Int; S, g=2)

Characterizes the spin at a given atom index within the crystal unit cell. S is an integer multiple of 1/2 and gives the spin angular momentum in units of ħ. g is the g-factor or tensor, such that an angular momentum dipole $s$ produces a magnetic moment $+ g s$ in units of the Bohr magneton.

SpinWaveTheory(sys, energy_ϵ::Float64=1e-8)

Constructs an object to perform linear spin wave theory. Use it with dispersion and dssf functions.

The optional parameter energy_ϵ adds a small positive shift to the diagonal of the dynamical matrix $D$ to avoid numerical issues with zero-energy quasi-particle modes.

System(crystal::Crystal, latsize, infos, mode; seed::Int)

Construct a System of spins for a given Crystal symmetry. The latsize parameter determines the number of unit cells in each lattice vector direction. The infos parameter is a list of SpinInfo objects, one for each symmetry-inequivalent sublattice.

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 Strength Renormalization.

An optional seed may be provided to achieve reproducible random number generation.

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

Units(energy, length=:angstrom)

Physical constants in units of reference energy and length scales. Possible lengths are [:angstrom, :nm] and possible energies 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)

# 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_trajectory uses the spin configuration contained in the System to generate a correlation data and accumulate it into sc. For static structure factors, this involves analyzing the spin-spin correlations of the spin configuration provided. For a dynamic structure factor, a trajectory is calculated using the given spin configuration as an initial condition. The spin-spin correlations are then calculated in time and accumulated into sc.

This function will change the state of sys when calculating dynamical structure factor data. To preserve the initial state of sys, it must be saved separately prior to calling add_sample!. Alternatively, the initial spin configuration may be copied into a new System and this new System can be passed to add_sample!.

available_energies(sc::SampledCorrelations; negative_energies=false)

Return the ω values for the energy index of a SampledCorrelations. By default, only returns values for non-negative energies, which corresponds to the default output of intensities. Set negative_energies to true to retrieve all ω values.

available_wave_vectors(sc::SampledCorrelations; bzsize=(1,1,1))

Returns all wave vectors for which sc contains exact values. bsize specifies the number of Brillouin zones to be included.

axes_bincenters(params::BinningParameters)

Returns tick marks which label the bins of the histogram described by BinningParameters by their bin centers.

The following alternative syntax can be used to compute bin centers for a single axis:

axes_bincenters(binstart,binend,binwidth)

broaden_energy(sc::SampledCorrelations, vals, kernel::Function; negative_energies=false)

Performs a real-space convolution along the energy axis of an array of intensities. Assumes the format of the intensities array corresponds to what would be returned by intensities_interpolated. kernel must be a function that takes two numbers: kernel(ω, ω₀), where ω is a frequency, and ω₀ is the center frequency of the kernel. Sunny provides lorentzian for the most common use case:

newvals = broaden_energy(sc, vals, (ω, ω₀) -> lorentzian(fwhm=0.2)(ω-ω₀))

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.

count_bins(binstart,binend,binwidth)

Returns the number of bins in the binning scheme implied by binstart, binend, and binwidth. To count the bins in a BinningParameters, use params.numbins.

This function defines how partial bins are handled, so it should be used preferentially over computing the number of bins manually.

dispersion(swt::SpinWaveTheory, qs)

Computes the spin excitation energy dispersion relations given a SpinWaveTheory and an array of wave vectors qs. Each element $q$ of qs must be a 3-vector in units of reciprocal lattice units. I.e., $qᵢ$ is given in $2π/|aᵢ|$ with $|aᵢ|$ the lattice constant of the original chemical lattice.

The first indices of the returned array correspond to those of qs. A final index, corresponding to mode, is added to these. Each entry of the array is an energy.

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

dssf(swt::SpinWaveTheory, qs)

Given a SpinWaveTheory object, computes the dynamical spin structure factor,

\[ 𝒮^{αβ}(𝐤, ω) = 1/(2πN)∫dt ∑_𝐫 \exp[i(ωt - 𝐤⋅𝐫)] ⟨S^α(𝐫, t)S^β(0, 0)⟩,\]

using the result from linear spin-wave theory,

\[ 𝒮^{αβ}(𝐤, ω) = ∑_n |A_n^{αβ}(𝐤)|^2 δ[ω-ω_n(𝐤)].\]

qs is an array of wave vectors of arbitrary dimension. Each element $q$ of qs must be a 3-vector in reciprocal lattice units (RLU), i.e., in the basis of reciprocal lattice vectors.

The first indices of the returned array correspond to those of qs. A final index, corresponding to mode, is added to these. Each entry of this array is a tensor (3×3 matrix) corresponding to the indices $α$ and $β$.

dynamical_correlations(sys::System; dt, nω, ωmax, 
                       observables=nothing, correlations=nothing)

Creates an empty SampledCorrelations object for calculating and storing dynamical structure factor intensities $𝒮(𝐪,ω)$. Call add_sample! to accumulate data for the given configuration of a spin system. Internally, this will run a dynamical trajectory and measure time correlations. The $𝒮(𝐪,ω)$ data can be retrieved by calling intensities_interpolated. Alternatively, instant_intensities_interpolated will integrate out $ω$ to obtain $𝒮(𝐪)$, optionally applying classical-to-quantum correction factors.

Three keywords are required to specify the dynamics used for the trajectory calculation.

• dt: The time step used for calculating the trajectory from which dynamic spin-spin correlations are calculated. The trajectories are calculated with an ImplicitMidpoint integrator.
• ωmax: The maximum energy, $ω$, that will be resolved. Note that allowed values of ωmax are constrained by the given dt, so Sunny will choose the smallest possible value that is no smaller than the specified ωmax.
• nω: The number of energy bins to calculated between 0 and ωmax.

Additional keyword options are the following:

• observables: Allows the user to specify custom observables. The observables must be given as a list of complex N×N matrices or LinearMaps. It's recommended to name each observable, for example: observables = [:A => a_observable_matrix, :B => b_map, ...]. By default, Sunny uses the 3 components of the dipole, :Sx, :Sy and :Sz.
• correlations: Specify which correlation functions are calculated, i.e. which matrix elements $αβ$ of $𝒮^{αβ}(q,ω)$ are calculated and stored. Specified with a vector of tuples. By default Sunny records all auto- and cross-correlations generated by all observables. To retain only the xx and xy correlations, one would set correlations=[(:Sx,:Sx), (:Sx,:Sy)] or correlations=[(1,1),(1,2)].

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

# Valid for a system with lengths in Å and energies in meV
units = Units(:meV)
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.

gaussian(; {fwhm, σ})

Returns the function exp(-x^2/2σ^2) / √(2π*σ^2). Exactly one of fwhm or σ must be specified, where fwhm = (2.355...) * σ denotes 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)]

instant_correlations(sys::System; process_trajectory=:none, observables=nothing, correlations=nothing)

Creates an empty SampledCorrelations object for calculating and storing instantaneous structure factor intensities $𝒮(𝐪)$. Call add_sample! to accumulate data for the given configuration of a spin system. Call instant_intensities_interpolated to retrieve averaged $𝒮(𝐪)$ data.

Important note: When dealing with continuous (non-Ising) spins, consider creating using dynamical_correlations instead of instant_correlations. The former will provide full $𝒮(𝐪,ω)$ data, from which $𝒮(𝐪)$ can be obtained by integrating out $ω$. During this integration step, Sunny can incorporate temperature- and $ω$-dependent classical-to-quantum correction factors to produce more accurate $𝒮(𝐪)$ estimates. See instant_intensities_interpolated for more information.

The following optional keywords are available:

• observables: Allows the user to specify custom observables. The observables must be given as a list of complex N×N matrices or LinearMaps. It's recommended to name each observable, for example: observables = [:A => a_observable_matrix, :B => b_map, ...]. By default, Sunny uses the 3 components of the dipole, :Sx, :Sy and :Sz.
• correlations: Specify which correlation functions are calculated, i.e. which matrix elements $αβ$ of $𝒮^{αβ}(q,ω)$ are calculated and stored. Specified with a vector of tuples. By default Sunny records all auto- and cross-correlations generated by all observables. To retain only the xx and xy correlations, one would set correlations=[(:Sx,:Sx), (:Sx,:Sy)] or correlations=[(1,1),(1,2)].

instant_intensities_interpolated(sc::SampledCorrelations, qs, formula::ClassicalIntensityFormula; kwargs...)

Return $𝒮(𝐪)$ intensities at wave vectors qs. The functionality is very similar to intensities_interpolated, except the returned array has dimensions identical to qs. If called on a SampledCorrelations with dynamical information, i.e., $𝒮(𝐪,ω)$, the $ω$ information is integrated out.

integrate_axes!(params::BinningParameters; axes)

Integrate over one or more axes of the histogram by setting the number of bins in that axis to 1. Examples:

integrate_axes!(params; axes = [2,3])
integrate_axes!(params; axes = 2)

integrated_gaussian(; {fwhm, σ})

Returns the function erf(x/√2σ)/2, which is the integral of gaussian over the range $[0, x]$. Exactly one of fwhm or σ must be specified, where fwhm = (2.355...) * σ denotes the full width at half maximum. Intended for use with intensities_binned.

integrated_lorentzian(; fwhm)

Returns the function atan(2x/Γ)/π, which is the integral of lorentzian over the range $[0, x]$, where Γ = fwhm is the full width at half maximum. Intended for use with intensities_binned.

dispersion, intensities = intensities_bands(swt::SpinWaveTheory, qs, formula::SpinWaveIntensityFormula)

Computes the scattering intensities at each energy band for each momentum transfer q in qs, according to linear spin wave theory and the given intensity formula. The formula must have a delta-function kernel, e.g.:

formula = intensity_formula(swt, :perp, formula; kernel=delta_function_kernel)

or else the bands will be broadened, and their intensity can not be computed.

The outputs will be arrays with indices identical to qs, with the last index giving the band index. dispersions reports the energy of each band, while intensities reports the scattering intensity.

intensity, counts = intensities_binned(sc::SampledCorrelations, params::BinningParameters, formula; integrated_kernel)

Given correlation data contained in a SampledCorrelations and BinningParameters describing the shape of a histogram, compute the intensity and normalization for each histogram bin using a given intensity_formula.

The BinningParameters are expected to accept (q,ω) in R.L.U. for the (possibly reshaped) crystal associated with sc.

This is an alternative to intensities_interpolated which bins the scattering intensities into a histogram instead of interpolating between them at specified qs values. See unit_resolution_binning_parameters for a reasonable default choice of BinningParameters which roughly emulates intensities_interpolated with interpolation = :round.

If a function integrated_kernel(Δω) is passed, it will be used as the CDF of a kernel function for energy broadening. For example, integrated_kernel = Δω -> atan(Δω/η)/pi (c.f. integrated_lorentzian implements Lorentzian broadening with parameter η. Energy-dependent energy broadening can be achieved by providing an integrated_kernel(ω,Δω) whose first argument is the energy transfer ω.

Currently, energy broadening is only supported if the BinningParameters are such that the first three axes are purely spatial and the last (energy) axis is [0,0,0,1].

intensities_broadened(swt::SpinWaveTheory, qs, ωs, formula)

Computes the scattering intensities at each (q,ω) according to Linear Spin Wave Theory and the given intensity formula. The required formula must have a non-delta-function kernel, e.g.:

formula = intensity_formula(swt, :perp; kernel = lorentzian(0.05))

or else the intensity at ωs which are not exactly on the dispersion curve can not be calculated.

The intensity is computed at each wave vector in qs and each energy in ωs. The output will be an array with indices identical to qs, with the last index matching ωs.

Note that qs is an array of wave vectors of arbitrary dimension. Each element $q$ of qs must be a 3-wavevector in absolute units.

intensities_interpolated(sc::SampledCorrelations, qs, formula:ClassicalIntensityFormula; interpolation=nothing, negative_energies=false)

The basic function for retrieving $𝒮(𝐪,ω)$ information from a SampledCorrelations. Maps an array of wave vectors qs to an array of structure factor intensities, including an additional energy index. The values of $ω$ associated with the energy index can be retrieved by calling available_energies. The three coordinates of each wave vector are measured in reciprocal lattice units, i.e., multiples of the reciprocal lattice vectors.

• interpolation: Since $𝒮(𝐪, ω)$ is calculated on a finite lattice, data is only available at discrete wave vectors. By default, Sunny will round a requested q to the nearest available wave vector. Linear interpolation can be applied by setting interpolation=:linear.
• negative_energies: If set to true, Sunny will return the periodic extension of the energy axis. Most users will not want this.

A custom intensity formula can be specifed by providing a function intensity = f(q,ω,correlations) and specifying which correlations it requires:

intensity_formula(f,sc::SampledCorrelations, required_correlations; kwargs...)

The function is intended to be specified using do notation. For example, this custom formula sums the off-diagonal correlations:

required = [(:Sx,:Sy),(:Sy,:Sz),(:Sx,:Sz)]
intensity_formula(sc,required,return_type = ComplexF64) do k, ω, off_diagonal_correlations
    sum(off_diagonal_correlations)
end

If your custom formula returns a type other than Float64, use the return_type keyword argument to flag this.

formula = intensity_formula(sc::SampledCorrelations)

Establish a formula for computing the intensity of the discrete scattering modes (q,ω) using the correlation data $𝒮^{αβ}(q,ω)$ stored in the SampledCorrelations. The formula returned from intensity_formula can be passed to intensities_interpolated or intensities_binned.

intensity_formula(sc,...; kT = Inf, formfactors = ...)

There are keyword arguments providing temperature and form factor corrections:

• kT: If a temperature is provided, the intensities will be rescaled by a temperature- and ω-dependent classical-to-quantum factor. kT should be specified when making comparisons with spin wave calculations or experimental data. If kT is not specified, infinite temperature (no correction) is assumed.
• formfactors: To apply form factor corrections, provide this keyword with a list of FormFactors, one for each symmetry-distinct site in the crystal. The order of FormFactors must correspond to the order of site symmetry classes, e.g., as they appear when printed in display(crystal).

formula = intensity_formula(swt::SpinWaveTheory; kernel = ...)

Establish a formula for computing the scattering intensity by diagonalizing the hamiltonian $H(q)$ using Linear Spin Wave Theory.

If kernel = delta_function_kernel, then the resulting formula can be used with intensities_bands.

If kernel is an energy broadening kernel function, then the resulting formula can be used with intensities_broadened. Energy broadening kernel functions can either be a function of Δω only, e.g.:

kernel = Δω -> ...

or a function of both the energy transfer ω and of Δω, e.g.:

kernel = (ω,Δω) -> ...

The integral of a properly normalized kernel function over all Δω is one.

intensity_formula([swt or sc], contraction_mode::Symbol)

Sunny has several built-in formulas that can be selected by setting contraction_mode to one of these values:

• :trace (default), which yields $\operatorname{tr} 𝒮(q,ω) = ∑_α 𝒮^{αα}(q,ω)$
• :perp, which contracts $𝒮^{αβ}(q,ω)$ with the dipole factor $δ_{αβ} - q_{α}q_{β}$, returning the unpolarized intensity.
• :full, which will return all elements $𝒮^{αβ}(𝐪,ω)$ without contraction.

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 correspond to the conventional unit cell defined by the lattice constants $(a, b, c)$ and the angles $(α, β, γ)$ in degrees. The inverse mapping is lattice_params.

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.

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_binned(sc::SampledCorrelations, radial_binning_parameters; formula
                     ω_binning_parameters, integrated_kernel = nothing, bzsize = nothing)

This function emulates the experimental situation of "powder averaging," where only the magnitude (and not the direction) of the momentum transfer is resolvable. The intensities are binned similarly to intensities_binned, but the histogram x-axis is |k| in absolute units, which is a nonlinear function of kx,ky,kz. The y-axis is energy.

Radial binning parameters are specified as tuples (start,end,bin_width), e.g. radial_binning_parameters = (0,6π,6π/55).

Energy broadening is supported in the same way as intensities_binned, and this function accepts the same kind of intensity_formula.

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; digits=4)

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

For use with LocalSampler.

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.

randomize_spins!(sys::System)

Randomizes all spins under appropriate the uniform distribution.

reciprocal_space_path(cryst::Crystal, qs, density)

Returns a pair (path, xticks). The path return value is a list of wavevectors that samples linearly between the provided wavevectors qs. The xticks return value can be used to label the special $𝐪$ values on the x-axis of a plot.

Special note about units: the wavevectors qs must be provided in reciprocal lattice units (RLU) for the given crystal, but the sampling density must be specified in the global frame. Specifically, the density is given as number of sample points per unit of radian inverse length, where the unit of length is the same as that used to specify the lattice vectors of the Crystal. The path will therefore include more samples between q-points that are further apart in absolute Fourier distance.

reciprocal_space_path_bins(sc,qs,density,args...;kwargs...)

Takes a list of wave vectors, qs in R.L.U., and builds a series of histogram BinningParameters whose first axis traces a path through the provided points. The second and third axes are integrated over according to the args and kwargs, which are passed through to slice_2D_binning_parameters.

Also returned is a list of marker indices corresponding to the input points, and a list of ranges giving the indices of each histogram x-axis within a concatenated histogram. The density parameter is given in samples per reciprocal lattice unit (R.L.U.).

reciprocal_space_shell(cryst::Crystal, radius, n)

Sample n points on the reciprocal space sphere with a given radius (units of inverse length).

Examples

# Sample wavevectors on the sphere at fixed density
reciprocal_space_shell(cryst, r, 4π*r^2*density)

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, dims)

Remove periodic interactions along the dimensions where dims 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{N}, counts::NTuple{3,Int}) where N

Creates a System identical to sys but repeated a given number of times in each dimension, specified by the tuple counts.

See also reshape_supercell.

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.

resize_supercell(sys::System{N}, latsize::NTuple{3,Int}) where N

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

Convenience function for:

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

See also reshape_supercell.

rotate_operator(A, R)

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

rotation_in_rlu(cryst::Crystal, axis, angle)

Returns a $3×3$ matrix that rotates wavevectors in reciprocal lattice units (RLU). The axis vector is a real-space direction in absolute units (but arbitrary magnitude), and the angle is in radians.

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

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

For a standard SpinInfo, these amplitudes will be interpreted in the eigenbasis of $𝒮̂ᶻ$. 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 appropriate Interaction Strength Renormalization. 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)
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_spiral_order!(sys; k, axis, S0)

Initializes the system with a spiral order described by the wavevector k, an axis of rotation axis, and an initial dipole direction S0 at the real-space origin. The wavevector is expected in repicrocal lattice units (RLU), while the direction vectors axis and S0 are expected in global Cartesian coordinates.

Example

# Spiral order for a wavevector propagating in the direction of the first
# reciprocal lattice vector (i.e., orthogonal to the lattice vectors ``𝐚_2``
# and ``𝐚_3``), repeating with a period of 3 lattice constants, and spiraling
# about the ``ẑ``-axis. The spin at the origin will point in the direction
# ``𝐒_0 = ŷ + ẑ``.  Here, ``(x̂, ŷ, ẑ)`` are the axes of Cartesian coordinate
# system in the global frame.
set_spiral_order!(sys; k=[1/3, 0, 0], axis=[0, 0, 1], S0=[0, 1, 1])

See also set_spiral_order_on_sublattice!.

set_spiral_order_on_sublattice!(sys, i; k, axis, S0)

Initializes sublattice i with a spiral order described by the wavevector k, an axis of rotation axis, and an initial dipole direction S0. The phase is selected such that the spin at sys.dipole[1,1,1,i] will point in the direction of S0. The wavevector is expected in repicrocal lattice units (RLU), while the direction vectors axis and S0 are expected in global Cartesian coordinates.

This function is not available for systems with reshaped unit cells.

See also set_spiral_order!.

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

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

slice_2D_binning_parameter(sc::SampledCorrelations, cut_from_q, cut_to_q, cut_bins::Int64, cut_width::Float64; plane_normal = [0,0,1],cut_height = cutwidth)

Creates BinningParameters which make a cut along one dimension of Q-space.

The x-axis of the resulting histogram consists of cut_bins-many bins ranging from cut_from_q to cut_to_q. The width of the bins in the transverse direciton is controlled by cut_width and cut_height.

The binning in the transverse directions is defined in the following way, which sets their normalization and orthogonality properties:

cut_covector = normalize(cut_to_q - cut_from_q)
transverse_covector = normalize(plane_normal × cut_covector)
cotransverse_covector = normalize(transverse_covector × cut_covector)

In other words, the axes are orthonormal with respect to the Euclidean metric.

If the cut is too narrow, there will be very few scattering vectors per bin, or the number per bin will vary substantially along the cut. If the output appears under-resolved, try increasing cut_width.

The four axes of the resulting histogram are:

1. Along the cut
2. Fist transverse Q direction
3. Second transverse Q direction
4. Energy

This function can be used without reference to a SampledCorrelations using this alternate syntax to manually specify the bin centers for the energy axis:

slice_2D_binning_parameter(ω_bincenters, cut_from, cut_to,...)

where ω_bincenters specifies the energy axis, and both cut_from and cut_to are arbitrary covectors, in any units.

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

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

unit_resolution_binning_parameters(sc::SampledCorrelations)

Create BinningParameters which place one histogram bin centered at each possible (q,ω) scattering vector of the crystal. This is the finest possible binning without creating bins with zero scattering vectors in them.

This function can be used without reference to a SampledCorrelations using an alternate syntax to manually specify the bin centers for the energy axis and the lattice size:

unit_resolution_binning_parameters(ω_bincenters,latsize,[reciprocal lattice vectors])

The last argument may be a 3x3 matrix specifying the reciprocal lattice vectors, or a Crystal.

Lastly, binning parameters for a single axis may be specifed by their bin centers:

(binstart,binend,binwidth) = unit_resolution_binning_parameters(bincenters::Vector{Float64})

@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, dims=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).
• dims: 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.

view_crystal(crystal::Crystal; refbonds=10, orthographic=false, ghost_radius=nothing, dims=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.
• dims: 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!.