Skip to content

Commit

Permalink
Add SpinW ports to docs
Browse files Browse the repository at this point in the history
  • Loading branch information
kbarros committed Sep 18, 2023
1 parent d7ca521 commit 78288a0
Show file tree
Hide file tree
Showing 2 changed files with 167 additions and 0 deletions.
64 changes: 64 additions & 0 deletions examples/spinw_ports/08_Kagome_AFM.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
## Kagome Antiferromagnet
#
# - Sunny port of the SpinW tutorial authored by Bjorn Fak and Sandor Toth,
# https://spinw.org/tutorials/08tutorial.
# - Authors: Harry Lane
# - Goal: Calculate the linear spin wave theory spectrum for the $\sqrt{3}
# \times \sqrt{3}$ order of a Kagome antiferromagnet.

# Load Packages

using Sunny, GLMakie

# Build a [`Crystal`](@ref) with $P\overline{3}$ space group and Cr$^{+}$ ions
# on each site.

a = b = 6.0 # (Å)
c = 40.0
latvecs = lattice_vectors(a, b, c, 90, 90, 120)
crystal = Crystal(latvecs, [[1/2,0,0]], 147; types=["Cr"])

# Build a [`System`](@ref) with antiferrogmanetic nearest neighbor exchange
# $J=1$.

S = 1
sys = System(crystal, (3,3,1), [SpinInfo(1; S, g=2)], :dipole)
J = 1.0
set_exchange!(sys, J, Bond(2,3,[0,0,0]))

# Initialize to the known magnetic structure, which is 120° order.

q = -[1/3, 1/3, 0]
axis = [0,0,1]
set_spiral_order_on_sublattice!(sys, 1; q, axis, S0=[cos(0),sin(0),0])
set_spiral_order_on_sublattice!(sys, 2; q, axis, S0=[cos(0),sin(0),0])
set_spiral_order_on_sublattice!(sys, 3; q, axis, S0=[cos(2π/3),sin(2π/3),0])
plot_spins(sys; ghost_radius=30, orthographic=true)

# Check energy. Each site participates in 4 bonds with energy J*S^2*cos(2π/3).
# Factor of 1/2 avoids double counting.

@assert energy_per_site(sys) (4/2)*J*S^2*cos(2π/3)

# Define a path in reciprocal space.

points_rlu = [[-1/2, 0, 0], [0, 0, 0], [1/2, 1/2, 0]]
density = 100
path, xticks = reciprocal_space_path(crystal, points_rlu, density)

# Calculate discrete intensities

swt = SpinWaveTheory(sys)
formula = intensity_formula(swt, :perp; kernel=delta_function_kernel)
disp, intensity = intensities_bands(swt, path, formula)

# Plot over a restricted color range from [0,1e-2]. Note that the intensities of
# the flat band at zero-energy are divergent.

fig = Figure()
ax = Axis(fig[1,1]; xlabel="Momentum (r.l.u.)", ylabel="Energy (meV)", xticks, xticklabelrotation=π/6)
ylims!(ax, -1e-1, 2.3)
for i in axes(disp)[2]
lines!(ax, 1:length(disp[:,i]), disp[:,i]; color=intensity[:,i], colorrange=(0,1e-2))
end
fig
103 changes: 103 additions & 0 deletions examples/spinw_ports/15_Ba3NbFe3Si2O14.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,103 @@
## Ba<sub>3</sub>NbFe<sub>3</sub>Si<sub>2</sub>O<sub>14</sub>
#
# - Sunny port of the SpinW tutorial authored by Toth et al.,
# https://spinw.org/tutorials/15tutorial.
# - Authors: Harry Lane
# - Goal: Calculate the linear spin wave theory spectrum for
# Ba<sub>3</sub>NbFe<sub>3</sub>Si<sub>2</sub>O<sub>14</sub>.

# Load Packages

using Sunny, GLMakie

# Build a [`Crystal`](@ref) for
# Ba<sub>3</sub>NbFe<sub>3</sub>Si<sub>2</sub>O<sub>14</sub> using the crystal
# structure from [Marty et al., Phys. Rev. Lett. **101**, 247201
# (2008)](http://dx.doi.org/10.1103/PhysRevLett.101.247201).

a = b = 8.539 # (Å)
c = 5.2414
latvecs = lattice_vectors(a, b, c, 90, 90, 120)
types = ["Fe","Nb","Ba","Si","O","O","O"]
positions = [[0.24964,0,0.5],[0,0,0],[0.56598,0,0],[2/3,1/3,0.5220],[2/3,1/3,0.2162],[0.5259,0.7024,0.3536],[0.7840,0.9002,0.7760]]
langasite = Crystal(latvecs, positions, 150; types)
crystal = subcrystal(langasite, "Fe")
view_crystal(crystal, 7)

# Create a [`System`](@ref) with a lattice size of $(1,1,7)$. The magnetic
# structure of Ba<sub>3</sub>NbFe<sub>3</sub>Si<sub>2</sub>O<sub>14</sub> was
# determined to have the ordering wavevector $𝐐=(0,0,1/7)$ and hence the
# magnetic unit cell has 7 sites. By passing an explicit `seed`, the system's
# random number generator will give repeatable results.

latsize = (1,1,7)
S = 5/2
seed = 5
sys = System(crystal, latsize, [SpinInfo(1; S, g=2)], :dipole; seed)

# Set exchange interactions as parametrized in [Loire et al., Phys. Rev. Lett.
# **106**, 207201 (2011)](http://dx.doi.org/10.1103/PhysRevLett.106.207201)

J₁ = 0.85
J₂ = 0.24
J₃ = 0.053
J₄ = 0.017
J₅ = 0.24
set_exchange!(sys, J₁, Bond(3, 2, [1,1,0]))
set_exchange!(sys, J₄, Bond(1, 1, [0,0,1]))
set_exchange!(sys, J₂, Bond(1, 3, [0,0,0]))

# The final two exchanges define the chirality of the magnetic structure. The
# crystal chirality, $\epsilon_T$, the chirality of each triangle, $ϵ_D$ and the
# sense of rotation of the spin helices along $c$, $ϵ_{H}$. The three
# chiralities are related by $ϵ_T=ϵ_D ϵ_H$. We now assign $J_3$ and $J_5$
# according to the crystal chirality.

ϵD = -1
ϵH = +1
ϵT = ϵD * ϵH

if ϵT == -1
set_exchange!(sys, J₃, Bond(2, 3, [-1,-1,1]))
set_exchange!(sys, J₅, Bond(3, 2, [1,1,1]))
elseif ϵT == 1
set_exchange!(sys, J₅, Bond(2, 3, [-1,-1,1]))
set_exchange!(sys, J₃, Bond(3, 2, [1,1,1]))
else
throw("Provide a valid chirality")
end

# Whilst Sunny provides tools to optimize the ground state automatically, in
# this case we already know the model ground state. Set the spiral magnetic
# order using [`set_spiral_order_on_sublattice!`](@ref). It takes an ordering
# wavevector `q`, an axis of rotation for the spins `axis`, and the initial spin
# `S0` for each sublattice.

q = [0, 0, 1/7]
axis = [0,0,1]
set_spiral_order_on_sublattice!(sys, 1; q, axis, S0=[1, 0, 0])
set_spiral_order_on_sublattice!(sys, 2; q, axis, S0=[-1/2, -sqrt(3)/2, 0])
set_spiral_order_on_sublattice!(sys, 3; q, axis, S0=[-1/2, +sqrt(3)/2, 0])

plot_spins(sys; color=[s[1] for s in sys.dipoles])

# Define a path in reciprocal space, $[0,1,-1+\xi]$ for $\xi = 0 \dots 3$.

points_rlu = [[0,1,-1],[0,1,-1+1],[0,1,-1+2],[0,1,-1+3]];
density = 100
path, xticks = reciprocal_space_path(crystal, points_rlu, density);

# Calculate broadened intensities

swt = SpinWaveTheory(sys)
γ = 0.15 # width in meV
broadened_formula = intensity_formula(swt, :perp; kernel=lorentzian(γ))
energies = collect(0:0.01:6) # 0 < ω < 6 (meV).
is = intensities_broadened(swt, path, energies, broadened_formula)

# Plot

fig = Figure()
ax = Axis(fig[1,1]; xlabel="Momentum (r.l.u.)", ylabel="Energy (meV)", xticks, xticklabelrotation=π/6)
heatmap!(ax, 1:size(is,1), energies, is, colorrange=(0,5))
fig

0 comments on commit 78288a0

Please sign in to comment.