Add precompilation/example scripts

FROM code.ornl.gov:4567/gravitas/containers/jupyter-notebook:24.07

FROM savannah.ornl.gov/gravitas/jupyter-notebook:24.12

COPY examples ./examples

RUN chmod -R 777 examples

# RUN conda update -y -n base -c conda-forge conda

RUN conda install -y -c conda-forge phonopy

USER root

RUN apt update -y && \

    apt-get install -y libopenblas0-serial && \

    apt install -y quantum-espresso

RUN conda install -y -c conda-forge phonopy qe

# # 1. Spin wave simulations of CoRh₂O₄

#

# This tutorial introduces Sunny through its features for performing

# conventional spin wave theory calculations. We consider the crystal CoRh₂O₄

# and reproduce the calculations of [Ge et al., Phys. Rev. B **96**,

# 064413 (2017)](https://doi.org/10.1103/PhysRevB.96.064413).

# ### Get Julia and Sunny

# 

# Sunny is implemented in Julia, which allows for interactive development (like

# Python or Matlab) while also providing high numerical efficiency (like C++ or

# Fortran). New Julia users should begin with our **[Getting

# Started](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia)**

# guide. Sunny requires Julia 1.10 or later.

#

# From the Julia prompt, load Sunny and also [GLMakie](https://docs.makie.org/)

# for graphics.

using Sunny, WGLMakie

# If these packages are not yet installed, Julia will offer to install them. If

# executing this tutorial gives an error, you may need to update Sunny and

# GLMakie from the [built-in package

# manager](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia#the-built-in-julia-package-manager).

# ### Units system

# The [`Units`](@ref Units) object selects reference energy and length scales,

# and uses these to provide physical constants. For example, `units.K` returns

# one kelvin as 0.086 meV, where the Boltzmann constant is implicit.

units = Units(:meV, :angstrom);

# ### Crystal cell

#

# A crystallographic cell may be loaded from a `.cif` file, or specified from

# atom positions and types.

#

# Start by defining the shape of the conventional chemical cell. CoRh₂O₄ has

# cubic spacegroup 227 (Fd-3m). Its lattice constants are 8.5 Å, and the cell

# angles are 90°. With this information, [`lattice_vectors`](@ref) constructs a

# 3×3 matrix `latvecs`. Columns of `latvecs` define the lattice vectors ``(𝐚_1,

# 𝐚_2, 𝐚_3)`` in the global Cartesian coordinate system. Conversely, columns

# of `inv(latvecs)` define the global Cartesian axes ``(\hat{x}, \hat{y},

# \hat{z})`` in components of the lattice vectors.

a = 8.5031 # (Å)

latvecs = lattice_vectors(a, a, a, 90, 90, 90)

# Construct a [`Crystal`](@ref) cell from spacegroup 227 in the ITA standard

# setting. Cobalt atoms belong to Wyckoff 8a, which is the diamond cubic

# lattice.

positions = [[1/8, 1/8, 1/8]]

cryst = Crystal(latvecs, positions, 227; types=["Co"])

# [`view_crystal`](@ref) launches an interface for interactive inspection and

# symmetry analysis.

view_crystal(cryst)

# ### Spin system

# A [`System`](@ref) will define the spin model. Each cobalt atom carries

# quantum spin ``s = 3/2``, with a ``g``-factor of 2. Specify this

# [`Moment`](@ref) data for cobalt atom 1. By symmetry, the same moment data

# also applies to cobalt atoms 2, 3, ... 7. The option `:dipole` indicates a

# traditional model type, for which quantum spin is modeled as a dipole

# expectation value.

sys = System(cryst, [1 => Moment(s=3/2, g=2)], :dipole)

# Ge et al. demonstrated that inelastic neutron scattering data for CoRh₂O₄ is

# well modeled by antiferromagnetic nearest neighbor exchange, `J = 0.63` meV.

# Call [`set_exchange!`](@ref) with the bond that connects atom 1 to atom 3, and

# has zero displacement between chemical cells. Consistent with the symmetries

# of spacegroup 227, this interaction will be propagated to all other

# nearest-neighbor bonds. Calling [`view_crystal`](@ref) with `sys` now shows

# the antiferromagnetic Heisenberg interactions as blue polkadot spheres.

J = +0.63 # (meV)

set_exchange!(sys, J, Bond(2, 3, [0, 0, 0]))

view_crystal(sys)

# ### Optimizing spins

# To search for the ground state, call [`randomize_spins!`](@ref) and

# [`minimize_energy!`](@ref) in sequence. For this problem, optimization

# converges rapidly to the expected Néel order. See this with

# [`plot_spins`](@ref), where spins are colored according to their global

# ``z``-component.

randomize_spins!(sys)

minimize_energy!(sys)

plot_spins(sys; color=[S[3] for S in sys.dipoles])

# The diamond lattice is bipartite, allowing each spin to perfectly anti-align

# with its 4 nearest-neighbors. Each of these 4 bonds contribute ``-J s^2`` to

# the total energy. Two sites participate in each bond, so the energy per site

# is ``-2 J s^2``. Check this by calling [`energy_per_site`](@ref).

@assert energy_per_site(sys) ≈ -2J*(3/2)^2

# ### Reshaping the magnetic cell

# The most compact magnetic cell for this Néel order is the primitive unit cell.

# Columns of the [`primitive_cell`](@ref) matrix provide the primitive lattice

# vectors as multiples of the conventional cubic lattice vectors.

shape = primitive_cell(cryst)

# Reduce the magnetic cell size using [`reshape_supercell`](@ref). Verify that

# the energy per site is unchanged after the reshaping the supercell.

sys_prim = reshape_supercell(sys, shape)

@assert energy_per_site(sys_prim) ≈ -2J*(3/2)^2

# Plotting `sys_prim` shows the two spins within the primitive cell, as well as

# the larger conventional cubic cell for context.

plot_spins(sys_prim; color=[S[3] for S in sys_prim.dipoles])

# ### Spin wave theory

# With this primitive cell, we will perform a [`SpinWaveTheory`](@ref)

# calculation of the structure factor ``\mathcal{S}(𝐪,ω)``. The measurement

# [`ssf_perp`](@ref) indicates projection of the spin structure factor

# ``\mathcal{S}(𝐪,ω)`` perpendicular to the direction of momentum transfer, as

# appropriate for unpolarized neutron scattering. The isotropic

# [`FormFactor`](@ref) for Co²⁺ dampens intensities at large ``𝐪``.

formfactors = [1 => FormFactor("Co2")]

measure = ssf_perp(sys_prim; formfactors)

swt = SpinWaveTheory(sys_prim; measure)

# Select [`lorentzian`](@ref) broadening with a full-width at half-maximum

# (FWHM) of 0.8 meV. 

kernel = lorentzian(fwhm=0.8)

# Define a [`q_space_path`](@ref) that connects high-symmetry points in

# reciprocal space. The ``𝐪``-points are given in reciprocal lattice units

# (RLU) for the _original_ cubic cell. For example, `[1/2, 1/2, 0]` denotes the

# sum of the first two reciprocal lattice vectors, ``𝐛_1/2 + 𝐛_2/2``. A total

# of 500 ``𝐪``-points will be sampled along the path.

qs = [[0, 0, 0], [1/2, 0, 0], [1/2, 1/2, 0], [0, 0, 0]]

path = q_space_path(cryst, qs, 500)

# Calculate single-crystal scattering [`intensities`](@ref) along this path, for

# energies between 0 and 6 meV. Use [`plot_intensities`](@ref) to visualize the

# result.

energies = range(0, 6, 300)

res = intensities(swt, path; energies, kernel)

plot_intensities(res; units, title="CoRh₂O₄ LSWT")

# Sometimes experimental data is only available as a powder average, i.e., as an

# average over all possible crystal orientations. Use [`powder_average`](@ref)

# to simulate these intensities. Each ``𝐪``-magnitude defines a spherical shell

# in reciprocal space. Consider 200 radii from 0 to 3 inverse angstroms, and

# collect 2000 random samples per spherical shell. As configured, this

# calculation completes in about two seconds. Had we used the conventional cubic

# cell, the calculation would be an order of magnitude slower.

radii = range(0, 3, 200) # (1/Å)

res = powder_average(cryst, radii, 2000) do qs

    intensities(swt, qs; energies, kernel)

end

plot_intensities(res; units, saturation=1.0, title="CoRh₂O₄ Powder Average")

# This result can be compared to experimental neutron scattering data

# from Fig. 5 of [Ge et al.](https://doi.org/10.1103/PhysRevB.96.064413)

# ```@raw html

# <img width="95%" src="https://raw.githubusercontent.com/SunnySuite/Sunny.jl/main/docs/src/assets/CoRh2O4_intensity.jpg">

# ```

# ### What's next?

#

# * For more spin wave calculations of this type, browse the [SpinW tutorials

#   ported to Sunny](@ref "SW01 - FM Heisenberg chain").

# * Spin wave theory neglects thermal fluctuations of the magnetic order. The

#   [next CoRh₂O₄ tutorial](@ref "2. Landau-Lifshitz dynamics of CoRh₂O₄ at

#   finite *T*") demonstrates how to sample spins in thermal equilibrium, and

#   measure dynamical correlations from the classical spin dynamics.

# * Sunny also offers features that go beyond the dipole approximation of a

#   quantum spin via the theory of SU(_N_) coherent states. This can be

#   especially useful for systems with strong single-ion anisotropy, as

#   demonstrated in the [FeI₂ tutorial](@ref "3. Multi-flavor spin wave

#   simulations of FeI₂").

# # 2. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T*

#

# In the [previous tutorial](@ref "1. Spin wave simulations of CoRh₂O₄"), we

# used spin wave theory to calculate the dynamical spin structure factor of

# CoRh₂O₄. Here, we perform a similar calculation using equilibrium samples from

# the Boltzmann distribution at finite ``T``. For each sampled spin

# configuration, we will simulate the classical Landau-Lifshitz spin dynamics

# and extract dynamical spin-spin correlations. After applying a

# classical-to-quantum correction factor, the resulting intensities can be

# compared to inelastic neutron scattering data.

# Construct the CoRh₂O₄ antiferromagnet as before. Energy minimization yields

# the expected Néel order.

using Sunny, WGLMakie

units = Units(:meV, :angstrom)

a = 8.5031 # (Å)

latvecs = lattice_vectors(a, a, a, 90, 90, 90)

cryst = Crystal(latvecs, [[1/8, 1/8, 1/8]], 227)

sys = System(cryst, [1 => Moment(s=3/2, g=2)], :dipole)

J = 0.63 # (meV)

set_exchange!(sys, J, Bond(2, 3, [0, 0, 0]))

randomize_spins!(sys)

minimize_energy!(sys)

plot_spins(sys; color=[S[3] for S in sys.dipoles])

# Use [`repeat_periodically`](@ref) to extend the system to 10×10×10 chemical

# unit cells. The ground state Néel order is retained. Increasing the system

# size further would reduce finite-size artifacts and increase momentum-space

# resolution, but would also make the simulations slower.

sys = repeat_periodically(sys, (10, 10, 10))

plot_spins(sys; color=[S[3] for S in sys.dipoles])

# ### Langevin dynamics for sampling

# We will be using a [`Langevin`](@ref) spin dynamics to thermalize the system.

# This dynamics is a variant of the Landau-Lifshitz equation that incorporates

# noise and dissipation terms, which are linked by a fluctuation-dissipation

# theorem. The temperature 16 K ≈ 1.38 meV is slightly above ordering for this

# model. The dimensionless `damping` magnitude sets a timescale for coupling to

# the implicit thermal bath; 0.2 is usually a good choice.

langevin = Langevin(; damping=0.2, kT=16*units.K)

# Use [`suggest_timestep`](@ref) to select an integration timestep. A

# dimensionless error tolerance of `1e-2` is usually a good choice. The

# suggested timestep will vary according to the magnetic configuration. It is

# reasonable to start from an energy-minimized configuration.

suggest_timestep(sys, langevin; tol=1e-2)

langevin.dt = 0.025;

# Now run a Langevin trajectory to sample spin configurations. Keep track of the

# energy per site at each time step.

energies = [energy_per_site(sys)]

for _ in 1:1000

    step!(sys, langevin)

    push!(energies, energy_per_site(sys))

end

# From the relaxed spin configuration, we can learn that `dt` was a little

# smaller than necessary; increasing it will make the remaining simulations

# faster.

suggest_timestep(sys, langevin; tol=1e-2)

langevin.dt = 0.042;

# Plot energy versus time using the Makie

# [`lines`](https://docs.makie.org/stable/reference/plots/lines) function. The

# plateau suggests that the system has reached thermal equilibrium.

lines(energies, color=:blue, figure=(size=(600,300),), axis=(xlabel="Timesteps", ylabel="Energy (meV)"))

# Plot the spins colored by their alignment with a reference spin at the origin.

# The field `sys.dipoles` is a 4D array storing the spin dipole data. The first

# three indices of label the chemical cell, while the fourth index labels an

# atom within the cell. Note that Julia arrays use 1-based indexing. Thermal

# fluctuations are apparent in the plot.

S0 = sys.dipoles[1,1,1,1]

plot_spins(sys; color=[S'*S0 for S in sys.dipoles])

# ### Static structure factor

# Use [`SampledCorrelationsStatic`](@ref) to estimate spatial correlations for

# configurations in classical thermal equilibrium. Measure [`ssf_perp`](@ref),

# which is appropriate for unpolarized neutron scattering. Include the

# [`FormFactor`](@ref) for Co2⁺. Each call to [`add_sample!`](@ref) will

# accumulate data for the current spin snapshot.

formfactors = [1 => FormFactor("Co2")]

measure = ssf_perp(sys; formfactors)

sc = SampledCorrelationsStatic(sys; measure)

add_sample!(sc, sys)    # Accumulate the newly sampled structure factor into `sf`

# Collect 20 additional samples. Perform 100 Langevin time-steps between

# measurements to approximately decorrelate the sample in thermal equilibrium.

for _ in 1:20

    for _ in 1:100

        step!(sys, langevin)

    end

    add_sample!(sc, sys)

end

# Use [`q_space_grid`](@ref) to define a slice of momentum space ``[H, K, 0]``,

# where ``H`` and ``K`` each range from -10 to 10 in RLU. This command produces

# a 200×200 grid of sample points.

grid = q_space_grid(cryst, [1, 0, 0], range(-10, 10, 200), [0, 1, 0], (-10, 10))

# Calculate and plot the instantaneous structure factor on the slice. Selecting

# `saturation = 1.0` sets the color saturation point to the maximum intensity

# value. This is reasonable because we are above the ordering temperature, and

# do not have sharp Bragg peaks.

res = intensities_static(sc, grid)

plot_intensities(res; saturation=1.0, title="Static Intensities at T = 16 K")

# ### Dynamical structure factor

# To collect statistics for the _dynamical_ structure factor intensities

# ``\mathcal{S}(𝐪,ω)`` at finite temperature, use

# [`SampledCorrelations`](@ref). It requires a range of `energies` to resolve,

# which will be associated with frequencies of the classical spin dynamics. The

# integration timestep `dt` can be somewhat larger than that used by the

# Langevin dynamics.

dt = 2*langevin.dt

energies = range(0, 6, 50)

sc = SampledCorrelations(sys; dt, energies, measure)

# Like before, use Langevin dynamics to sample spin configurations from thermal

# equilibrium. Now, however, each call to [`add_sample!`](@ref) will run a

# classical spin dynamics trajectory and measure dynamical correlations. To make

# the tutorial run quickly, we average over just 5 trajectories. To make a

# publication quality figure, this number should be significantly increased for

# better statistics.

for _ in 1:5

    for _ in 1:100

        step!(sys, langevin)

    end

    add_sample!(sc, sys)

end

# Select points that define a piecewise-linear path through reciprocal space,

# and a sampling density.

qs = [[3/4, 3/4,   0],

      [  0,   0,   0],

      [  0, 1/2, 1/2],

      [1/2,   1,   0],

      [  0,   1,   0],

      [1/4,   1, 1/4],

      [  0,   1,   0],

      [  0,  -4,   0]]

path = q_space_path(cryst, qs, 500)

# Calculate and plot the intensities along this path.

res = intensities(sc, path; energies, langevin.kT)

plot_intensities(res; units, title="Intensities at 16 K")

# ### Powder averaged intensity

# Define spherical shells in reciprocal space via their radii, in absolute units

# of 1/Å. For each shell, calculate and average the intensities at 350

# ``𝐪``-points

radii = range(0, 3.5, 200) # (1/Å)

res = powder_average(cryst, radii, 350) do qs

    intensities(sc, qs; energies, langevin.kT)

end

plot_intensities(res; units, title="Powder Average at 16 K")

# # 4. Generalized spin dynamics of FeI₂ at finite *T*

# The [previous FeI₂ tutorial](@ref "3. Multi-flavor spin wave simulations of

# FeI₂") used multi-flavor spin wave theory to calculate the dynamical spin

# structure factor. This tutorial performs an analogous calculation at finite

# temperature using the [classical dynamics of SU(_N_) coherent

# states](https://doi.org/10.1103/PhysRevB.106.054423).

#

# Compared to spin wave theory, classical spin dynamics in real-space is

# typically much slower, and is limited in ``𝐪``-space resolution. The

# approach, however, allows for thermal fluctuations, can be used to explore

# [finite temperature phases](https://doi.org/10.1103/PhysRevB.109.014427), and

# enables the study of [highly non-equilibrium

# processes](https://doi.org/10.1103/PhysRevB.106.235154).

#

# The structure of this tutorial largely follows the [previous study of CoRh₂O₄

# at finite *T*](@ref "2. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T*").

# The main difference is that CoRh₂O₄ can be well described with `:dipole` mode,

# whereas FeI₂ is best modeled using `:SUN` mode, owing to its strong easy-axis

# anisotropy.

#

# Construct the FeI₂ system as [previously described](@ref "3. Multi-flavor spin

# wave simulations of FeI₂").

using Sunny, WGLMakie

units = Units(:meV, :angstrom)

a = b = 4.05012

c = 6.75214

latvecs = lattice_vectors(a, b, c, 90, 90, 120)

cryst = Crystal(latvecs, [[0,0,0]], 164; types=["Fe"])

sys = System(cryst, [1 => Moment(s=1, g=2)], :SUN)

J1pm   = -0.236

J1pmpm = -0.161

J1zpm  = -0.261

J2pm   = 0.026

J3pm   = 0.166

J′0pm  = 0.037

J′1pm  = 0.013

J′2apm = 0.068

J1zz   = -0.236

J2zz   = 0.113

J3zz   = 0.211

J′0zz  = -0.036

J′1zz  = 0.051

J′2azz = 0.073

J1xx = J1pm + J1pmpm

J1yy = J1pm - J1pmpm

J1yz = J1zpm

set_exchange!(sys, [J1xx 0.0 0.0; 0.0 J1yy J1yz; 0.0 J1yz J1zz], Bond(1,1,[1,0,0]))

set_exchange!(sys, [J2pm 0.0 0.0; 0.0 J2pm 0.0; 0.0 0.0 J2zz], Bond(1,1,[1,2,0]))

set_exchange!(sys, [J3pm 0.0 0.0; 0.0 J3pm 0.0; 0.0 0.0 J3zz], Bond(1,1,[2,0,0]))

set_exchange!(sys, [J′0pm 0.0 0.0; 0.0 J′0pm 0.0; 0.0 0.0 J′0zz], Bond(1,1,[0,0,1]))

set_exchange!(sys, [J′1pm 0.0 0.0; 0.0 J′1pm 0.0; 0.0 0.0 J′1zz], Bond(1,1,[1,0,1]))

set_exchange!(sys, [J′2apm 0.0 0.0; 0.0 J′2apm 0.0; 0.0 0.0 J′2azz], Bond(1,1,[1,2,1]))

D = 2.165

set_onsite_coupling!(sys, S -> -D*S[3]^2, 1)

# ### Thermalization

# To study thermal fluctuations in real-space, use a large system size with

# 16×16×4 copies of the chemical cell.

sys = resize_supercell(sys, (16, 16, 4))

# Direct optimization via [`minimize_energy!`](@ref) is susceptible to trapping

# in a local minimum. An alternative approach is to simulate the system using

# [`Langevin`](@ref) spin dynamics. This requires a bit more set-up, but allows

# sampling from thermal equilibrium at any target temperature. Select the

# temperature 2.3 K ≈ 0.2 meV. This temperature is small enough to magnetically

# order, but large enough so that the dynamics can readily overcome local energy

# barriers and annihilate defects.

langevin = Langevin(; damping=0.2, kT=2.3*units.K)

# Use [`suggest_timestep`](@ref) to select an integration timestep for the error

# tolerance `tol=1e-2`. Initializing `sys` to some low-energy configuration

# usually works well.

randomize_spins!(sys)

minimize_energy!(sys; maxiters=10)

suggest_timestep(sys, langevin; tol=1e-2)

langevin.dt = 0.03;

# Run a Langevin trajectory for 10,000 time-steps and plot the spins. The

# magnetic order is present, but may be difficult to see.

for _ in 1:10_000

    step!(sys, langevin)

end

plot_spins(sys; color=[S[3] for S in sys.dipoles])

# Verify the expected two-up, two-down spiral magnetic order by calling

# [`print_wrapped_intensities`](@ref). A single propagation wavevector ``±𝐤``

# dominates the static intensity in ``\mathcal{S}(𝐪)``, indicating the expected

# 2 up, 2 down magnetic spiral order. A smaller amount of intensity is spread

# among many other wavevectors due to thermal fluctuations.

print_wrapped_intensities(sys)

# Thermalization has not substantially altered the suggested `dt`.

suggest_timestep(sys, langevin; tol=1e-2)

# ### Structure factor in the paramagnetic phase

# The remainder of this tutorial will focus on the paramagnetic phase.

# Re-thermalize the system to the temperature of 3.5 K ≈ 0.30 meV.

langevin.kT = 3.5 * units.K

for _ in 1:10_000

    step!(sys, langevin)

end

# The suggested timestep has increased slightly. Following this suggestion will

# make the simulations a bit faster.

suggest_timestep(sys, langevin; tol=1e-2)

langevin.dt = 0.040;

# Collect dynamical spin structure factor data using

# [`SampledCorrelations`](@ref). This procedure involves sampling spin

# configurations from thermal equilibrium and using the [spin dynamics of

# SU(_N_) coherent states](https://arxiv.org/abs/2204.07563) to estimate

# dynamical correlations. With proper classical-to-quantum corrections, the

# associated structure factor intensities ``S^{αβ}(q,ω)`` can be compared with

# finite-temperature inelastic neutron scattering data. Incorporate the

# [`FormFactor`](@ref) appropriate to Fe²⁺. 

dt = 2*langevin.dt

energies = range(0, 7.5, 120)

formfactors = [1 => FormFactor("Fe2"; g_lande=3/2)]

sc = SampledCorrelations(sys; dt, energies, measure=ssf_perp(sys; formfactors))

# The function [`add_sample!`](@ref) will collect data by running a dynamical

# trajectory starting from the current system configuration. 

add_sample!(sc, sys)

# To collect additional data, it is required to re-sample the spin configuration

# from the thermal distribution. Statistical error is reduced by fully

# decorrelating the spin configurations between calls to `add_sample!`.

for _ in 1:2

    for _ in 1:1000               # Enough steps to decorrelate spins

        step!(sys, langevin)

    end

    add_sample!(sc, sys)

end

# Perform an intensity calculation for two special ``𝐪``-points in reciprocal

# lattice units (RLU). A classical-to-quantum rescaling of normal mode

# occupations will be performed according to the temperature `kT`. The large

# statistical noise could be reduced by averaging over more thermal samples.

res = intensities(sc, [[0, 0, 0], [0.5, 0.5, 0.5]]; energies, langevin.kT)

fig = lines(res.energies, res.data[:, 1]; axis=(xlabel="meV", ylabel="Intensity"), label="(0,0,0)")

lines!(res.energies, res.data[:, 2]; label="(π,π,π)")

axislegend()

fig

# Next, we will measure intensities along a [`q_space_path`](@ref) that connects

# high symmetry points. Because this is a real-space calculation, data is only

# available for discrete ``𝐪`` modes, with resolution that scales inversely to

# linear system size. Intensities at ``ω = 0`` dominate, so to enhance

# visibility, we restrict the color range empirically.

qs = [[0,   0, 0],  # List of wave vectors that define a path

      [1,   0, 0],

      [0,   1, 0],

      [1/2, 0, 0],

      [0,   1, 0],

      [0,   0, 0]] 

qpath = q_space_path(cryst, qs, 500)

res = intensities(sc, qpath; energies, langevin.kT)

plot_intensities(res; colorrange=(0.0, 1.0), title="Intensities at T = 2.3 K")

# One can also view the intensity along a [`q_space_grid`](@ref) for a fixed

# energy value. Alternatively, use [`intensities_static`](@ref) to integrate

# over all available energies.

grid = q_space_grid(cryst, [1, 0, 0], range(-1.5, 1.5, 300), [0, 1, 0], (-1.5, 1.5); orthogonalize=true)

res = intensities(sc, grid; energies=[3.5], langevin.kT)

plot_intensities(res; title="Intensity slice at ω = 3.5 meV")