sunny 0.2.0

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

from savannah.ornl.gov/gravitas/jupyter-notebook:25.2

#run apt-get update -y && \

#    apt-get upgrade -y

user jovyan

ENV JULIA_CPU_TARGET="generic;tigerlake,clone_all;znver3,clone_all"

run julia -e 'using Pkg; Pkg.Registry.update(); Pkg.add(["Sunny", "WGLMakie"]); Pkg.precompile()'

ADD jupyter_lab_config.py /home/$NB_USER/.jupyter/

#ENV JULIA_CPU_TARGET="generic;tigerlake,clone_all;znver3,clone_all"

ENV JULIA_CPU_TARGET="generic;sandybridge,-xsaveopt,clone_all;haswell,-rdrnd,base(1);x86-64-v4,-rdrnd,base(1)"

run julia -e 'using Pkg; Pkg.Registry.update(); \

    Pkg.add([Pkg.PackageSpec(;name="JLLWrappers", version="1.7.0"), \

             Pkg.PackageSpec(;name="Sunny", version="0.7.5"), \

             Pkg.PackageSpec(;name="WGLMakie"), \

             Pkg.PackageSpec(;url="https://github.com/quantumsteve/Bonito.jl",rev="jupyter_base_url")]); \

    Pkg.precompile()'

run chmod -R 777 /opt/julia

#ADD jupyter_lab_config.py /home/$NB_USER/.jupyter/

#ADD jupyter_notebook_config.py /home/$NB_USER/.jupyter/

user root

COPY examples ./examples

run chmod -R 777 examples

%% Cell type:markdown id: tags:

# 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).

%% Cell type:markdown id: tags:

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

%% Cell type:code id: tags:

``` julia

using Sunny, WGLMakie

@assert pkgversion(Sunny) >= v"0.7.4"

@assert pkgversion(Sunny) >= v"0.7.5"

```

%% Cell type:markdown id: tags:

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

%% Cell type:markdown id: tags:

### Units system

%% Cell type:markdown id: tags:

The `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.

%% Cell type:code id: tags:

``` julia

units = Units(:meV, :angstrom);

```

%% Cell type:markdown id: tags:

### 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` 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.

%% Cell type:code id: tags:

``` julia

a = 8.5031 # (Å)

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

```

%% Cell type:markdown id: tags:

Construct a `Crystal` cell from spacegroup 227 in the ITA standard

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

lattice.

%% Cell type:code id: tags:

``` julia

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

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

```

%% Cell type:markdown id: tags:

`view_crystal` launches an interface for interactive inspection and

symmetry analysis.

%% Cell type:code id: tags:

``` julia

view_crystal(cryst)

```

%% Cell type:markdown id: tags:

### Spin system

%% Cell type:markdown id: tags:

A `System` will define the spin model. Each cobalt atom carries

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

`Moment` 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.

%% Cell type:code id: tags:

``` julia

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

```

%% Cell type:markdown id: tags:

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!` 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` with `sys` now shows

the antiferromagnetic Heisenberg interactions as blue polkadot spheres.

%% Cell type:code id: tags:

``` julia

J = +0.63 # (meV)

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

view_crystal(sys)

```

%% Cell type:markdown id: tags:

### Optimizing spins

%% Cell type:markdown id: tags:

To search for the ground state, call `randomize_spins!` and

`minimize_energy!` in sequence. For this problem, optimization

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

`plot_spins`, where spins are colored according to their global

$z$-component.

%% Cell type:code id: tags:

``` julia

randomize_spins!(sys)

minimize_energy!(sys)

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

```

%% Cell type:markdown id: tags:

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

%% Cell type:code id: tags:

``` julia

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

```

%% Cell type:markdown id: tags:

### Reshaping the magnetic cell

%% Cell type:markdown id: tags:

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

Columns of the `primitive_cell` matrix provide the primitive lattice

vectors as multiples of the conventional cubic lattice vectors.

%% Cell type:code id: tags:

``` julia

shape = primitive_cell(cryst)

```

%% Cell type:markdown id: tags:

Reduce the magnetic cell size using `reshape_supercell`. Verify that

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

%% Cell type:code id: tags:

``` julia

sys_prim = reshape_supercell(sys, shape)

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

```

%% Cell type:markdown id: tags:

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

the larger conventional cubic cell for context.

%% Cell type:code id: tags:

``` julia

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

```

%% Cell type:markdown id: tags:

### Spin wave theory

%% Cell type:markdown id: tags:

With this primitive cell, we will perform a `SpinWaveTheory`

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

`ssf_perp` 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` for Co²⁺ dampens intensities at large $𝐪$.

%% Cell type:code id: tags:

``` julia

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

measure = ssf_perp(sys_prim; formfactors)

swt = SpinWaveTheory(sys_prim; measure)

```

%% Cell type:markdown id: tags:

Select `lorentzian` broadening with a full-width at half-maximum

(FWHM) of 0.8 meV.

%% Cell type:code id: tags:

``` julia

kernel = lorentzian(fwhm=0.8)

```

%% Cell type:markdown id: tags:

Define a `q_space_path` 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.

%% Cell type:code id: tags:

``` julia

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

path = q_space_path(cryst, qs, 500)

```

%% Cell type:markdown id: tags:

Calculate single-crystal scattering `intensities` along this path, for

energies between 0 and 6 meV. Use `plot_intensities` to visualize the

result.

%% Cell type:code id: tags:

``` julia

energies = range(0, 6, 300)

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

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

```

%% Cell type:markdown id: tags:

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

average over all possible crystal orientations. Use `powder_average`

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.

%% Cell type:code id: tags:

``` julia

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")

```

%% Cell type:markdown id: tags:

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)

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

%% Cell type:markdown id: tags:

### What's next?

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

  ported to Sunny.

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

  next CoRh₂O₄ tutorial 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.

%% Cell type:markdown id: tags:

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

In the previous tutorial, 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.

%% Cell type:markdown id: tags:

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

the expected Néel order.

%% Cell type:code id: tags:

``` julia

using Sunny, WGLMakie

@assert pkgversion(Sunny) >= v"0.7.4"

@assert pkgversion(Sunny) >= v"0.7.5"

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])

```

%% Cell type:markdown id: tags:

Use `repeat_periodically` 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.

%% Cell type:code id: tags:

``` julia

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

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

```

%% Cell type:markdown id: tags:

### Langevin dynamics for sampling

%% Cell type:markdown id: tags:

We will be using a `Langevin` 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.

%% Cell type:code id: tags:

``` julia

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

```

%% Cell type:markdown id: tags:

Use `suggest_timestep` 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.

%% Cell type:code id: tags:

``` julia

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

langevin.dt = 0.025;

```

%% Cell type:markdown id: tags:

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

energy per site at each time step.

%% Cell type:code id: tags:

``` julia

energies = [energy_per_site(sys)]

for _ in 1:1000

    step!(sys, langevin)

    push!(energies, energy_per_site(sys))

end

```

%% Cell type:markdown id: tags:

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

smaller than necessary; increasing it will make the remaining simulations

faster.

%% Cell type:code id: tags:

``` julia

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

langevin.dt = 0.042;

```

%% Cell type:markdown id: tags:

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.

%% Cell type:code id: tags:

``` julia

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

```

%% Cell type:markdown id: tags:

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.

%% Cell type:code id: tags:

``` julia

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

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

```

%% Cell type:markdown id: tags:

### Static structure factor

%% Cell type:markdown id: tags:

Use `SampledCorrelationsStatic` to estimate spatial correlations for

configurations in classical thermal equilibrium. Measure `ssf_perp`,

which is appropriate for unpolarized neutron scattering. Include the

`FormFactor` for Co2⁺. Each call to `add_sample!` will

accumulate data for the current spin snapshot.

%% Cell type:code id: tags:

``` julia

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`

```

%% Cell type:markdown id: tags:

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

measurements to approximately decorrelate the sample in thermal equilibrium.

%% Cell type:code id: tags:

``` julia

for _ in 1:20

    for _ in 1:100

        step!(sys, langevin)

    end

    add_sample!(sc, sys)

end

```

%% Cell type:markdown id: tags:

Use `q_space_grid` 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.

%% Cell type:code id: tags:

``` julia

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

```

%% Cell type:markdown id: tags:

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.

%% Cell type:code id: tags:

``` julia

res = intensities_static(sc, grid)

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

```

%% Cell type:markdown id: tags:

### Dynamical structure factor

%% Cell type:markdown id: tags:

To collect statistics for the _dynamical_ structure factor intensities

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

`SampledCorrelations`. 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.

%% Cell type:code id: tags:

``` julia

dt = 2*langevin.dt

energies = range(0, 6, 50)

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

```

%% Cell type:markdown id: tags:

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

equilibrium. Now, however, each call to `add_sample!` 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.

%% Cell type:code id: tags:

``` julia

for _ in 1:5

    for _ in 1:100

        step!(sys, langevin)

    end

    add_sample!(sc, sys)

end

```

%% Cell type:markdown id: tags:

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

and a sampling density.

%% Cell type:code id: tags:

``` julia

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)

```

%% Cell type:markdown id: tags:

Calculate and plot the intensities along this path.

%% Cell type:code id: tags:

``` julia

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

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

```

%% Cell type:markdown id: tags:

### Powder averaged intensity

%% Cell type:markdown id: tags:

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

%% Cell type:code id: tags:

``` julia

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")

```

%% Cell type:markdown id: tags:

# 3. Multi-flavor spin wave simulations of FeI₂

This tutorial illustrates various advanced features such as symmetry analysis,

energy minimization, and spin wave theory with multi-flavor bosons.

Our context will be FeI₂, an effective spin-1 material with strong single-ion

anisotropy. Quadrupolar fluctuations give rise to a single-ion bound state

that is observable in neutron scattering, and cannot be described by a

dipole-only model. We will use the linear spin wave theory of SU(3) coherent

states (i.e. 2-flavor bosons) to model the magnetic spectrum of FeI₂. The

original study was performed in [Bai et al., Nature Physics **17**, 467–472

(2021)](https://doi.org/10.1038/s41567-020-01110-1).

<img src="https://raw.githubusercontent.com/SunnySuite/Sunny.jl/main/docs/src/assets/FeI2_crystal.jpg" style="float: left;" width="400">

The Fe atoms are arranged in stacked triangular layers. The effective spin

Hamiltonian takes the form,

$$

\mathcal{H}=\sum_{(i,j)} 𝐒_i ⋅ J_{ij} 𝐒_j - D\sum_i \left(S_i^z\right)^2,

$$

where the exchange matrices $J_{ij}$ between bonded sites $(i,j)$ include

competing ferromagnetic and antiferromagnetic interactions. This model also

includes a strong easy axis anisotropy, $D > 0$.

%% Cell type:markdown id: tags:

Load packages.

%% Cell type:code id: tags:

``` julia

using Sunny, WGLMakie

@assert pkgversion(Sunny) >= v"0.7.4"

@assert pkgversion(Sunny) >= v"0.7.5"

```

%% Cell type:markdown id: tags:

Construct the chemical cell of FeI₂ by specifying the lattice vectors and the

full set of atoms.

%% Cell type:code id: tags:

``` julia

units = Units(:meV, :angstrom)

a = b = 4.05012  # Lattice constants for triangular lattice (Å)

c = 6.75214      # Spacing between layers (Å)

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

positions = [[0, 0, 0], [1/3, 2/3, 1/4], [2/3, 1/3, 3/4]]

types = ["Fe", "I", "I"]

cryst = Crystal(latvecs, positions; types)

```

%% Cell type:markdown id: tags:

Observe that the space group 'P -3 m 1' (164) has been inferred, as well as

point group symmetries. Only the Fe atoms are magnetic, so we focus on them

with `subcrystal`. Importantly, the new crystal retains the symmetry

information for spacegroup 164.

%% Cell type:code id: tags:

``` julia

cryst = subcrystal(cryst, "Fe")

view_crystal(cryst)

```

%% Cell type:markdown id: tags:

### Symmetry analysis

The command `print_symmetry_table` provides a list of all the

symmetry-allowed interactions out to 8 Å.

%% Cell type:code id: tags:

``` julia

print_symmetry_table(cryst, 8.0)

```

%% Cell type:markdown id: tags:

The allowed $g$-tensor is expressed as a 3×3 matrix in the free coefficients

`A`, `B`, ... The allowed single-ion anisotropy is expressed as a linear

combination of Stevens operators. The latter correspond to polynomials of the

spin operators, as we will describe below.

The allowed exchange interactions are given as 3×3 matrices for representative

bonds. The notation `Bond(i, j, n)` indicates a bond between atom indices `i`

and `j`, with cell offset `n`. Note that the order of the pair $(i, j)$ is

significant if the exchange tensor contains antisymmetric

Dzyaloshinskii–Moriya (DM) interactions.

The bonds can be visualized in the `view_crystal` interface. By default,

`Bond(1, 1, [1,0,0])` is toggled on, to show the 6 nearest-neighbor Fe-Fe

bonds on a triangular lattice layer. Toggling `Bond(1, 1, [0,0,1])` shows the

Fe-Fe bond between layers, etc.

%% Cell type:markdown id: tags:

### Defining the spin model

%% Cell type:markdown id: tags:

Construct a `System` with spin $s=1$ and $g=2$ for the Fe ions.

Recall that quantum spin-1 is, in reality, a linear combination of basis

states $|m⟩$ with well-defined angular momentum $m = -1, 0, 1$. FeI₂ has a

strong easy-axis anisotropy, which stabilizes a single-ion bound state of zero

angular momentum, $|m=0⟩$. Such a bound state is inaccessible to traditional

spin wave theory, which works with dipole expectation values of fixed

magnitude. This physics is, however, well captured with a theory of SU(_N_)

coherent states, where $N = 2S+1 = 3$ is the number of levels. Activate this

generalized theory by selecting `:SUN` mode instead of `:dipole` mode.

An optional `seed` for random number generation can be used to to make the

calculation exactly reproducible.

%% Cell type:code id: tags:

``` julia

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

```

%% Cell type:markdown id: tags:

Set the exchange interactions for FeI₂ following the fits of [Bai et

al.](https://doi.org/10.1038/s41567-020-01110-1)

%% Cell type:code id: tags:

``` julia

J1pm   = -0.236 # (meV)

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