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 ("C 1 2/m 1" setting).

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, "C 1 2/m 1"; types)
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)
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.5069821896015302 0.4499934111736808 -0.735170041097373; 0.5498268278198403 0.48802197566983113 0.6778827410938318; 0.6638219029501892 -0.7478906879775881 -6.698358364433462e-12], [0.019418397373178128 0.06954879691804414 0.9973895380896625; 0.26821812771512105 0.960648180478007 -0.07220878969550726; -0.9631624794438645 0.2689201334812022 1.0145501275969563e-11], [0.6967909876676388 -0.5971877654622452 0.39730226816287784; -0.30166755886613383 0.25854550124710163 0.9176878051454279; -0.6507526438155317 -0.7592897974865037 -2.2072404442720218e-10], [0.9564360149358555 -0.13085216474332595 0.2609748269768682; 0.25856617764781703 -0.035375020737857435 -0.9653456063423058; 0.13554955223171467 0.9907705682395809 2.1183461113370849e-10], [0.5642586650401666 -0.7237147688456205 0.3973022681558963; -0.2442892303439594 0.3133238966625044 0.9176878051484505; -0.7886285125869121 -0.6148699611584183 -2.1388529264011614e-10], [0.9113746600812773 -0.3182583363250159 0.2609748269879097; 0.24638413718438887 -0.08603904500442602 -0.9653456063393206; 0.3296833115364509 0.9440915814127142 2.025515923356077e-10]], StaticArraysCore.SVector{3, Float64}[[0.5069821896015302, 0.4499934111736808, -0.735170041097373] [0.019418397373178128, 0.06954879691804414, 0.9973895380896625] … [0.5642586650401666, -0.7237147688456205, 0.3973022681558963] [0.9113746600812773, -0.3182583363250159, 0.2609748269879097]; [0.5498268278198403, 0.48802197566983113, 0.6778827410938318] [0.26821812771512105, 0.960648180478007, -0.07220878969550726] … [-0.2442892303439594, 0.3133238966625044, 0.9176878051484505] [0.24638413718438887, -0.08603904500442602, -0.9653456063393206]; [0.6638219029501892, -0.7478906879775881, -6.698358364433462e-12] [-0.9631624794438645, 0.2689201334812022, 1.0145501275969563e-11] … [-0.7886285125869121, -0.6148699611584183, -2.1388529264011614e-10] [0.3296833115364509, 0.9440915814127142, 2.025515923356077e-10]], 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.2140975054602642, 2.138207278835357e-12, 0.8929512472478204], [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)