Skip to content

Commit

Permalink
Docs fine tuning (#305)
Browse files Browse the repository at this point in the history
Also fixing indexing error in get_swt_formfactor.
  • Loading branch information
kbarros authored Aug 30, 2024
1 parent 4093635 commit 50075d5
Show file tree
Hide file tree
Showing 17 changed files with 237 additions and 210 deletions.
84 changes: 42 additions & 42 deletions examples/01_LSWT_CoRh2O4.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
# # 1. Spin wave simulations of CoRh₂O₄
#
# This tutorial introduces Sunny through its features for performing
# conventional spin wave theory calculations. For concreteness, we consider the
# crystal CoRh₂O₄ and reproduce the calculations of [Ge et al., Phys. Rev. B 96,
# conventional spin wave theory calculations. We consider the crystal CoRh₂O₄
# and reproduce the calculations of [Ge et al., Phys. Rev. B 96,
# 064413](https://doi.org/10.1103/PhysRevB.96.064413).

# ### Get Julia and Sunny
Expand All @@ -13,29 +13,28 @@
# 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 both `Sunny` and `GLMakie`. The latter is needed
# From the Julia prompt, load Sunny and also [GLMakie](https://docs.makie.org/)
# for graphics.

using Sunny, GLMakie

# If these packages are not yet installed, Julia will offer to install them for
# you. If executing this script gives an error, you may need to `update` Sunny
# and GLMakie from the [built-in package
# 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

# The [`Units`](@ref Sunny.Units) object selects reference energy and length
# scales. This is achieved by providing physical constants. For example,
# `units.K` would provide one kelvin as 0.086 meV, with the Boltzmann constant
# implicit.
# 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 can be specified
# from atom positions and types.
# 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
Expand Down Expand Up @@ -72,13 +71,13 @@ view_crystal(cryst)

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

# Previous work demonstrated that inelastic neutron scattering data for CoRh₂O₄
# is well described with a single antiferromagnetic nearest neighbor exchange,
# `J = 0.63` meV. Use [`set_exchange!`](@ref) with the bond that connects atom 1
# to atom 3, and has zero displacement between chemical cells. Applying the
# symmetries of spacegroup 227, Sunny will propagate this interaction to the
# other nearest-neighbor bonds. Calling [`view_crystal`](@ref) with `sys` now
# shows the antiferromagnetic Heisenberg interactions as blue polkadot spheres.
# 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(1, 3, [0, 0, 0]))
Expand All @@ -105,21 +104,21 @@ plot_spins(sys; color=[S[3] for S in sys.dipoles])

# ### Reshaping the magnetic cell

# The same Néel order can also be described with a magnetic cell that consists
# of the 2 cobalt atoms in the primitive cell. Columns of the 3×3 `shape` matrix
# below are the primitive lattice vectors in units of the conventional, cubic
# lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. Use [`reshape_supercell`](@ref) to
# construct a system with this shape, and verify that the energy per site is
# unchanged.
# The most compact magnetic cell for this Néel order is a primitive unit cell.
# Reduce the magnetic cell size using [`reshape_supercell`](@ref), where columns
# of the `shape` matrix are primitive lattice vectors as multiples of the
# conventional cubic lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. One could
# alternatively use `shape = cryst.latvecs \ cryst.prim_latvecs`. Verify that
# the energy per site is unchanged after the reshaping the supercell.

shape = [0 1 1;
1 0 1;
1 1 0] / 2
sys_prim = reshape_supercell(sys, shape)
@assert energy_per_site(sys_prim) -2J*(3/2)^2

# Plotting the spins of `sys_prim` shows the primitive cell as a gray wireframe
# inside the conventional cubic cell.
# 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])

Expand Down Expand Up @@ -150,20 +149,21 @@ kernel = lorentzian(fwhm=0.8)
qs = [[0, 0, 0], [1/2, 0, 0], [1/2, 1/2, 0], [0, 0, 0]]
path = q_space_path(cryst, qs, 500)

# Calculate the single-crystal scattering [`intensities`](@ref)` along the path,
# for 300 energy points between 0 and 6 meV. Use [`plot_intensities`](@ref) to
# visualize the result.
# 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)

# To directly compare with the available experimental data, perform a
# [`powder_average`](@ref) over all possible crystal orientations. Consider 200
# ``𝐪`` magnitudes ranging from 0 to 3 inverse angstroms. Each magnitude
# defines spherical shell in reciprocal space, to be sampled with `2000`
# ``𝐪``-points. The calculation completes in just a couple seconds because the
# magnetic cell size is small.
# 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
Expand All @@ -180,14 +180,14 @@ plot_intensities(res; units, saturation=1.0)

# ### What's next?
#
# * For more spin wave calculations of this traditional type, one can browse the
# [SpinW tutorials ported to Sunny](@ref "SW01 - FM Heisenberg chain").
# * Spin wave theory neglects thermal fluctuations of the magnetic order. Our
# [next 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.
# * 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 our [tutorial on FeI₂](@ref "3. Multi-flavor spin wave
# demonstrated in the [FeI₂ tutorial](@ref "3. Multi-flavor spin wave
# simulations of FeI₂").
45 changes: 23 additions & 22 deletions examples/03_LSWT_SU3_FeI2.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
# # 3. Multi-flavor spin wave simulations of FeI₂
#
# This tutorial illustrates various powerful features in Sunny, including
# symmetry analysis, energy minimization, and spin wave theory with multi-flavor
# bosons.
# 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
Expand Down Expand Up @@ -42,10 +41,10 @@ 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)

# Observe that Sunny inferred the space group, 'P -3 m 1' (164) and labeled the
# atoms according to their point group symmetries. Only the Fe atoms are
# magnetic, so we focus on them with [`subcrystal`](@ref). Importantly, this
# operation preserves the spacegroup symmetries.
# 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`](@ref). Importantly, this operation preserves the
# spacegroup symmetries.

cryst = subcrystal(cryst, "Fe")
view_crystal(cryst)
Expand Down Expand Up @@ -143,14 +142,15 @@ set_onsite_coupling!(sys, S -> -D*S[3]^2, 1)

# ### Finding the ground state

# This model has been carefully designed so that energy minimization produces
# the physically correct magnetic ordering. Using [`set_dipole!`](@ref), this
# magnetic structure can be entered manually. Sunny also provides tools to
# search for an unknown magnetic order, as we will now demonstrate.
# These model parameters have already been fitted so that energy minimization
# yields the physically correct ground state. Knowing this, we could set the
# magnetic configuration manually by calling [`set_dipole!`](@ref) on each site
# in the system. Another approach, as we now demonstrate, is to use the built-in
# optimization tools to search for the ground-state in an automated way.
#
# To minimize bias in the search, use [`resize_supercell`](@ref) to create a
# relatively large system of 4×4×4 chemical cells. Randomize all spins (as SU(3)
# coherent states) and minimize the energy.
# relatively large system of 4×4×4 chemical cells. Randomize all spins
# (represented as SU(3) coherent states) and minimize the energy.

sys = resize_supercell(sys, (4, 4, 4))
randomize_spins!(sys)
Expand All @@ -164,11 +164,12 @@ minimize_energy!(sys)
plot_spins(sys; color=[S[3] for S in sys.dipoles])

# To better understand the spin configuration, we could inspect the static
# structure factor ``\mathcal{S}(𝐪)`` in the 3D space of momenta ``𝐪``. For
# this, Sunny provides [`SampledCorrelationsStatic`](@ref). Here, however, we
# will use [`print_wrapped_intensities`](@ref), which gives static intensities
# averaged over the individual Bravais sublattices (in effect, all ``𝐪``
# intensities are periodically wrapped to the first Brillouin zone).
# structure factor ``\mathcal{S}(𝐪)`` in the 3D space of momenta ``𝐪``. The
# general tool for this analysis is [`SampledCorrelationsStatic`](@ref). For the
# present purposes, however, it is more convenient to use
# [`print_wrapped_intensities`](@ref), which reports ``\mathcal{S}(𝐪)`` with
# periodical wrapping of all commensurate ``𝐪`` wavevectors into the first
# Brillouin zone.

print_wrapped_intensities(sys)

Expand Down Expand Up @@ -204,14 +205,14 @@ plot_spins(sys_min; color=[S[3] for S in sys_min.dipoles], ghost_radius=12)
# SU(3) coherent states, this calculation will require a multi-flavor boson
# generalization of the usual spin wave theory.

swt = SpinWaveTheory(sys_min; measure=ssf_perp(sys_min))

# Calculate and plot the spectrum along a momentum-space path that connects a
# sequence of high-symmetry ``𝐪``-points. These Sunny commands are identical to
# those described in [`our previous CoRh₂O₄ tutorial`](@ref "1. Spin wave
# simulations of CoRh₂O₄").
# sequence of high-symmetry ``𝐪``-points. This interface abstracts over the
# underlying calculator type.

qs = [[0,0,0], [1,0,0], [0,1,0], [1/2,0,0], [0,1,0], [0,0,0]]
path = q_space_path(cryst, qs, 500)
swt = SpinWaveTheory(sys_min; measure=ssf_perp(sys_min))
res = intensities_bands(swt, path)
plot_intensities(res; units)

Expand Down
65 changes: 33 additions & 32 deletions examples/04_GSD_FeI2.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

# 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. Here we perform an analogous calculation at finite
# 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).
#
Expand All @@ -13,14 +13,14 @@
# enables the study of [highly non-equilibrium
# processes](https://doi.org/10.1103/PhysRevB.106.235154).
#
# The structure of this tutorial largely follows our [previous study of CoRh₂O₄
# 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₂ has a strong easy-axis anisotropy that introduces a single-ion
# bound state and necessitates the use of `:SUN` mode.
# whereas FeI₂ is best modeled using `:SUN` mode, owing to its strong easy-axis
# anisotropy.
#
# Construct the FeI₂ system as described in the [previous tutorial](@ref "3.
# Multi-flavor spin wave simulations of FeI₂").
# Construct the FeI₂ system as [previously described](@ref "3. Multi-flavor spin
# wave simulations of FeI₂").

using Sunny, GLMakie

Expand Down Expand Up @@ -64,10 +64,12 @@ set_onsite_coupling!(sys, S -> -D*S[3]^2, 1)

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

# Previously, we used [`minimize_energy!`](@ref) to find a local energy minimum.
# Here, we will instead use [`Langevin`](@ref) dynamics to relax the system into
# thermal equilibrium. The temperature 2.3 K ≈ 0.2 meV is within the ordered
# phase, but large enough so that the dynamics can overcome local energy
# 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)
Expand All @@ -89,42 +91,42 @@ for _ in 1:10_000
end
plot_spins(sys; color=[S[3] for S in sys.dipoles])

# The two-up, two-down spiral order can be verified by calling
# Verify the expected two-up, two-down spiral magnetic order by calling
# [`print_wrapped_intensities`](@ref). A single propagation wavevector ``±𝐤``
# provides most of the static intensity in ``\mathcal{S}(𝐪)``. A smaller amount
# of intensity is spread among many other wavevectors due to thermal
# fluctuations.
# 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)

# Calling [`suggest_timestep`](@ref) shows that thermalization has not
# substantially altered the suggested `dt`.
# Thermalization has not substantially altered the suggested `dt`.

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

# ### Structure factor in the paramagnetic phase

# Now we will re-thermalize the system to a temperature of 3.5 K ≈ 0.30 meV,
# which is 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

# At this higher temperature, the suggested timestep has increased slightly.
# 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;

# Now collect dynamical spin structure factor data using a
# [`SampledCorrelations`](@ref) object. This will involve sampling spin
# configurations from thermal equilibrium and integrating a [classical spin
# dynamics for SU(_N_) coherent states](https://arxiv.org/abs/2204.07563).
# Normal modes appearing in the classical dynamics can be quantized to yield
# magnetic excitations. The associated structure factor intensities
# ``S^{αβ}(q,ω)`` can be compared with inelastic neutron scattering data.
# Incorporate the [`FormFactor`](@ref) appropriate to Fe²⁺.
# 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)
Expand All @@ -147,11 +149,10 @@ for _ in 1:2
add_sample!(sc, sys)
end

# Measure intensities along a path connecting high-symmetry ``𝐪``-points,
# specified 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.
# 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)")
Expand Down
Loading

0 comments on commit 50075d5

Please sign in to comment.