Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Update FeI2 doc #112

Merged
merged 2 commits into from
Aug 1, 2023
Merged

Update FeI2 doc #112

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d6d8ccd
Update FeI2 doc
kbarros Jul 30, 2023
78097f7
Tweaks
kbarros Aug 1, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
346 changes: 84 additions & 262 deletions examples/fei2_tutorial.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,84 +157,94 @@ set_onsite_coupling!(sys, -D*S[3]^2, 1)
# operators with [`print_stevens_expansion`](@ref).

# # Calculating structure factor intensities

# In the remainder of this tutorial, we will examine Sunny's tools for
# calculating structure factors using generalized SU(_N_) classical dynamics.
# This will involve the sampling of spin configurations from the Boltzmann
# distribution at a finite temperature. Spin-spin correlations measured
# dynamically then provide an estimate of the structure factor
# $\mathcal{S}^{\alpha\beta}(\mathbf{q},\omega)$.
#
# ## Simulated annealing
# calculating the dynamical structure factor using a generalization of linear
# spin wave theory that captures quasi-particle excitations that include coupled
# dipolar and quadrupolar fluctuations.
#
# The [Langevin dynamics of SU(_N_) coherent
# states](https://arxiv.org/abs/2209.01265) can be used to sample spin
# configurations in thermal equlibrium. Begin by constructing a
# [`Langevin`](@ref) object.

Δt = 0.05/D # Single-ion anisotropy is the strongest interaction, so 1/D is
# a natural dynamical time-scale (in units with ħ=1).
λ = 0.1 # Dimensionless magnitude of coupling to thermal bath
langevin = Langevin(Δt; kT=0, λ);

# Attempt to find a low-energy spin configuration by lowering the temperature kT
# from 2 to 0 using 20,000 Langevin time-steps.
# A review of the quantum theory is provided in https://arxiv.org/abs/1307.7731.
# The same theory can alternatively be viewed as the quantization of a classical
# dynamics on SU(_N_) coherent states, https://arxiv.org/abs/2106.14125.
# Although outside the scope of this tutorial, Sunny provides support for
# simulating the classical SU(_N_) spin dynamics via a [`Langevin`](@ref)
# equation that generalizes the stochastic Landau-Lifshitz equation,
# https://arxiv.org/abs/2209.01265. This Langevin equation can be used to study
# equilibrium thermal fluctuations of classical SU(_N_) coherent states, or
# their non-equilibrium dynamics. Examples of the latter include relaxation of a
# spin glass, or the driven-dissipative dynamics of topological magnetic
# textures. To learn more, see the other Sunny documentation examples.

# ## Finding the ground state

# Begin with a random configuration and use [`minimize_energy!`](@ref) to find a
# configuration of the SU(3) coherent states (i.e. spin dipoles and quadrupoles)
# that locally minimizes energy.

randomize_spins!(sys)
for kT in range(2, 0, 20_000)
langevin.kT = kT
step!(sys, langevin)
end
minimize_energy!(sys);

# Because the quench was relatively fast, it is expected to find defects in the
# magnetic order. These can be visualized.
# The expected ground state for FeI$_2$ is an antiferrogmanetic striped phase
# with a period of four spins (two up, two down). Visualizing the result of
# optimization, however, may indicate the system got stuck in a local minimum
# with defects.

plot_spins(sys; arrowlength=2.5, linewidth=0.75, arrowsize=1.5)

# If we had used a slower annealing procedure, involving 100,000 or more
# Langevin time-steps, it would very likely find the correct ground state.
# Instead, for purposes of illustration, let's analyze the imperfect spin
# configuration currently stored in `sys`.
#
# An experimental probe of magnetic order order is the 'instantaneous' or
# 'static' structure factor intensity, available via
# [`instant_correlations`](@ref) and related functions. To infer periodicities
# of the magnetic supercell, however, it is sufficient to look at the structure
# factor weights of spin sublattices individually, _without_ phase averaging.
# This information is provided by [`print_wrapped_intensities`](@ref) (see the
# API documentation for a physical interpretation).
# A better understanding of the magnetic ordering can often be obtained by
# moving to Fourier space. The 'instant' structure factor $𝒮(𝐪)$ is an
# experimental observable. To investigate $𝒮(𝐪)$ as true 3D data, Sunny
# provides [`instant_correlations`](@ref) and related functions. Here, however,
# we will use the lighter weight function [`print_wrapped_intensities`](@ref) to
# get a quick understanding of the periodicities present in the spin
# configuration.

print_wrapped_intensities(sys)

# The above is consistent with known results. The zero-field energy-minimizing
# magnetic structure of FeI$_2$ is single-$q$. If annealing were perfect, then
# spontaneous symmetry breaking would select one of $±q = [0, -1/4, 1/4]$,
# $[1/4, 0, 1/4]$, or $[-1/4,1/4,1/4]$.

# Let's break the symmetry by hand. Given a list of $q$ modes, Sunny can suggest
# a magnetic supercell with appropriate periodicity. The result is printed in
# units of the crystal lattice vectors.
# The precise output may vary with Sunny version due to, e.g., different
# floating point roundoff effects. Very likely, however, the result will be
# approximately consistent with the known zero-field energy-minimizing magnetic
# structure of FeI$_2$, which is single-$Q$. Mathematically, spontaneous
# symmetry breaking should select one of $±Q = [0, -1/4, 1/4]$, $[1/4, 0, 1/4]$,
# or $[-1/4,1/4,1/4]$, associated with the three-fold rotational symmetry of the
# crystal spacegroup. In practice, however, one will frequently encounter
# competing "domains" associated with the three possible orientations of the
# ground state.

# If the desired ground state is already known, as with FeI$_2$, it could be
# entered by hand using [`polarize_spin!`](@ref). Alternatively, in the case of
# FeI$_2$, we could repeatedly employ the above randomization and minimization
# procedure until a defect-free configuration is found. Some systems will have
# more complicated ground states, which can be much more challenging to find.
# For this, Sunny provides experimental support for powerful simulated annealing
# via [parallel tempering](https://en.wikipedia.org/wiki/Parallel_tempering),
# but that is outside the scope of this tutorial.

# Here, let's break the three-fold symmetry of FeI$_2$ by hand. Given one or
# more desired $Q$ modes, Sunny can suggest a magnetic supercell with
# appropriate periodicity. Let's arbitrarily select one of the three possible
# ordering wavevectors, $Q = [0, -1/4, 1/4]$. Sunny suggest a corresponding
# magnetic supercell in units of the crystal lattice vectors.

suggest_magnetic_supercell([[0, -1/4, 1/4]], sys.latsize)

# The function [`reshape_geometry`](@ref) allows an arbitrary reshaping of the
# system. After selecting the supercell geometry, it becomes much easier to find
# the energy-minimizing spin configuration.
# system's supercell. We select the supercell appropriate to the broken-symmetry
# ground-state, which makes optimization much easier.

sys_supercell = reshape_geometry(sys, [1 0 0; 0 1 -2; 0 1 2])

langevin.kT = 0
for i in 1:10_000
step!(sys_supercell, langevin)
end

plot_spins(sys_supercell; arrowlength=2.5, linewidth=0.75, arrowsize=1.5)
sys_min = reshape_geometry(sys, [1 0 0; 0 1 -2; 0 1 2])
randomize_spins!(sys_min)
minimize_energy!(sys_min)
plot_spins(sys_min; arrowlength=2.5, linewidth=0.75, arrowsize=1.5)

# ## Linear spin wave theory
#
# Now that we have found the ground state for a magnetic supercell, we can
# immediately proceed to perform zero-temperature calculations using linear spin
# wave theory. We begin by instantiating a `SpinWaveTheory` type using the
# supercell.

swt = SpinWaveTheory(sys_supercell);
swt = SpinWaveTheory(sys_min);

# The dispersion relation can be calculated by providing `dispersion` with a
# `SpinWaveTheory` and a list of wave vectors. For each wave vector,
Expand Down Expand Up @@ -267,12 +277,11 @@ formula = intensity_formula(swt; kernel = delta_function_kernel)
# The `delta_function_kernel` specifies that we want the energy and intensity of each
# band individually.

disp, intensity = intensities_bands(swt, path; formula = formula)
disp, intensity = intensities_bands(swt, path; formula)

fig = Figure()
ax = Axis(fig[1,1]; xlabel="𝐪", ylabel="Energy (meV)",
xticks=(markers, labels), xticklabelrotation=π/6,
)
xticks=(markers, labels), xticklabelrotation=π/6)
ylims!(ax, 0.0, 7.5)
xlims!(ax, 1, size(disp, 1))
for i in axes(disp)[2]
Expand All @@ -291,11 +300,11 @@ broadened_formula = intensity_formula(swt; kernel = lorentzian(γ))
energies = collect(0.0:0.01:7.5) # Energies to calculate
is = intensities_broadened(swt, path, energies, broadened_formula)

heatmap(1:size(is, 1), energies, is;
axis=(xlabel = "(H,0,0)", ylabel="Energy (meV)", xticks=(markers, labels),
xticklabelrotation=π/8,
),
)
fig = Figure()
ax = Axis(fig[1,1], xlabel="(H,0,0)", ylabel="Energy (meV)", xticks=(markers, labels),
xticklabelrotation=π/8)
heatmap!(ax, 1:size(is, 1), energies, is)
fig

# The existence of a lower-energy, single-ion bound state is in qualitative
# agreement with the experimental data in [Bai et
Expand All @@ -315,203 +324,16 @@ disp, Sαβ = dssf(swt, [[0, 0, 0]]);

# `disp` is identical to the output that is obtained from `dispersion` and
# contains the energy of each mode at the specified wave vectors. `Sαβ` contains
# a 3x3 matrix for each of these modes. The matrix elements of `Sαβ` correspond
# to correlations of the different spin components (ordered x, y, z). For
# example, the full set of matrix elements for the first mode may be obtained as
# follows,

Sαβ[1]

# and the `xx` matrix element is

Sαβ[1][1,1]

# ## Dynamical structure factors with classical dynamics
# Linear spin wave calculations are very useful for getting quick, high-quality,
# results at zero temperature. Moreover, these results are obtained in the
# thermodynamic limit. Classical dynamics may also be used to produce similar
# results, albeit at a higher computational cost and on a finite sized lattice.
# The classical approach nonetheless provides a number of complementary
# advantages: it is possible perform simulations at finite temperature while
# retaining nonlinearities; out-of-equilibrium behavior may be examined
# directly; and it is straightforward to incorporate inhomogenties, chemical or
# otherwise.

# Because classical simulations are conducted on a finite-sized lattice,
# obtaining acceptable resolution in momentum space requires the use of a larger
# system size. We can now resize the magnetic supercell to a much larger
# simulation volume, provided as multiples of the original unit cell.

sys_large = resize_periodically(sys_supercell, (16,16,4))
plot_spins(sys_large; arrowlength=2.5, linewidth=0.75, arrowsize=1.5)

# Apply Langevin dynamics to thermalize the system to a target temperature.

kT = 0.5 * meV_per_K # 0.5K in units of meV
langevin.kT = kT

for _ in 1:10_000
step!(sys_large, langevin)
end

# To estimate the dynamic structure factor, we can collect spin-spin correlation
# data by first generating an initial condition at thermal equilibrium and then
# integrating the [Hamiltonian dynamics of SU(_N_) coherent
# states](https://arxiv.org/abs/2204.07563). Samples are accumulated into a
# `SampledCorrelations`, which is initialized by calling
# [`dynamical_correlations`](@ref). `dynamical_correlations` takes a `System`
# and three keyword parameters that determine the ω information that will be
# available: an integration step size, the number of ωs to resolve, and the
# maximum ω to resolve. For the time step, twice the value used for the Langevin
# integrator is usually a good choice.

sc = dynamical_correlations(sys_large; Δt=2Δt, nω=120, ωmax=7.5)
# a $3×3$ matrix for each of these modes. For example, `Sαβ[1]` are the
# dynamical structure factor intensities for the first mode, and `Sαβ[1][2,3]`
# is associated with the $α=ŷ$ and $β=ẑ$ spin components.

# `sc` currently contains no data. A sample can be accumulated into it by
# calling [`add_sample!`](@ref).

add_sample!(sc, sys_large)

# Additional samples can be added after generating new spin configurations:

for _ in 1:2
for _ in 1:1000 # Fewer steps needed in equilibrium
step!(sys_large, langevin)
end
add_sample!(sc, sys_large) # Accumulate the sample into `sc`
end

# ## Accessing structure factor data
# The basic functions for accessing intensity data are [`intensities_interpolated`](@ref)
# and [`intensities_binned`](@ref). Both functions accept an [`intensity_formula`](@ref)
# to specify how to combine the correlations recorded in the `SampledCorrelations` into
# intensity data. By default, a formula computing the unpolarized intensity is used,
# but alternative formulas can be specified.
# ## What's next?
#
# By way of example, we will use a formula which computes the trace of the structure
# factor and applies a classical-to-quantum temperature-dependent rescaling `kT`.

formula = intensity_formula(sc, :trace; kT = kT)

# Using the formula, we plot single-$q$ slices at (0,0,0) and (π,π,π):

qs = [[0, 0, 0], [0.5, 0.5, 0.5]]
is = intensities_interpolated(sc, qs; interpolation = :round, formula = formula)

fig = lines(ωs(sc), is[1,:]; axis=(xlabel="meV", ylabel="Intensity"), label="(0,0,0)")
lines!(ωs(sc), is[2,:]; label="(π,π,π)")
axislegend()
fig

# For real calculations, one often wants to apply further corrections and more accurate formulas
# Here, we apply [`FormFactor`](@ref) corrections appropriate for `Fe2` magnetic ions,
# and a dipole polarization correction `:perp`.

formfactors = [FormFactor(1, "Fe2"; g_lande=3/2)]
new_formula = intensity_formula(sc, :perp; kT = kT, formfactors = formfactors)

# Frequently one wants to extract energy intensities along lines that connect
# special wave vectors. The function [`connected_path`](@ref) linearly samples
# between provided $q$-points, with a given sample density.

points = [[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]]
density = 40
path, markers = connected_path(sc, points, density);

# Since scattering intensities are only available at a certain discrete ``(Q,\omega)``
# points, the intensity on the path can be calculated by interpolating between these
# discrete points:

is = intensities_interpolated(sc, path;
interpolation = :linear, # Interpolate between available wave vectors
formula = new_formula
)
is = broaden_energy(sc, is, (ω, ω₀)->lorentzian(ω-ω₀, 0.05)) # Add artificial broadening
labels = string.(points)
heatmap(1:size(is,1), ωs(sc), is;
axis = (
ylabel = "meV",
xticks = (markers, labels),
xticklabelrotation=π/8,
xticklabelsize=12,
)
)

# Whereas [`intensities_interpolated`](@ref) either rounds or linearly interpolates
# between the discrete ``(Q,\omega)`` points Sunny calculates correlations at, [`intensities_binned`](@ref)
# performs histogram binning analgous to what is done in experiments.
# The graphs produced by each method are similar.
cut_width = 0.3
density = 15
paramsList, markers, ranges = connected_path_bins(sc,points,density,cut_width)

total_bins = ranges[end][end]
energy_bins = paramsList[1].numbins[4]
is = zeros(Float64,total_bins,energy_bins)
integrated_kernel = integrated_lorentzian(0.05) # Lorentzian broadening
for k in 1:length(paramsList)
h,c = intensities_binned(sc,paramsList[k];formula = new_formula,integrated_kernel = integrated_kernel)
is[ranges[k],:] = h[:,1,1,:] ./ c[:,1,1,:]
end

heatmap(1:size(is,1), ωs(sc), is;
axis = (
ylabel = "meV",
xticks = (markers, labels),
xticklabelrotation=π/8,
xticklabelsize=12,
)
)


# Often it is useful to plot cuts across multiple wave vectors but at a single
# energy.

npoints = 60
qvals = range(-2, 2, length=npoints)
qs = [[a, b, 0] for a in qvals, b in qvals]

is = intensities_interpolated(sc, qs; formula = new_formula,interpolation = :linear);

ωidx = 30
hm = heatmap(is[:,:,ωidx]; axis=(title="ω=$(ωs(sc)[ωidx]) meV", aspect=true))
Colorbar(hm.figure[1,2], hm.plot)
hidedecorations!(hm.axis); hidespines!(hm.axis)
hm

# Note that Brillouin zones appear 'skewed'. This is a consequence of the fact
# that Sunny measures $q$-vectors as multiples of reciprocal lattice vectors,
# which are not orthogonal. It is often useful to express our wave
# vectors in terms of an orthogonal basis, where each basis element is specified
# as a linear combination of reciprocal lattice vectors. For our crystal, with
# reciprocal vectors $a^*$, $b^*$ and $c^*$, we can define an orthogonal basis
# by taking $\hat{a}^* = 0.5(a^* + b^*)$, $\hat{b}^*=a^* - b^*$, and
# $\hat{c}^*=c^*$. Below, we map `qs` to wavevectors `ks` in the new coordinate
# system and get their intensities.

A = [0.5 1 0;
0.5 -1 0;
0 0 1]
ks = [A*q for q in qs]

is_ortho = intensities_interpolated(sc, ks; formula = new_formula, interpolation = :linear)

hm = heatmap(is_ortho[:,:,ωidx]; axis=(title="ω=$(ωs(sc)[ωidx]) meV", aspect=true))
Colorbar(hm.figure[1,2], hm.plot)
hidedecorations!(hm.axis); hidespines!(hm.axis)
hm

# Finally, we note that instantaneous structure factor data, ``𝒮(𝐪)``, can be
# obtained from a dynamic structure factor with [`instant_intensities_interpolated`](@ref).

is_static = instant_intensities_interpolated(sc, ks; formula = new_formula, interpolation = :linear)

hm = heatmap(is_static; axis=(title="Instantaneous Structure Factor", aspect=true))
Colorbar(hm.figure[1,2], hm.plot)
hidedecorations!(hm.axis); hidespines!(hm.axis)
hm
# Sunny provides much more than just linear spin wave theory. Please explore the
# other examples in the documentation for more details. Sunny's capability to
# simulate the classical [`Langevin`](@ref) dynamics of SU(_N_) coherent states
# allows to incorporate finite-temperature spin fluctuations into the dynamical
# structure factor calculation, or to study systems under highly non-equilibrium
# conditions. Sunny also supports the construction of inhomogeneous systems
# (e.g., systems with quenched disorder) via [`to_inhomogeneous`](@ref).
Loading