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.

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

Add SpinW tutorials 12, 13, 14 #296

Merged
merged 5 commits into from
Aug 23, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
2 changes: 1 addition & 1 deletion examples/02_LSWT_CoRh2O4.jl
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,7 @@ swt = SpinWaveTheory(sys_prim; measure=ssf_perp(sys_prim))
# path that connects high-symmetry points in reciprocal space.

qs = [[0, 0, 0], [1/2, 0, 0], [1/2, 1/2, 0], [0, 0, 0]]
path = q_space_path(cryst, qs, 400)
path = q_space_path(cryst, qs, 500)

# Select [`lorentzian`](@ref) broadening with a full-width at half-maximum
# (FWHM) of 0.8 meV. Use [`ssf_perp`](@ref) to calculate unpolarized scattering
Expand Down
2 changes: 1 addition & 1 deletion examples/03_LLD_CoRh2O4.jl
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,7 @@ qs = [[3/4, 3/4, 0],
[1/4, 1, 1/4],
[ 0, 1, 0],
[ 0, -4, 0]]
qpts = q_space_path(cryst, qs, 1000)
qpts = q_space_path(cryst, qs, 500)

# Calculate ``I(𝐪, ω)`` intensities along this path and plot.

Expand Down
2 changes: 1 addition & 1 deletion examples/07_Dipole_Dipole.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ plot_spins(sys_prim; ghost_radius=8, color=[:red, :blue, :yellow, :purple])

qs = [[0,0,0], [0,1,0], [1,1/2,0], [1/2,1/2,1/2], [3/4,3/4,0], [0,0,0]]
labels = ["Γ", "X", "W", "L", "K", "Γ"]
path = q_space_path(cryst, qs, 400; labels)
path = q_space_path(cryst, qs, 500; labels)

measure = ssf_trace(sys_prim)
swt = SpinWaveTheory(sys_prim; measure)
Expand Down
3 changes: 1 addition & 2 deletions examples/08_Momentum_Conventions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ for _ in 1:nsamples
end
add_sample!(sc, sys)
end
path = q_space_path(cryst, [[0,0,-1/2], [0,0,+1/2]], 300)
path = q_space_path(cryst, [[0,0,-1/2], [0,0,+1/2]], 400)
res1 = intensities(sc, path; energies=:available, kT)

# Calculate the same quantity with linear spin wave theory at ``T = 0``. Because
Expand All @@ -76,7 +76,6 @@ res1 = intensities(sc, path; energies=:available, kT)
sys_small = resize_supercell(sys, (1,1,1))
minimize_energy!(sys_small)
swt = SpinWaveTheory(sys_small; measure=ssf_trace(sys_small))
path = q_space_path(cryst, [[0,0,-1/2], [0,0,+1/2]], 400)
res2 = intensities_bands(swt, path)

# This model system has a single magnon band with dispersion ``ϵ(𝐪) = 1 - D/B
Expand Down
10 changes: 6 additions & 4 deletions examples/spinw_tutorials/SW11_La2CuO4.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
# # SW11 - La₂CuO₄
#
# This is a Sunny port of [SpinW Tutorial
# 11](https://spinw.org/tutorials/15tutorial), originally authored by Sandor
# 11](https://spinw.org/tutorials/11tutorial), originally authored by Sandor
# Toth. It calculates the spin wave spectrum of La₂CuO₄.

# Load packages
Expand All @@ -14,7 +14,9 @@ using Sunny, GLMakie

units = Units(:meV, :angstrom)
latvecs = lattice_vectors(1, 1, 10, 90, 90, 90)
cryst = Crystal(latvecs, [[0, 0, 0]])
positions = [[0, 0, 0]]
types = ["Cu"]
cryst = Crystal(latvecs, positions; types)
view_crystal(cryst; dims=2)

# Build a spin system using the exchange parameters from [R. Coldea, Phys. Rev.
Expand Down Expand Up @@ -45,9 +47,9 @@ energies = range(0, 320, 400)
swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
res = intensities(swt, path; energies, kernel=gaussian(fwhm=35))
res.energies .*= 1.18
plot_intensities(res)
plot_intensities(res; units)

# Plot instantaneous itensities, integrated over ω.

res = intensities_instant(swt, path)
plot_intensities(res; colorrange=(0,20))
plot_intensities(res; colorrange=(0,20), units)
67 changes: 67 additions & 0 deletions examples/spinw_tutorials/SW12_Triangular_easy_plane.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
# # SW12 - Triangular lattice with easy plane
#
# This is a Sunny port of [SpinW Tutorial
# 12](https://spinw.org/tutorials/12tutorial), originally authored by Sandor
# Toth. It calculates the spin wave dispersion of a triangular lattice model
# with antiferromagnetic interactions and easy-plane single-ion anisotropy.

# Load packages

using Sunny, GLMakie

# Build a triangular lattice with arbitrary lattice constant of 3 Å.

latvecs = lattice_vectors(3, 3, 4, 90, 90, 120)
cryst = Crystal(latvecs, [[0, 0, 0]])

# Build a system with exchange +1 meV along nearest neighbor bonds.

S = 3/2
J1 = +1.0
sys = System(cryst, (3,3,1), [SpinInfo(1; S, g=2)], :dipole)
set_exchange!(sys, J1, Bond(1, 1, [1, 0, 0]))

# Set an easy-axis anisotropy operator ``+D S_z^2`` using
# [`set_onsite_coupling!`](@ref). **Important note**: When introducing a
# single-ion anisotropy in `:dipole` mode, Sunny will automatically include a
# classical-to-quantum correction factor, as described in the document
# [Interaction Strength Renormalization](@ref). The present single-ion
# anisotropy is quadratic in the spin operators, so the prescription is to
# rescale the interaction strength as ``D → (1 - 1/2S) D``. For purposes of this
# tutorial, our aim is to reproduce the SpinW result, which requires to "undo"
# Sunny's classical-to-quantum rescaling factor. Another way to achieve the same
# thing is to select system mode `:dipole_large_S` instead of `:dipole`.

undo_classical_to_quantum_rescaling = 1 / (1 - 1/2S)
D = 0.2 * undo_classical_to_quantum_rescaling
set_onsite_coupling!(sys, S -> D*S[3]^2, 1)
randomize_spins!(sys)
minimize_energy!(sys)
plot_spins(sys; dims=2)

# Plot the spin wave spectrum for a path through ``𝐪``-space.

qs = [[0, 0, 0], [1, 1, 0]]
path = q_space_path(cryst, qs, 400)
swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
res = intensities_bands(swt, path)
plot_intensities(res)

# To select a specific linear combination of spin structure factor (SSF)
# components in global Cartesian coordinates, one can use [`ssf_custom`](@ref).
# Here we calculate and plot the real part of ``\mathcal{S}^{zz}(𝐪, ω)``.

measure = ssf_custom(sys) do q, ssf
return real(ssf[3, 3])
end
swt = SpinWaveTheory(sys; measure)
res = intensities_bands(swt, path)
plot_intensities(res)

# It's also possible to get data for the full 3×3 SSF. For example, this is the
# SSF for the 7th energy band, at the 10th ``𝐪``-point along the path.

measure = ssf_custom((q, ssf) -> ssf, sys)
swt = SpinWaveTheory(sys; measure)
res = intensities_bands(swt, path)
res.data[7, 10]
94 changes: 94 additions & 0 deletions examples/spinw_tutorials/SW13_LiNiPO4.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
# # SW13 - LiNiPO₄
#
# This is a Sunny port of [SpinW Tutorial
# 13](https://spinw.org/tutorials/13tutorial), originally authored by Sandor
# Toth. It calculates the spin wave spectrum of LiNiPO₄.

# Load packages

using Sunny, GLMakie

# Build an orthorhombic lattice and populate the Ni atoms according to
# spacegroup 62 (Pnma).

units = Units(:meV, :angstrom)
a = 10.02
b = 5.86
c = 4.68
latvecs = lattice_vectors(a, b, c, 90, 90, 90)
positions = [[1/4, 1/4, 0]]
types = ["Ni"]
cryst = Crystal(latvecs, positions, 62, setting=""; types)
view_crystal(cryst)

# Create a system with exchange parameters taken from [T. Jensen, et al., PRB
# **79**, 092413 (2009)](https://doi.org/10.1103/PhysRevB.79.092413). The
# corrected anisotropy values are taken from the thesis of T. Jensen. The mode
# `:dipole_large_S` avoids a [classical-to-quantum rescaling factor](@ref
# "Interaction Strength Renormalization") of anisotropy strengths, as needed for
# consistency with the original fits.

S = 3/2
sys = System(cryst, (1,1,1), [SpinInfo(1, S=1, g=2)], :dipole_large_S)
Jbc = 1.036
Jb = 0.6701
Jc = -0.0469
Jac = -0.1121
Jab = 0.2977
Da = 0.1969
Db = 0.9097
set_exchange!(sys, Jbc, Bond(2, 3, [0, 0, 0]))
set_exchange!(sys, Jc, Bond(1, 1, [0, 0, -1]))
set_exchange!(sys, Jb, Bond(1, 1, [0, 1, 0]))
set_exchange!(sys, Jab, Bond(1, 2, [0, 0, 0]))
set_exchange!(sys, Jab, Bond(3, 4, [0, 0, 0]))
set_exchange!(sys, Jac, Bond(3, 1, [0, 0, 0]))
set_exchange!(sys, Jac, Bond(4, 2, [0, 0, 0]))
set_onsite_coupling!(sys, S -> Da*S[1]^2 + Db*S[2]^2, 1)

# Energy minimization yields a co-linear order along the ``c`` axis.

randomize_spins!(sys)
minimize_energy!(sys)
plot_spins(sys; color=[s[3] for s in sys.dipoles])

# Calculate the spectrum along path [ξ, 1, 0]

swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
qs = [[0, 1, 0], [2, 1, 0]]
path = q_space_path(cryst, qs, 400)
res = intensities_bands(swt, path)
fig = Figure(size=(768, 300))
plot_intensities!(fig[1, 1], res; units);

# There are two physical bands with nonvanishing intensity. To extract these
# intensity curves, we must filter out the additional bands with zero intensity.
# One way is to sort the data along dimension 1 (the band index) with the
# comparison operator "is intensity less than ``10^{-12}``". Doing so moves the
# physical bands to the end of the array axis. Call the [Makie `lines!`
# function](https://docs.makie.org/stable/reference/plots/lines) to make a
# custom plot.

data_sorted = sort(res.data; dims=1, by= >(1e-12))
ax = Axis(fig[1, 2], xlabel="Momentum (r.l.u.)", ylabel="Intensity",
xticks=res.qpts.xticks, xticklabelrotation=π/6)
lines!(ax, data_sorted[end, :]; label="Lower band")
lines!(ax, data_sorted[end-1, :]; label="Upper band")
axislegend(ax)
fig

# Make the same plots along path [0, 1, ξ]

qs = [[0, 1, 0], [0, 1, 2]]
path = q_space_path(cryst, qs, 400)
res = intensities_bands(swt, path)
fig = Figure(size=(768, 300))
plot_intensities!(fig[1, 1], res; units)

data_sorted = sort(res.data; dims=1, by=x->abs(x)>1e-12)
ax = Axis(fig[1, 2], xlabel="Momentum (r.l.u.)", ylabel="Intensity",
xticks=res.qpts.xticks, xticklabelrotation=π/6)
lines!(ax, data_sorted[end, :]; label="Lower band")
lines!(ax, data_sorted[end-1, :]; label="Upper band")
axislegend(ax)
fig
69 changes: 69 additions & 0 deletions examples/spinw_tutorials/SW14_YVO3.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,69 @@
# # SW14 - YVO₃
#
# This is a Sunny port of [SpinW Tutorial
# 14](https://spinw.org/tutorials/14tutorial), originally authored by Sandor
# Toth. It calculates the spin wave spectrum of YVO₃.

# Load packages

using Sunny, GLMakie

# Build an orthorhombic lattice and populate the V atoms according to the
# pseudocubic unit cell, doubled along the c-axis. The listed spacegroup of the
# YVO₃ chemical cell is international number 62. It has been observed, however,
# that the exchange interactions break this symmetry. For this reason, disable
# all symmetry analysis by selecting spacegroup 1 (P1).

units = Units(:meV, :angstrom)
a = 5.2821 / sqrt(2)
b = 5.6144 / sqrt(2)
c = 7.5283
latvecs = lattice_vectors(a, b, c, 90, 90, 90)
positions = [[0, 0, 0], [0, 0, 1/2]]
types = ["V", "V"]
cryst = Crystal(latvecs, positions, 1; types)

# Create a system following the model of [C. Ulrich, et al. PRL **91**, 257202
# (2003)](https://doi.org/10.1103/PhysRevLett.91.257202). The mode
# `:dipole_large_S` avoids a [classical-to-quantum rescaling factor](@ref
# "Interaction Strength Renormalization") of anisotropy strengths, as needed for
# consistency with the original fits.

sys = System(cryst, (2,2,1), [SpinInfo(1, S=1/2, g=2), SpinInfo(2, S=1/2, g=2)], :dipole_large_S)
Jab = 2.6
Jc = 3.1
δ = 0.35
K1 = 0.90
K2 = 0.97
d = 1.15
Jc1 = [-Jc*(1+δ)+K2 0 -d; 0 -Jc*(1+δ) 0; +d 0 -Jc*(1+δ)]
Jc2 = [-Jc*(1-δ)+K2 0 +d; 0 -Jc*(1-δ) 0; -d 0 -Jc*(1-δ)]
set_exchange!(sys, Jab, Bond(1, 1, [1, 0, 0]))
set_exchange!(sys, Jab, Bond(2, 2, [1, 0, 0]))
set_exchange!(sys, Jc1, Bond(1, 2, [0, 0, 0]))
set_exchange!(sys, Jc2, Bond(2, 1, [0, 0, 1]))
set_exchange!(sys, Jab, Bond(1, 1, [0, 1, 0]))
set_exchange!(sys, Jab, Bond(2, 2, [0, 1, 0]))
set_onsite_coupling!(sys, S -> -K1*S[1]^2, 1)
set_onsite_coupling!(sys, S -> -K1*S[1]^2, 2)

# When using spacegroup P1, there is no symmetry-propagation of interactions
# because all bonds are considered inequivalent. One can visualize the
# interactions in the system by clicking the toggles in the
# [`view_crystal`](@ref) GUI.

view_crystal(sys)

# Energy minimization yields a Néel order with canting

randomize_spins!(sys)
minimize_energy!(sys)
plot_spins(sys)

# Plot the spin wave spectrum along a path

swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
qs = [[0.75, 0.75, 0], [0.5, 0.5, 0], [0.5, 0.5, 1]]
path = q_space_path(cryst, qs, 400)
res = intensities_bands(swt, path)
plot_intensities(res; units)
4 changes: 2 additions & 2 deletions examples/spinw_tutorials/SW18_Distorted_kagome.jl
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ energy_per_site(sys2) # < -0.7834 meV
# Define a path in q-space

qs = [[0,0,0], [1,0,0]]
path = q_space_path(cryst, qs, 512)
path = q_space_path(cryst, qs, 400)

# Calculate intensities for the incommensurate spiral phase using
# [`SpiralSpinWaveTheory`](@ref). It is necessary to provide the original `sys`,
Expand All @@ -93,7 +93,7 @@ plot_intensities(res; units)
radii = range(0, 2, 100) # (1/Å)
energies = range(0, 6, 200)
kernel = gaussian(fwhm=0.05)
res = powder_average(cryst, radii, 200) do qs
res = powder_average(cryst, radii, 400) do qs
intensities(swt, qs; energies, kernel)
end
plot_intensities(res; units)
2 changes: 1 addition & 1 deletion examples/spinw_tutorials/SW19_Different_Ions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ plot_spins(sys)

swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
qs = [[0,0,0], [1,0,0]]
path = q_space_path(cryst, qs, 512)
path = q_space_path(cryst, qs, 400)

# Plot three types of pair correlation intensities

Expand Down
Loading