Download this example as Julia file or Jupyter notebook.

SW18 - Distorted kagome

This is a Sunny port of SpinW Tutorial 18, originally authored by Goran Nilsen and Sandor Toth. This tutorial illustrates spin wave calculations for KCu₃As₂O₇(OD)₃. The Cu ions are arranged in a distorted kagome lattice, and exhibit an incommensurate helical magnetic order, as described in G. J. Nilsen, et al., Phys. Rev. B 89, 140412 (2014). The model follows Toth and Lake, J. Phys.: Condens. Matter 27, 166002 (2015).

using Sunny, GLMakie

Build the distorted kagome crystal, with spacegroup 12.

units = Units(:meV, :angstrom)
latvecs = lattice_vectors(10.2, 5.94, 7.81, 90, 117.7, 90)
positions = [[0, 0, 0], [1/4, 1/4, 0]]
types = ["Cu1", "Cu2"]
cryst = Crystal(latvecs, positions, 12; types, setting="b1")
view_crystal(cryst)

Define the interactions.

moments = [1 => Moment(s=1/2, g=2), 3 => Moment(s=1/2, g=2)]
sys = System(cryst, moments, :dipole, seed=0)
J   = -2
Jp  = -1
Jab = 0.75
Ja  = -J/.66 - Jab
Jip = 0.01
set_exchange!(sys, J, Bond(1, 3, [0, 0, 0]))
set_exchange!(sys, Jp, Bond(3, 5, [0, 0, 0]))
set_exchange!(sys, Ja, Bond(3, 4, [0, 0, 0]))
set_exchange!(sys, Jab, Bond(1, 2, [0, 0, 0]))
set_exchange!(sys, Jip, Bond(3, 4, [0, 0, 1]))

Use minimize_spiral_energy! to optimize the generalized spiral order. This determines the propagation wavevector k, and fits the spin values within the unit cell. One must provide a fixed axis perpendicular to the polarization plane. For this system, all interactions are rotationally invariant, and the axis vector is arbitrary. In other cases, a good axis will frequently be determined from symmetry considerations.

axis = [0, 0, 1]
randomize_spins!(sys)
k = minimize_spiral_energy!(sys, axis; k_guess=randn(3))
plot_spins(sys; ndims=2)

If successful, the optimization process will find one two propagation wavevectors, ±k_ref, with opposite chiralities. In this system, the spiral_energy_per_site is independent of chirality.

k_ref = [0.785902495, 0.0, 0.107048756]
k_ref_alt = [1, 0, 1] - k_ref
@assert isapprox(k, k_ref; atol=1e-6) || isapprox(k, k_ref_alt; atol=1e-6)
@assert spiral_energy_per_site(sys; k, axis) ≈ -0.78338383838

Check the energy with a real-space calculation using a large magnetic cell. First, we must determine a lattice size for which k becomes approximately commensurate.

suggest_magnetic_supercell([k_ref]; tol=1e-3)

Possible magnetic supercell in multiples of lattice vectors:

    [14 0 1; 0 1 0; 0 0 2]

for the rationalized wavevectors:

    [[11/14, 0, 3/28]]

Resize the system as suggested, and perform a real-space calculation. Working with a commensurate wavevector increases the energy slightly. The precise value might vary from run-to-run due to trapping in a local energy minimum.

new_shape = [14 0 1; 0 1 0; 0 0 2]
sys2 = reshape_supercell(sys, new_shape)
randomize_spins!(sys2)
minimize_energy!(sys2)
energy_per_site(sys2)

-0.7791966901275341

Return to the original system (with a single chemical cell) and construct SpinWaveTheorySpiral for calculations on the incommensurate spiral phase.

measure = ssf_perp(sys; apply_g=false)
swt = SpinWaveTheorySpiral(sys; measure, k, axis)

SpinWaveTheorySpiral(SpinWaveTheory(System([Dipole mode], Supercell (1×1×1)×6, Energy per site -0.3888), Sunny.SWTDataDipole(StaticArraysCore.SMatrix{3, 3, Float64, 9}[[0.9020943508786993 -0.3436939983510774 0.26096018395569187; -0.24386058119897325 0.09290981943523284 0.965349564867469; -0.35603061532688157 -0.93447429121938 5.0611908610256205e-9], [0.897789908020667 0.19003454417474283 0.39731618778567274; 0.38870387513544663 0.08227667063771513 -0.9176817786811832; -0.20708109163062025 0.9783237815207554 8.337602332466484e-10], [-0.016304612228017434 -0.07035946213553063 0.9973884427383811; 0.225160731899302 0.9716410605355589 0.07222391772747692; -0.974185200272856 0.22575029472701363 -4.272637459928541e-8], [0.22919386649422163 -0.6379735853755939 -0.7351597621772009; -0.24855522687613898 0.6918671658035171 -0.6778938885073088; 0.9411112956255386 0.33809692285795545 4.1576892212572746e-8], [-0.029957873352754866 -0.06571773195509688 0.9973884426495319; 0.41370778608165526 0.9075403975946613 0.07222391895445592; -0.9099166959450716 0.41479103948904494 -4.2900198887581894e-8], [0.09787947010812502 -0.6707903855978159 -0.7351597567337905; -0.10614788814945 0.7274561799602088 -0.6778938944105576; 0.9895212151201027 0.14438755080419396 4.259922265509661e-8]], StaticArraysCore.SVector{3, Float64}[[0.9020943508786993, -0.3436939983510774, 0.26096018395569187] [0.897789908020667, 0.19003454417474283, 0.39731618778567274] … [-0.029957873352754866, -0.06571773195509688, 0.9973884426495319] [0.09787947010812502, -0.6707903855978159, -0.7351597567337905]; [-0.24386058119897325, 0.09290981943523284, 0.965349564867469] [0.38870387513544663, 0.08227667063771513, -0.9176817786811832] … [0.41370778608165526, 0.9075403975946613, 0.07222391895445592] [-0.10614788814945, 0.7274561799602088, -0.6778938944105576]; [-0.35603061532688157, -0.93447429121938, 5.0611908610256205e-9] [-0.20708109163062025, 0.9783237815207554, 8.337602332466484e-10] … [-0.9099166959450716, 0.41479103948904494, -4.2900198887581894e-8] [0.9895212151201027, 0.14438755080419396, 4.259922265509661e-8]], Sunny.StevensExpansion[Sunny.StevensExpansion(0, [0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), Sunny.StevensExpansion(0, [0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), Sunny.StevensExpansion(0, [0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), Sunny.StevensExpansion(0, [0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), Sunny.StevensExpansion(0, [0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), Sunny.StevensExpansion(0, [0.0], [0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])], [0.7071067811865476, 0.7071067811865476, 0.7071067811865476, 0.7071067811865476, 0.7071067811865476, 0.7071067811865476]), MeasureSpec, 1.0e-8), [0.21409750520425255, 4.571731887113249e-10, 0.8929512469536419], [0.0, 0.0, 1.0])

Plot intensities for a path through $𝐪$-space.

qs = [[0,0,0], [1,0,0]]
path = q_space_path(cryst, qs, 400)
res = intensities_bands(swt, path)
plot_intensities(res; units)

Plot the powder-averaged intensities

radii = range(0, 2, 100) # (1/Å)
energies = range(0, 6, 200)
kernel = gaussian(fwhm=0.05)
res = powder_average(cryst, radii, 400) do qs
    intensities(swt, qs; energies, kernel)
end
plot_intensities(res; units)