From b582fb182a704c6ef46301684759052e14b3323f Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Sat, 24 Aug 2024 17:56:37 -0600 Subject: [PATCH] Update tutorials 1 through 4 (#299) Also: * Disable `minimize_energy!` warnings which were unreliable. * Simplify `latsize` printing of system ("reshaped" only relevant if unit cell has changed). --- README.md | 13 +- docs/src/structure-factor.md | 2 +- examples/01_LSWT_CoRh2O4.jl | 192 +++++++++++++++++++++ examples/01_LSWT_SU3_FeI2.jl | 322 ----------------------------------- examples/02_LLD_CoRh2O4.jl | 179 +++++++++++++++++++ examples/02_LSWT_CoRh2O4.jl | 106 ------------ examples/03_LLD_CoRh2O4.jl | 158 ----------------- examples/03_LSWT_SU3_FeI2.jl | 253 +++++++++++++++++++++++++++ examples/04_GSD_FeI2.jl | 297 +++++++++++++------------------- examples/05_MC_Ising.jl | 44 +++-- src/Optimization.jl | 7 +- src/System/System.jl | 7 +- 12 files changed, 778 insertions(+), 802 deletions(-) create mode 100644 examples/01_LSWT_CoRh2O4.jl delete mode 100644 examples/01_LSWT_SU3_FeI2.jl create mode 100644 examples/02_LLD_CoRh2O4.jl delete mode 100644 examples/02_LSWT_CoRh2O4.jl delete mode 100644 examples/03_LLD_CoRh2O4.jl create mode 100644 examples/03_LSWT_SU3_FeI2.jl diff --git a/README.md b/README.md index be48c0207..7de184ad1 100644 --- a/README.md +++ b/README.md @@ -25,16 +25,15 @@ Sunny is a Julia package for modeling atomic-scale magnetism using classical spin dynamics with quantum corrections. It provides powerful tools to estimate dynamical structure factor intensities, $\mathcal{S}(𝐪,ω)$, enabling quantitative comparison with experimental scattering data, e.g., neutrons or x-rays. -A unique feature of Sunny is its treatment of spins as [SU(_N_) coherent states](https://doi.org/10.48550/arXiv.2106.14125). Each quantum spin-_S_ state is a superposition of $N=2S+1$ levels, and evolves under unitary, SU(_N_) transformations. When entanglement is neglected, the formalism allows to generalize the Landau-Lifshitz dynamics of spin dipoles to a dynamics of spin multipoles. The theory becomes especially useful for modeling materials with strong single-ion anisotropy effects (see our [FeI₂ tutorial](https://sunnysuite.github.io/Sunny.jl/dev/examples/01_LSWT_SU3_FeI2.html)). In the future, Sunny will also support explicit spin-orbit coupling, and 'units' of locally entangled spins. +A unique feature of Sunny is its treatment of spins as [SU(_N_) coherent states](https://arxiv.org/abs/2106.14125). This theory generalizes Landau-Lifshitz spin dynamics to a [nonlinear Schrödinger equation](https://arxiv.org/abs/2204.07563), which retains $N=2S+1$ levels for each quantum spin-_S_ state. The theory is effective for models with strong single-ion anisotropy (see our **[FeI₂ tutorial](https://sunnysuite.github.io/Sunny.jl/dev/examples/03_LSWT_SU3_FeI2.html)**) and for [localized "units" of strongly entangled spins](https://arxiv.org/abs/2405.16315). Efficient simulation is enabled by several algorithmic developments as listed on our [publications page](https://github.com/SunnySuite/Sunny.jl/wiki/Sunny-literature). -At low-temperatures, Sunny supports the usual linear spin wave theory for spin dipoles, and its ['multi-boson' generalization](https://doi.org/10.48550/arXiv.1307.7731). At finite temperatures, the full classical dynamics (with quantum correction factors) may be preferable to capture strongly nonlinear fluctuations. The [coupling of SU(_N_) spins to a thermal bath](https://doi.org/10.48550/arXiv.2209.01265) also makes possible the study of various non-equilibrium dynamics, e.g., thermal transport, pump-probe experiments, and spin-glass relaxation. - -Sunny provides a number of tools to facilitate the specification and solution of spin Hamiltonians. This includes spacegroup symmetry analysis, powerful Monte Carlo sampling algorithms, and interactive 3D visualization. Efficient simulation is made possible by several algorithmic developments. See our [Sunny publications](https://github.com/SunnySuite/Sunny.jl/wiki/Sunny-literature) for more information. +At low-temperatures, Sunny supports the usual linear spin wave theory for spin dipoles, and its ['multi-boson' generalization](https://arxiv.org/abs/1307.7731). At finite temperatures, Sunny can calculate the dynamical structure factor using classical spin dynamics with quantum corrections. Langevin coupling to a thermal bath additionally makes possible the study of various non-equilibrium dynamics, e.g., thermal transport, pump-probe experiments, and spin-glass relaxation. +Sunny provides a number of tools to facilitate the specification and solution of spin Hamiltonians. This includes spacegroup symmetry analysis, powerful Monte Carlo sampling algorithms, and interactive 3D visualization. ## Try it out! -New Julia users should begin with our **[Getting Started](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia)** guide. A showcase for Sunny is the **[FeI₂ tutorial](https://sunnysuite.github.io/Sunny.jl/dev/examples/01_LSWT_SU3_FeI2.html)**. This compound has a strong easy-axis anisotropy which gives rise to a single-ion bound state, and serves to exemplify the power of simulating SU(_N_) coherent states. Regarding the traditional linear spin wave theory of dipoles, see our [adaptations of the SpinW tutorials](https://sunnysuite.github.io/Sunny.jl/dev/examples/spinw/SW08_Kagome_AFM.html). +New Julia users should begin with our **[Tutorials](https://sunnysuite.github.io/Sunny.jl/dev/examples/01_LSWT_CoRh2O4)**. For traditional linear spin wave theory, see also the **[SpinW ports](https://sunnysuite.github.io/Sunny.jl/dev/examples/spinw/SW01_FM_Heseinberg_chain.html)**. Sunny is evolving rapidly. [Version History](https://sunnysuite.github.io/Sunny.jl/dev/versions/) lists the new features and breaking changes. To install a specific version of Sunny, say `v0.x`, use the command `add Sunny@0.x`. @@ -53,9 +52,7 @@ Sunny is inspired by SpinW, especially regarding symmetry analysis, model specif | Ewald summation for dipole-dipole | ❌ | ❌ | ✅ | | Programming language | C++ | Matlab | [Julia](https://julialang.org/) | -The classical SU(_N_) spin dynamics in Sunny generalizes the Landau-Lifshitz equation for $S > 1/2$ quantum spins. Linearizing and quantizing SU(_N_) dynamics yields a generalization of spin wave theory involving multi-flavor bosons. - -Codes like [Spirit](https://github.com/spirit-code/spirit) and [Vampire](https://vampire.york.ac.uk/) focus less on capturing quantum effects, but may be good options for classical dynamics of pure dipoles. +Codes like [Spirit](https://github.com/spirit-code/spirit) and [Vampire](https://vampire.york.ac.uk/) focus less on capturing quantum effects, but may be better options for the dynamics of classical dipoles, e.g., in the micromagnetics context. ## Join our community diff --git a/docs/src/structure-factor.md b/docs/src/structure-factor.md index 098f4900e..77251b8ad 100644 --- a/docs/src/structure-factor.md +++ b/docs/src/structure-factor.md @@ -338,7 +338,7 @@ sublattice $j$. Calculating the structure factor involves several steps, with various possible settings. Sunny provides tools to facilitate this calculation and to extract -information from the results. For details, please see our [tutorials](@ref "2. +information from the results. For details, please see our [tutorials](@ref "1. Spin wave simulations of CoRh₂O₄") as well as the complete [Library API](@ref). Sunny will calculate the structure factor in dimensionless, intensive units, diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl new file mode 100644 index 000000000..d7f02efc1 --- /dev/null +++ b/examples/01_LSWT_CoRh2O4.jl @@ -0,0 +1,192 @@ +# # 1. Spin wave simulations of CoRh₂O₄ +# +# This tutorial introduces Sunny through its features for performing +# conventional spin wave theory calculations. For concreteness, we consider the +# crystal CoRh₂O₄ and reproduce the calculations of [Ge et al., Phys. Rev. B 96, +# 064413](https://doi.org/10.1103/PhysRevB.96.064413). + +# ### Get Julia and Sunny +# +# Sunny is implemented in Julia, which allows for interactive development (like +# Python or Matlab) while also providing high numerical efficiency (like C++ or +# Fortran). New Julia users should begin with our **[Getting +# Started](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia)** +# guide. Sunny requires Julia 1.10 or later. +# +# From the Julia prompt, load both `Sunny` and `GLMakie`. The latter is needed +# for graphics. + +using Sunny, GLMakie + +# If these packages are not yet installed, Julia will offer to install them for +# you. If executing this script gives an error, you may need to `update` Sunny +# and GLMakie from the [built-in package +# manager](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia#the-built-in-julia-package-manager). + +# ### Units + +# The [`Units`](@ref) object defines physical constants for conversions. Select +# meV as the default energy scale and angstrom as the default length scale. + +units = Units(:meV, :angstrom); + +# ### Crystal cell +# +# A crystallographic cell may be loaded from a `.cif` file, or can be specified +# from atom positions and types. +# +# Start by defining the shape of the conventional chemical cell. CoRh₂O₄ has +# cubic spacegroup 227 (Fd-3m). Its lattice constants are 8.5 Å, and the cell +# angles are 90°. With this information, [`lattice_vectors`](@ref) constructs a +# 3×3 matrix `latvecs`. Columns of `latvecs` define the lattice vectors ``(𝐚_1, +# 𝐚_2, 𝐚_3)`` in the global Cartesian coordinate system. Conversely, columns +# of `inv(latvecs)` define the global Cartesian axes ``(\hat{x}, \hat{y}, +# \hat{z})`` in components of the lattice vectors. + +a = 8.5031 # (Å) +latvecs = lattice_vectors(a, a, a, 90, 90, 90) + +# Construct the [`Crystal`](@ref) cell from the spacegroup number 227 and one +# representative atom of each occupied Wyckoff. In the standard setting of +# spacegroup 227, position `[0, 0, 0]` belongs to Wyckoff 8a, which is the +# diamond cubic crystal. + +positions = [[0, 0, 0]] +cryst = Crystal(latvecs, positions, 227; types=["Co"], setting="1") + +# [`view_crystal`](@ref) launches an interface for interactive inspection and +# symmetry analysis. + +view_crystal(cryst) + +# ### Spin system + +# A [`System`](@ref) will define the spin model. This requires +# [`SpinInfo`](@ref) information for one representative atom per +# symmetry-distinct site. The cobalt atoms have quantum spin ``S = 3/2``. The +# ``g``-factor defines the magnetic moment ``μ = g 𝐒`` in units of the Bohr +# magneton. The option `:dipole` indicates a traditional model type, for which +# quantum spin is modeled as a dipole expectation value. + +latsize = (1, 1, 1) +S = 3/2 +sys = System(cryst, latsize, [SpinInfo(1; S, g=2)], :dipole) + +# Previous work demonstrated that inelastic neutron scattering data for CoRh₂O₄ +# is well described with a single antiferromagnetic nearest neighbor exchange, +# `J = 0.63` meV. Use [`set_exchange!`](@ref) with the bond that connects atom 1 +# to atom 3, and has zero displacement between chemical cells. Applying the +# symmetries of spacegroup 227, Sunny will propagate this interaction to the +# other nearest-neighbor bonds. Calling [`view_crystal`](@ref) with `sys` now +# shows the antiferromagnetic Heisenberg interactions as blue polkadot spheres. + +J = +0.63 # (meV) +set_exchange!(sys, J, Bond(1, 3, [0, 0, 0])) +view_crystal(sys) + +# ### Optimizing spins + +# To search for the ground state, call [`randomize_spins!`](@ref) and +# [`minimize_energy!`](@ref) in sequence. For this problem, optimization +# converges rapidly to the expected Néel order. See this with +# [`plot_spins`](@ref), where spins are colored according to their global +# ``z``-component. + +randomize_spins!(sys) +minimize_energy!(sys) +plot_spins(sys; color=[s[3] for s in sys.dipoles]) + +# The diamond lattice is bipartite, allowing each spin to perfectly anti-align +# with its 4 nearest-neighbors. Each of these 4 bonds contribute ``-JS^2`` to +# the total energy. Two sites participate in each bond, so the energy per site +# is ``-2JS^2``. Check this by calling [`energy_per_site`](@ref). + +@assert energy_per_site(sys) ≈ -2J*S^2 + +# ### Reshaping the magnetic cell + +# The same Néel order can also be described with a magnetic cell that consists +# of the 2 cobalt atoms in the primitive cell. Columns of the 3×3 `shape` matrix +# below are the primitive lattice vectors in units of the conventional, cubic +# lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. Use [`reshape_supercell`](@ref) to +# construct a system with this shape, and verify that the energy per site is +# unchanged. + +shape = [0 1 1; + 1 0 1; + 1 1 0] / 2 +sys_prim = reshape_supercell(sys, shape) +@assert energy_per_site(sys_prim) ≈ -2J*S^2 + +# Plotting the spins of `sys_prim` shows the primitive cell as a gray wireframe +# inside the conventional cubic cell. + +plot_spins(sys_prim; color=[s[3] for s in sys_prim.dipoles]) + +# ### Spin wave theory + +# With this primitive cell, we will perform a [`SpinWaveTheory`](@ref) +# calculation of the structure factor ``\mathcal{S}(𝐪,ω)``. The measurement +# [`ssf_perp`](@ref) indicates projection of the spin structure factor +# perpendicular to the direction of momentum transfer. This measurement is +# appropriate for unpolarized neutron scattering. + +swt = SpinWaveTheory(sys_prim; measure=ssf_perp(sys_prim)) + +# Define a [`q_space_path`](@ref) that connects high-symmetry points in +# reciprocal space. The ``𝐪``-points are given in reciprocal lattice units +# (RLU) for the _original_ cubic cell. For example, `[1/2, 1/2, 0]` denotes the +# sum of the first two reciprocal lattice vectors, ``𝐛_1/2 + 𝐛_2/2``. A total +# of 500 ``𝐪``-points will be sampled along the path. + +qs = [[0, 0, 0], [1/2, 0, 0], [1/2, 1/2, 0], [0, 0, 0]] +path = q_space_path(cryst, qs, 500) + +# Select [`lorentzian`](@ref) broadening with a full-width at half-maximum +# (FWHM) of 0.8 meV. The isotropic [`FormFactor`](@ref) for Co²⁺ dampens +# intensities at large ``𝐪``. + +kernel = lorentzian(fwhm=0.8) +formfactors = [FormFactor("Co2")]; + +# Calculate the single-crystal scattering [`intensities`](@ref)` along the path, +# for 300 energy points between 0 and 6 meV. Use [`plot_intensities`](@ref) to +# visualize the result. + +energies = range(0, 6, 300) +res = intensities(swt, path; energies, kernel, formfactors) +plot_intensities(res; units) + +# To directly compare with the available experimental data, perform a +# [`powder_average`](@ref) over all possible crystal orientations. Consider 200 +# ``𝐪`` magnitudes ranging from 0 to 3 inverse angstroms. Each magnitude +# defines spherical shell in reciprocal space, to be sampled with `2000` +# ``𝐪``-points. The calculation completes in just a couple seconds because the +# magnetic cell size is small. + +radii = range(0, 3, 200) # (1/Å) +res = powder_average(cryst, radii, 2000) do qs + intensities(swt, qs; energies, kernel, formfactors) +end +plot_intensities(res; units, saturation=1.0) + +# This result can be compared to experimental neutron scattering data +# from Fig. 5 of [Ge et al.](https://doi.org/10.1103/PhysRevB.96.064413) +# ```@raw html +# +# ``` + + +# ### What's next? +# +# * For more spin wave calculations of this traditional type, one can browse the +# [SpinW tutorials ported to Sunny](@ref "SW01 - FM Heisenberg chain"). +# * Spin wave theory neglects thermal fluctuations of the magnetic order. Our +# [next tutorial](@ref "2. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T*") +# demonstrates how to sample spins in thermal equilibrium, and measure +# dynamical correlations from the classical spin dynamics. +# * Sunny also offers features that go beyond the dipole approximation of a +# quantum spin via the theory of SU(_N_) coherent states. This can be +# especially useful for systems with strong single-ion anisotropy, as +# demonstrated in our [tutorial on FeI₂](@ref "3. Multi-flavor spin wave +# simulations of FeI₂"). diff --git a/examples/01_LSWT_SU3_FeI2.jl b/examples/01_LSWT_SU3_FeI2.jl deleted file mode 100644 index cdcd22afb..000000000 --- a/examples/01_LSWT_SU3_FeI2.jl +++ /dev/null @@ -1,322 +0,0 @@ -# # 1. Multi-flavor spin wave simulations of FeI₂ (Showcase) -# -# FeI₂ is an effective spin-1 material with strong single-ion anisotropy. -# Quadrupolar fluctuations give rise to a single-ion bound state that cannot be -# described by a dipole-only model. This tutorial illustrates how to use the -# linear spin wave theory of SU(3) coherent states (i.e. 2-flavor bosons) to -# model the magnetic behavior in FeI₂. The original study was performed in [Bai -# et al., Nature Physics 17, 467–472 -# (2021)](https://doi.org/10.1038/s41567-020-01110-1). -# -# ```@raw html -# -# ``` -# -# The Fe atoms are arranged in stacked triangular layers. The effective spin -# Hamiltonian takes the form, -# -# ```math -# \mathcal{H}=\sum_{(i,j)} 𝐒_i ⋅ J_{ij} 𝐒_j - D\sum_i \left(S_i^z\right)^2, -# ``` -# -# where the set of exchange matrices ``J_{ij}`` between bonded sites ``(i,j)`` -# includes competing ferromagnetic and antiferromagnetic interactions. This -# model also includes a strong easy axis anisotropy, ``D > 0``. -# -# We will formulate this Hamiltonian in Sunny and then calculate its dynamic -# structure factor. - -# ## Get Julia and Sunny -# -# Sunny is implemented in Julia. This is a relatively new programming language -# that allows for interactive development (like Python or Matlab) while also -# providing high numerical efficiency (like C++ or Fortran). New Julia users may -# wish to take a look at our [Getting Started with -# Julia](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia) -# guide. Sunny requires Julia 1.9 or later. -# -# From the Julia prompt, load `Sunny` and `GLMakie` for 3D graphics. - -using Sunny, GLMakie - -# If these packages are not yet installed, Julia should offer to install them -# using its built-in package management system. If old versions are installed, -# you may need to update them to run this tutorial. - -# ## Crystals -# -# A [`Crystal`](@ref) describes the crystallographic unit cell and will usually -# be loaded from a `.cif` file. Here, we instead build a crystal by listing all -# atoms and their types. - -units = Units(:meV, :angstrom) -a = b = 4.05012 # Lattice constants for triangular lattice -c = 6.75214 # Spacing in the z-direction - -latvecs = lattice_vectors(a, b, c, 90, 90, 120) # A 3x3 matrix of lattice vectors that - ## define the conventional unit cell -positions = [[0, 0, 0], [1/3, 2/3, 1/4], [2/3, 1/3, 3/4]] # Positions of atoms in fractions - ## of lattice vectors -types = ["Fe", "I", "I"] -FeI2 = Crystal(latvecs, positions; types) - -# Observe that Sunny inferred the space group, 'P -3 m 1' (164) and labeled the -# atoms according to their point group symmetries. - -# Only the Fe atoms are magnetic, so we discard the I ions using -# [`subcrystal`](@ref). - -cryst = subcrystal(FeI2, "Fe") - -# Importantly, `cryst` retains the spacegroup symmetry of the full FeI₂ crystal. -# This information will be used, for example, to propagate exchange interactions -# between symmetry-equivalent bonds. -# -# In a running Julia environment, the crystal can be viewed interactively using -# [`view_crystal`](@ref). - -view_crystal(cryst) - -# ## Symmetry analysis -# -# The command [`print_symmetry_table`](@ref) provides a list of all the -# symmetry-allowed interactions up to a cutoff distance. - -print_symmetry_table(cryst, 8.0) - -# The allowed ``g``-tensor is expressed as a 3×3 matrix in the free coefficients -# `A`, `B`, ... The allowed single-ion anisotropy is expressed as a linear -# combination of Stevens operators. The latter correspond to polynomials of the -# spin operators, as we will describe below. -# -# The allowed exchange interactions are given as a 3×3 matrix for representative -# bonds. The notation `Bond(i, j, n)` indicates a bond between atom indices `i` -# and `j`, with cell offset `n`. In the general case, it will be necessary to -# associate atom indices with their positions in the unit cell; these can be -# viewed with `display(cryst)`. Note that the order of the pair ``(i, j)`` is -# significant if the exchange tensor contains antisymmetric -# Dzyaloshinskii–Moriya (DM) interactions. -# -# In the case of FeI₂, `Bond(1, 1, [1,0,0])` is one of the 6 nearest-neighbor -# Fe-Fe bonds on a triangular lattice layer, and `Bond(1, 1, [0,0,1])` is an -# Fe-Fe bond between layers. - -# ## Building a spin System - -# In constructing a spin [`System`](@ref), we must provide several additional -# details about the spins. - -sys = System(cryst, (4, 4, 4), [SpinInfo(1, S=1, g=2)], :SUN, seed=2) - -# This system includes ``4×4×4`` unit cells, i.e. 64 Fe atoms, each with spin -# ``S=1`` and a ``g``-factor of 2. Quantum mechanically, spin ``S=1`` involves a -# superposition of ``2S+1=3`` distinct angular momentum states. In `:SUN` mode, -# this superposition will be modeled explicitly using the formalism of SU(3) -# coherent states, which captures both dipolar and quadrupolar fluctuations. For -# the more traditional dipole dynamics, use `:dipole` mode instead. - -# Next we will use [`set_exchange!`](@ref) to assign interaction to bonds. Sunny -# will automatically propagate each interaction to all symmetry-equivalent bonds -# in the unit cell. The FeI₂ interactions below follow [Bai et -# al](https://doi.org/10.1038/s41567-020-01110-1). - -J1pm = -0.236 -J1pmpm = -0.161 -J1zpm = -0.261 -J2pm = 0.026 -J3pm = 0.166 -J′0pm = 0.037 -J′1pm = 0.013 -J′2apm = 0.068 - -J1zz = -0.236 -J2zz = 0.113 -J3zz = 0.211 -J′0zz = -0.036 -J′1zz = 0.051 -J′2azz = 0.073 - -J1xx = J1pm + J1pmpm -J1yy = J1pm - J1pmpm -J1yz = J1zpm - -set_exchange!(sys, [J1xx 0.0 0.0; - 0.0 J1yy J1yz; - 0.0 J1yz J1zz], Bond(1,1,[1,0,0])) -set_exchange!(sys, [J2pm 0.0 0.0; - 0.0 J2pm 0.0; - 0.0 0.0 J2zz], Bond(1,1,[1,2,0])) -set_exchange!(sys, [J3pm 0.0 0.0; - 0.0 J3pm 0.0; - 0.0 0.0 J3zz], Bond(1,1,[2,0,0])) -set_exchange!(sys, [J′0pm 0.0 0.0; - 0.0 J′0pm 0.0; - 0.0 0.0 J′0zz], Bond(1,1,[0,0,1])) -set_exchange!(sys, [J′1pm 0.0 0.0; - 0.0 J′1pm 0.0; - 0.0 0.0 J′1zz], Bond(1,1,[1,0,1])) -set_exchange!(sys, [J′2apm 0.0 0.0; - 0.0 J′2apm 0.0; - 0.0 0.0 J′2azz], Bond(1,1,[1,2,1])) - -# The function [`set_onsite_coupling!`](@ref) assigns a single-ion anisotropy. -# The argument can be constructed using [`spin_matrices`](@ref) or -# [`stevens_matrices`](@ref). Here we use Julia's anonymous function syntax to -# assign an easy-axis anisotropy along the direction ``\hat{z}``. - -D = 2.165 -set_onsite_coupling!(sys, S -> -D*S[3]^2, 1) - -# # Calculating structure factor intensities - -# In the remainder of this tutorial, we will examine Sunny's tools for -# calculating the dynamical structure factor using a [multi-boson -# generalization](https://arxiv.org/abs/1307.7731) of linear spin wave theory -# (LSWT). This theory describes non-interacting quasi-particle excitations that -# hybridize dipolar and quadrupolar modes. - -# ## 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) -minimize_energy!(sys) - -# A positive number above indicates that the procedure has converged to a local -# energy minimum. The configuration, however, may still have defects. This can -# be checked by visualizing the spins, colored according to their -# ``z``-components. - -plot_spins(sys; color=[s[3] for s in sys.dipoles]) - -# A different understanding of the magnetic ordering can be obtained by moving -# to Fourier space. The 'instantaneous' structure factor ``𝒮(𝐪)`` is an -# experimental observable. To investigate ``𝒮(𝐪)`` as true 3D data, Sunny -# provides [`SampledCorrelationsStatic`](@ref). Here, however, we will use -# [`print_wrapped_intensities`](@ref), which gives average intensities for the -# individual Bravais sublattices (in effect, all wavevectors are wrapped to the -# first Brillouin zone). - -print_wrapped_intensities(sys) - -# The result will likely be approximately consistent with the known zero-field -# energy-minimizing magnetic structure of FeI₂, which is single-``Q`` (two-up, -# two-down antiferromagnetic order). 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 nature, 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₂, it could be -# entered by hand using [`set_dipole!`](@ref). Alternatively, in the case of -# FeI₂, 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₂ 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 suggests a corresponding magnetic -# supercell in units of the crystal lattice vectors. - -suggest_magnetic_supercell([[0, -1/4, 1/4]]) - -# The system returned by [`reshape_supercell`](@ref) is smaller, and is sheared -# relative to the original system. This makes it much easier to find the global -# energy minimum. - -sys_min = reshape_supercell(sys, [1 0 0; 0 2 1; 0 -2 1]) -randomize_spins!(sys_min) -minimize_energy!(sys_min); - -# Plot the system again, now including "ghost" spins out to 12Å - -plot_spins(sys_min; color=[s[3] for s in sys_min.dipoles], ghost_radius=12) - -# ## Linear spin wave theory -# -# Now that we have found the ground state for a magnetic supercell, we can -# perform zero-temperature calculations using linear spin wave theory. - -# The function [`q_space_path`](@ref) will linearly sample a path between the -# provided ``q``-points in reciprocal lattice units (RLU). Here, we use a total -# of 500 wavevectors. - -qs = [[0,0,0], [1,0,0], [0,1,0], [1/2,0,0], [0,1,0], [0,0,0]] -path = q_space_path(cryst, qs, 500) - -# Construct a [`SpinWaveTheory`](@ref) object for the magnetic supercell and -# calculate scattering intensities with [`intensities_bands`](@ref). The -# measurement [`ssf_perp`](@ref) will project the dynamical spin structure -# factor onto the space perpendicular to the momentum transfer ``𝐪``, which is -# appropriate for an unpolarized neutron beam. - -swt = SpinWaveTheory(sys_min; measure=ssf_perp(sys_min)) -res = intensities_bands(swt, path) -plot_intensities(res; units) - -# To make comparisons with inelastic neutron scattering (INS) data, one can -# employ empirical broadening. Select [`lorentzian`](@ref) broadening, with a -# full-width at half-maximum of 0.3 meV. We will calculate intensities for 300 -# discrete energies between 0 and 10 meV. - -kernel = lorentzian(fwhm=0.3) -energies = range(0, 10, 300); # 0 < ω < 10 (meV) - -# A real FeI₂ sample will exhibit spontaneous breaking of its 3-fold rotational -# symmetry about the ``ẑ``-axis. We use [`domain_average`](@ref) to effectively -# average the broadened [`intensities`](@ref) calculations over the three -# possible domain orientations. In practice, this involves rotating the -# ``𝐪``-points by 0°, 120°, and 240° angles. - -rotations = [([0,0,1], n*(2π/3)) for n in 0:2] -weights = [1, 1, 1] -res = domain_average(cryst, path; rotations, weights) do path_rotated - intensities(swt, path_rotated; energies, kernel) -end -plot_intensities(res; units, colormap=:viridis) - -# This result can be directly compared to experimental neutron scattering data -# from [Bai et al.](https://doi.org/10.1038/s41567-020-01110-1) -# ```@raw html -# -# ``` -# -# (The publication figure accidentally used a non-standard coordinate system to -# label the wave vectors.) -# -# To get this agreement, the use of SU(3) coherent states is essential. In other -# words, we needed a theory of multi-flavored bosons. The lower band has large -# quadrupolar character, and arises from the strong easy-axis anisotropy of -# FeI₂. By setting `mode = :SUN`, the calculation captures this coupled -# dipole-quadrupole dynamics. -# -# An interesting exercise is to repeat the same study, but using `mode = -# :dipole` instead of `:SUN`. That alternative choice would constrain the -# coherent state dynamics to the space of dipoles only. -# -# ## What's next? -# -# The multi-boson linear spin wave theory, applied above, can be understood as -# the quantization of a certain generalization of the Landau-Lifshitz spin -# dynamics. Rather than dipoles, this dynamics takes places on the space of -# [SU(_N_) coherent states](https://arxiv.org/abs/2106.14125). -# -# The full SU(_N_) coherent state dynamics, with appropriate quantum correction -# factors, can be useful to model finite temperature scattering data. In -# particular, it captures certain anharmonic effects due to thermal -# fluctuations. See our [generalized spin dynamics tutorial](@ref "4. -# Generalized spin dynamics of FeI₂ at finite *T*"). -# -# The classical dynamics is also a good starting point to study non-equilibrium -# phenomena. Empirical noise and damping terms can be used to model [coupling to -# a thermal bath](https://arxiv.org/abs/2209.01265). This yields a Langevin -# dynamics of SU(_N_) coherent states. Our [dynamical SU(_N_) quench](@ref "6. -# Dynamical quench into CP² skyrmion liquid") tutorial illustrates how a -# temperature quench can give rise to novel liquid phase of CP² skyrmions. diff --git a/examples/02_LLD_CoRh2O4.jl b/examples/02_LLD_CoRh2O4.jl new file mode 100644 index 000000000..ae043f432 --- /dev/null +++ b/examples/02_LLD_CoRh2O4.jl @@ -0,0 +1,179 @@ +# # 2. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T* +# +# In the [previous tutorial](@ref "1. Spin wave simulations of CoRh₂O₄"), we +# used spin wave theory to calculate the dynamical spin structure factor of +# CoRh₂O₄. Here, we perform a similar calculation using equilibrium samples from +# the Boltzmann distribution at finite ``T``. For each sampled spin +# configuration, we will simulate the classical Landau-Lifshitz spin dynamics +# and extract dynamical spin-spin correlations. After applying a +# classical-to-quantum correction factor, the resulting intensities can be +# compared to inelastic neutron scattering data. + +# Construct the system as in the [previous tutorial](@ref "1. Spin wave +# simulations of CoRh₂O₄"). For this antiferromagnetic model on the diamond +# cubic lattice, the ground state is unfrustrated Néel order. + +using Sunny, GLMakie + +units = Units(:meV, :angstrom) +a = 8.5031 # (Å) +latvecs = lattice_vectors(a, a, a, 90, 90, 90) +cryst = Crystal(latvecs, [[0,0,0]], 227, setting="1") + +sys = System(cryst, (1,1,1), [SpinInfo(1; S=3/2, g=2)], :dipole) +J = 0.63 # (meV) +set_exchange!(sys, J, Bond(1, 3, [0,0,0])) +randomize_spins!(sys) +minimize_energy!(sys) +plot_spins(sys; color=[s[3] for s in sys.dipoles]) + +# Use [`resize_supercell`](@ref) to build a new system with a lattice of +# 10×10×10 chemical unit cells. The ground state Néel order is retained. +# Increasing the system size further would reduce finite-size artifacts, and +# would increase momentum-space resolution, but would make the simulations +# slower. + +sys = resize_supercell(sys, (10, 10, 10)) +plot_spins(sys; color=[s[3] for s in sys.dipoles]) + +# ### Langevin dynamics for sampling + +# We will be using a [`Langevin`](@ref) spin dynamics to thermalize the system. +# This dynamics is a variant of the Landau-Lifshitz equation that incorporates +# noise and dissipation terms, which are linked by a fluctuation-dissipation +# theorem. The temperature 6 K ≈ 1.38 meV is slightly above ordering for this +# model. The dimensionless `damping` magnitude sets a timescale for coupling to +# the implicit thermal bath; 0.2 is usually a good choice. + +langevin = Langevin(; damping=0.2, kT=16*units.K) + +# Use [`suggest_timestep`](@ref) to select an integration timestep. A +# dimensionless error tolerance of `1e-2` is usually a good choice. The +# suggested timestep will vary according to the magnetic configuration. It is +# reasonable to start from an energy-minimized configuration. + +suggest_timestep(sys, langevin; tol=1e-2) +langevin.dt = 0.025; + +# Now run a Langevin trajectory to sample spin configurations. Keep track of the +# energy per site at each time step. + +energies = [energy_per_site(sys)] +for _ in 1:1000 + step!(sys, langevin) + push!(energies, energy_per_site(sys)) +end + +# From the relaxed spin configuration, we can learn that `dt` was a little +# smaller than necessary; increasing it will make the remaining simulations +# faster. + +suggest_timestep(sys, langevin; tol=1e-2) +langevin.dt = 0.042; + +# Plot energy versus time using the [Makie `lines` +# function](https://docs.makie.org/stable/reference/plots/lines). The plateau +# suggests that the system has reached thermal equilibrium. + +lines(energies, color=:blue, figure=(size=(600,300),), axis=(xlabel="Timesteps", ylabel="Energy (meV)")) + +# Plot the spins colored by their alignment with a reference spin at the origin. +# The field `sys.dipoles` is a 4D array storing the spin dipole data. The first +# three indices of label the chemical cell, while the fourth index labels an +# atom within the cell. Note that Julia arrays use 1-based indexing. Thermal +# fluctuations are apparent in the plot. + +s_ref = sys.dipoles[1,1,1,1] +plot_spins(sys; color=[s'*s_ref for s in sys.dipoles]) + +# ### Static structure factor + +# Use [`SampledCorrelationsStatic`](@ref) to estimate spatial correlations for +# configurations in classical thermal equilibrium. Each call to +# [`add_sample!`](@ref) will accumulate data for the current spin snapshot. + +sc = SampledCorrelationsStatic(sys; measure=ssf_perp(sys)) +add_sample!(sc, sys) # Accumulate the newly sampled structure factor into `sf` + +# Collect 20 additional samples. Perform 100 Langevin time-steps between +# measurements to approximately decorrelate the sample in thermal equilibrium. + +for _ in 1:20 + for _ in 1:100 + step!(sys, langevin) + end + add_sample!(sc, sys) +end + +# Use [`q_space_grid`](@ref) to define a slice of momentum space ``[H, K, 0]``, +# where ``H`` and ``K`` each range from -10 to 10 in RLU. This command produces +# a 200×200 grid of sample points. + +grid = q_space_grid(cryst, [1, 0, 0], range(-10, 10, 200), [0, 1, 0], (-10, 10)) + +# Calculate and plot the instantaneous structure factor on the slice by +# integrating over all energy values ω. We employ the appropriate +# [`FormFactor`](@ref) for Co2⁺. Selecting `saturation = 1.0` sets the color +# saturation point to the maximum intensity value. This is reasonable because we +# are above the ordering temperature, and do not have sharp Bragg peaks. + +formfactors = [FormFactor("Co2")] +res = intensities_static(sc, grid; formfactors) +plot_intensities(res; saturation=1.0) + +# ### Dynamical structure factor + +# To collect statistics for the _dynamical_ structure factor intensities +# ``\mathcal{S}(𝐪,ω)`` at finite temperature, use +# [`SampledCorrelations`](@ref). It requires a range of `energies` to resolve, +# which will be associated with frequencies of the classical spin dynamics. The +# integration timestep `dt` can be somewhat larger than that used by the +# Langevin dynamics. + +dt = 2*langevin.dt +energies = range(0, 6, 50) +sc = SampledCorrelations(sys; dt, energies, measure=ssf_perp(sys)) + +# Like before, use Langevin dynamics to sample spin configurations from thermal +# equilibrium. Now, however, each call to [`add_sample!`](@ref) will run a +# classical spin dynamics trajectory and measure dynamical correlations. To make +# the tutorial run quickly, we average over just 5 trajectories. To make a +# publication quality figure, this number should be significantly increased for +# better statistics. + +for _ in 1:5 + for _ in 1:100 + step!(sys, langevin) + end + add_sample!(sc, sys) +end + +# Select points that define a piecewise-linear path through reciprocal space, +# and a sampling density. + +qs = [[3/4, 3/4, 0], + [ 0, 0, 0], + [ 0, 1/2, 1/2], + [1/2, 1, 0], + [ 0, 1, 0], + [1/4, 1, 1/4], + [ 0, 1, 0], + [ 0, -4, 0]] +path = q_space_path(cryst, qs, 500) + +# Calculate and plot the intensities along this path. + +res = intensities(sc, path; energies, langevin.kT) +plot_intensities(res; units) + +# ### Powder averaged intensity + +# Define spherical shells in reciprocal space via their radii, in absolute units +# of 1/Å. For each shell, calculate and average the intensities at 350 +# ``𝐪``-points + +radii = range(0, 3.5, 200) # (1/Å) +res = powder_average(cryst, radii, 350) do qs + intensities(sc, qs; energies, formfactors, langevin.kT) +end +plot_intensities(res; units) diff --git a/examples/02_LSWT_CoRh2O4.jl b/examples/02_LSWT_CoRh2O4.jl deleted file mode 100644 index 6520a8269..000000000 --- a/examples/02_LSWT_CoRh2O4.jl +++ /dev/null @@ -1,106 +0,0 @@ -# # 2. Spin wave simulations of CoRh₂O₄ -# -# This tutorial illustrates the conventional spin wave theory of dipoles. We -# consider a simple model of the diamond-cubic crystal CoRh₂O₄, with parameters -# extracted from [Ge et al., Phys. Rev. B 96, -# 064413](https://doi.org/10.1103/PhysRevB.96.064413). - -using Sunny, GLMakie - -# Construct a diamond [`Crystal`](@ref) in the conventional (non-primitive) -# cubic unit cell. Sunny will populate all eight symmetry-equivalent sites when -# given the international spacegroup number 227 ("Fd-3m") and the appropriate -# setting. For this spacegroup, there are two conventional translations of the -# unit cell, and it is necessary to disambiguate through the `setting` keyword -# argument. (On your own: what happens if `setting` is omitted?) - -units = Units(:meV, :angstrom) -a = 8.5031 # (Å) -latvecs = lattice_vectors(a, a, a, 90, 90, 90) -cryst = Crystal(latvecs, [[0,0,0]], 227, setting="1") - -# In a running Julia environment, the crystal can be viewed interactively using -# [`view_crystal`](@ref). - -view_crystal(cryst) - -# Construct a [`System`](@ref) with quantum spin ``S=3/2`` constrained to the -# space of dipoles. Including an antiferromagnetic nearest neighbor interaction -# `J` will favor Néel order. To optimize this magnetic structure, it is -# sufficient to employ a magnetic lattice consisting of a single crystal unit -# cell, `latsize=(1,1,1)`. Passing an explicit random number `seed` will ensure -# repeatable results. - -latsize = (1, 1, 1) -S = 3/2 -J = 0.63 # (meV) -sys = System(cryst, latsize, [SpinInfo(1; S, g=2)], :dipole; seed=0) -set_exchange!(sys, J, Bond(1, 3, [0,0,0])) - -# In the ground state, each spin is exactly anti-aligned with its 4 -# nearest-neighbors. Because every bond contributes an energy of ``-JS^2``, the -# energy per site is ``-2JS^2``. In this calculation, a factor of 1/2 avoids -# double-counting the bonds. Due to lack of frustration, direct energy -# minimization is successful in finding the ground state. - -randomize_spins!(sys) -minimize_energy!(sys) - -@assert energy_per_site(sys) ≈ -2J*S^2 - -# Plotting the spins confirms the expected Néel order. Note that the overall, -# global rotation of dipoles is arbitrary. - -s0 = sys.dipoles[1,1,1,1] -plot_spins(sys; color=[s'*s0 for s in sys.dipoles]) - -# For numerical efficiency, it is helpful to work with the smallest possible -# magnetic supercell; in this case, it is the primitive cell. The columns of the -# 3×3 `shape` matrix define the lattice vectors of the primitive cell as -# multiples of the conventional, cubic lattice vectors. After transforming the -# system with [`reshape_supercell`](@ref), the energy per site remains the same. -shape = [0 1 1; - 1 0 1; - 1 1 0] / 2 -sys_prim = reshape_supercell(sys, shape) -@assert energy_per_site(sys_prim) ≈ -2J*S^2 -plot_spins(sys_prim; color=[s'*s0 for s in sys_prim.dipoles]) - -# Now estimate ``𝒮(𝐪,ω)`` with [`SpinWaveTheory`](@ref). - -swt = SpinWaveTheory(sys_prim; measure=ssf_perp(sys_prim)) - -# For the "single crystal" result, we use [`q_space_path`](@ref) to construct a -# 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, 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 -# intensities. The isotropic [`FormFactor`](@ref) for Cobalt(2+) dampens -# intensities at large ``𝐪``. - -kernel = lorentzian(fwhm=0.8) -formfactors = [FormFactor("Co2")] -energies = range(0, 6, 300) -res = intensities(swt, path; energies, kernel, formfactors) -plot_intensities(res; units) - -# We use [`powder_average`](@ref) to average intensities over all possible -# crystal orientations. Perform this calculation for 200 momentum magnitudes, -# ranging from 0 to 3 inverse angstroms. Each ``𝐪``-magnitude defines a -# spherical shell in reciprocal space. Sample it with `2000` wavevectors of -# quasi-uniform distribution. - -radii = range(0, 3, 200) # (1/Å) -res = powder_average(cryst, radii, 2000) do qs - intensities(swt, qs; energies, kernel, formfactors) -end -plot_intensities(res; units, saturation=1.0) - -# This result can be compared to experimental neutron scattering data -# from Fig. 5 of [Ge et al.](https://doi.org/10.1103/PhysRevB.96.064413) -# ```@raw html -# -# ``` diff --git a/examples/03_LLD_CoRh2O4.jl b/examples/03_LLD_CoRh2O4.jl deleted file mode 100644 index 098442b95..000000000 --- a/examples/03_LLD_CoRh2O4.jl +++ /dev/null @@ -1,158 +0,0 @@ -# # 3. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T* - -using Sunny, GLMakie - -# ### System construction - -# Construct the system as in the previous [CoRh₂O₄ tutorial](@ref "2. Spin wave -# simulations of CoRh₂O₄"). After optimization, the system will be in an -# unfrustrated antiferromagnetic ground state. - -a = 8.5031 # (Å) -latvecs = lattice_vectors(a, a, a, 90, 90, 90) -cryst = Crystal(latvecs, [[0,0,0]], 227, setting="1") - -units = Units(:meV, :angstrom) -latsize = (2, 2, 2) -S = 3/2 -J = 0.63 # (meV) -sys = System(cryst, latsize, [SpinInfo(1; S, g=2)], :dipole; seed=0) -set_exchange!(sys, J, Bond(1, 3, [0,0,0])) -randomize_spins!(sys) -minimize_energy!(sys) - -# Use [`resize_supercell`](@ref) to build a new system with a lattice of -# 10×10×10 unit cells. The desired Néel order is retained. - -sys = resize_supercell(sys, (10, 10, 10)) -@assert energy_per_site(sys) ≈ -2J*S^2 - -# We will be using a [`Langevin`](@ref) spin dynamics to thermalize the system. -# This dynamics involves a dimensionless `damping` magnitude and target -# temperature `kT` for thermal fluctuations. - -kT = 16 * units.K # 16 K ≈ 1.38 meV, slightly below ordering temperature -langevin = Langevin(; damping=0.2, kT) - -# Use [`suggest_timestep`](@ref) to select an integration timestep for the given -# error tolerance, e.g. `tol=1e-2`. The spin configuration in `sys` should -# ideally be relaxed into thermal equilibrium, but the current, energy-minimized -# configuration will also work reasonably well. - -suggest_timestep(sys, langevin; tol=1e-2) -langevin.dt = 0.025; - -# Now run a dynamical trajectory to sample spin configurations. Keep track of -# the energy per site at each time step. - -energies = [energy_per_site(sys)] -for _ in 1:1000 - step!(sys, langevin) - push!(energies, energy_per_site(sys)) -end - -# Now that the spin configuration has relaxed, we can learn that `dt` was a -# little smaller than necessary; increasing it will make the remaining -# simulations faster. - -suggest_timestep(sys, langevin; tol=1e-2) -langevin.dt = 0.042; - -# The energy per site has converged, which suggests that the system has reached -# thermal equilibrium. - -plot(energies, color=:blue, figure=(size=(600,300),), axis=(xlabel="Timesteps", ylabel="Energy (meV)")) - -# Thermal fluctuations are apparent in the spin configuration. - -S_ref = sys.dipoles[1,1,1,1] -plot_spins(sys; color=[s'*S_ref for s in sys.dipoles]) - -# ### Static structure factor for the classical distribution - -# Use [`SampledCorrelationsStatic`](@ref) to estimate spatial correlations for -# configurations in classical thermal equilibrium. Each call to -# [`add_sample!`](@ref) will accumulate data for the current spin snapshot. - -sc = SampledCorrelationsStatic(sys; measure=ssf_perp(sys)) -add_sample!(sc, sys) # Accumulate the newly sampled structure factor into `sf` - -# Collect 20 additional samples. About 100 Langevin time-steps between -# measurements is sufficient to approximately decorrelate each sample for the -# thermal equilibrium. - -for _ in 1:20 - for _ in 1:100 - step!(sys, langevin) - end - add_sample!(sc, sys) -end - -# Use [`q_space_grid`](@ref) to define a slice of momentum space ``[H, K, 0]``, -# where ``H`` and ``K`` each range from -10 to 10 in RLU. This command produces -# a 200×200 grid of sample points. - -grid = q_space_grid(cryst, [1, 0, 0], range(-10, 10, 200), [0, 1, 0], (-10, 10)) - -# Calculate and plot the instantaneous structure factor on the slice by -# integrating over all energy values ω. We employ the appropriate -# [`FormFactor`](@ref) for Co2⁺. - -formfactors = [FormFactor("Co2")] -res = intensities_static(sc, grid; formfactors) -plot_intensities(res; saturation=0.995) - - -# ### Dynamical structure factor - -# To collect statistics for the _dynamical_ structure factor intensities -# ``I(𝐪,ω)`` at finite temperature, use instead [`SampledCorrelations`](@ref). -# It requires a range of `energies` to resolve, which will be associated with -# frequencies of the classical spin dynamics. The integration timestep `dt` can -# be somewhat larger than that used by the Langevin dynamics. - -dt = 2*langevin.dt -energies = range(0, 6, 50) -sc = SampledCorrelations(sys; dt, energies, measure=ssf_perp(sys)) - -# Again use Langevin dynamics to sample spin configurations from thermal -# equilibrium. For each sample, use [`add_sample!`](@ref) to run a classical -# spin dynamics trajectory and measure dynamical correlations. Here we average -# over just 5 samples, but this number could be increased for better statistics. - -for _ in 1:5 - for _ in 1:100 - step!(sys, langevin) - end - add_sample!(sc, sys) -end - -# Select points that define a piecewise-linear path through reciprocal space, -# and a sampling density. - -qs = [[3/4, 3/4, 0], - [ 0, 0, 0], - [ 0, 1/2, 1/2], - [1/2, 1, 0], - [ 0, 1, 0], - [1/4, 1, 1/4], - [ 0, 1, 0], - [ 0, -4, 0]] -qpts = q_space_path(cryst, qs, 500) - -# Calculate ``I(𝐪, ω)`` intensities along this path and plot. - -res = intensities(sc, qpts; energies, formfactors, kT) -plot_intensities(res; units, saturation=0.85, colormap=:viridis) - -# ### Powder averaged intensity - -# Define spherical shells in reciprocal space via their radii, in absolute units -# of 1/Å. For each shell, calculate and average the intensities at 100 -# ``𝐪``-points, sampled approximately uniformly. - -radii = range(0, 3.5, 200) # (1/Å) -res = powder_average(cryst, radii, 350) do qs - intensities(sc, qs; energies, formfactors, kT) -end -plot_intensities(res; units, colormap=:viridis) diff --git a/examples/03_LSWT_SU3_FeI2.jl b/examples/03_LSWT_SU3_FeI2.jl new file mode 100644 index 000000000..528dbff97 --- /dev/null +++ b/examples/03_LSWT_SU3_FeI2.jl @@ -0,0 +1,253 @@ +# # 3. Multi-flavor spin wave simulations of FeI₂ +# +# This tutorial illustrates various powerful features in Sunny, including +# symmetry analysis, energy minimization, and spin wave theory with multi-flavor +# bosons. +# +# Our context will be FeI₂, an effective spin-1 material with strong single-ion +# anisotropy. Quadrupolar fluctuations give rise to a single-ion bound state +# that is observable in neutron scattering, and cannot be described by a +# dipole-only model. This tutorial illustrates how to use the linear spin wave +# theory of SU(3) coherent states (i.e. 2-flavor bosons) to model the magnetic +# spectrum of FeI₂. The original study was performed in [Bai et al., Nature +# Physics 17, 467–472 (2021)](https://doi.org/10.1038/s41567-020-01110-1). +# +# ```@raw html +# +# ``` +# +# The Fe atoms are arranged in stacked triangular layers. The effective spin +# Hamiltonian takes the form, +# +# ```math +# \mathcal{H}=\sum_{(i,j)} 𝐒_i ⋅ J_{ij} 𝐒_j - D\sum_i \left(S_i^z\right)^2, +# ``` +# +# where the exchange matrices ``J_{ij}`` between bonded sites ``(i,j)`` include +# competing ferromagnetic and antiferromagnetic interactions. This model also +# includes a strong easy axis anisotropy, ``D > 0``. + +# Load packages. + +using Sunny, GLMakie + +# Construct the chemical cell of FeI₂ by specifying the lattice vectors and the +# full set of atoms. + +units = Units(:meV, :angstrom) +a = b = 4.05012 # Lattice constants for triangular lattice (Å) +c = 6.75214 # Spacing between layers (Å) +latvecs = lattice_vectors(a, b, c, 90, 90, 120) +positions = [[0, 0, 0], [1/3, 2/3, 1/4], [2/3, 1/3, 3/4]] +types = ["Fe", "I", "I"] +cryst = Crystal(latvecs, positions; types) + +# Observe that Sunny inferred the space group, 'P -3 m 1' (164) and labeled the +# atoms according to their point group symmetries. Only the Fe atoms are +# magnetic, so we focus on them with [`subcrystal`](@ref). Importantly, this +# operation preserves the spacegroup symmetries. + +cryst = subcrystal(cryst, "Fe") +view_crystal(cryst) + +# ### Symmetry analysis +# +# The command [`print_symmetry_table`](@ref) provides a list of all the +# symmetry-allowed interactions out to 8 Å. + +print_symmetry_table(cryst, 8.0) + +# The allowed ``g``-tensor is expressed as a 3×3 matrix in the free coefficients +# `A`, `B`, ... The allowed single-ion anisotropy is expressed as a linear +# combination of Stevens operators. The latter correspond to polynomials of the +# spin operators, as we will describe below. +# +# The allowed exchange interactions are given as 3×3 matrices for representative +# bonds. The notation `Bond(i, j, n)` indicates a bond between atom indices `i` +# and `j`, with cell offset `n`. Note that the order of the pair ``(i, j)`` is +# significant if the exchange tensor contains antisymmetric +# Dzyaloshinskii–Moriya (DM) interactions. +# +# The bonds can be visualized in the `view_crystal` interface. By default, +# `Bond(1, 1, [1,0,0])` is toggled on, to show the 6 nearest-neighbor Fe-Fe +# bonds on a triangular lattice layer. Toggling `Bond(1, 1, [0,0,1])` shows the +# Fe-Fe bond between layers, etc. + +# ### Defining the spin model + +# Construct a system [`System`](@ref) with spin-1 and g=2 for the Fe ions. +# +# Recall that a quantum spin-1 state is, in general, a linear combination of +# basis states ``|m⟩`` with well-defined angular momentum ``m = -1, 0, 1``. FeI₂ +# has a strong easy-axis anisotropy, which stabilizes a single-ion bound state +# of zero angular momentum, ``|m=0⟩``. Such a bound state is inaccessible to +# traditional spin wave theory, which works with dipole expectation values of +# fixed magnitude. This physics can, however, be captured with a theory of +# SU(_N_) coherent states, where ``N = 2S+1 = 3`` is the number of levels. We +# will therefore select `:SUN` mode instead of `:dipole` mode. +# +# Selecting an optional random number `seed` will make the calculations exactly +# reproducible. + +sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=1, g=2)], :SUN, seed=2) + +# Set the exchange interactions for FeI₂ following the fits of [Bai et +# al.](https://doi.org/10.1038/s41567-020-01110-1) + +J1pm = -0.236 # (meV) +J1pmpm = -0.161 +J1zpm = -0.261 +J2pm = 0.026 +J3pm = 0.166 +J′0pm = 0.037 +J′1pm = 0.013 +J′2apm = 0.068 + +J1zz = -0.236 +J2zz = 0.113 +J3zz = 0.211 +J′0zz = -0.036 +J′1zz = 0.051 +J′2azz = 0.073 + +J1xx = J1pm + J1pmpm +J1yy = J1pm - J1pmpm +J1yz = J1zpm + +set_exchange!(sys, [J1xx 0.0 0.0; + 0.0 J1yy J1yz; + 0.0 J1yz J1zz], Bond(1,1,[1,0,0])) +set_exchange!(sys, [J2pm 0.0 0.0; + 0.0 J2pm 0.0; + 0.0 0.0 J2zz], Bond(1,1,[1,2,0])) +set_exchange!(sys, [J3pm 0.0 0.0; + 0.0 J3pm 0.0; + 0.0 0.0 J3zz], Bond(1,1,[2,0,0])) +set_exchange!(sys, [J′0pm 0.0 0.0; + 0.0 J′0pm 0.0; + 0.0 0.0 J′0zz], Bond(1,1,[0,0,1])) +set_exchange!(sys, [J′1pm 0.0 0.0; + 0.0 J′1pm 0.0; + 0.0 0.0 J′1zz], Bond(1,1,[1,0,1])) +set_exchange!(sys, [J′2apm 0.0 0.0; + 0.0 J′2apm 0.0; + 0.0 0.0 J′2azz], Bond(1,1,[1,2,1])) + +# The function [`set_onsite_coupling!`](@ref) assigns a single-ion anisotropy. +# The argument can be constructed using [`spin_matrices`](@ref) or +# [`stevens_matrices`](@ref). Here we use Julia's anonymous function syntax to +# assign an easy-axis anisotropy along the direction ``\hat{z}``. + +D = 2.165 # (meV) +set_onsite_coupling!(sys, S -> -D*S[3]^2, 1) + +# ### Finding the ground state + +# This model has been carefully designed so that energy minimization produces +# the physically correct magnetic ordering. Using [`set_dipole!`](@ref), this +# magnetic structure can be entered manually. Sunny also provides tools to +# search for an unknown magnetic order, as we will now demonstrate. +# +# To minimize bias in the search, begin with a large system consisting of 4×4×4 +# chemical cells. Randomize all spins (as SU(3) coherent states) and minimize +# the energy. + +sys = resize_supercell(sys, (4, 4, 4)) +randomize_spins!(sys) +minimize_energy!(sys) + +# A positive number above indicates that the procedure has converged to a local +# energy minimum. The configuration, however, may still have defects. This can +# be checked by visualizing the expected spin dipoles, colored according to +# their ``z``-components. + +plot_spins(sys; color=[s[3] for s in sys.dipoles]) + +# To better understand the spin configuration, we could inspect the static +# structure factor ``\mathcal{S}(𝐪)`` in the 3D space of momenta ``𝐪``. For +# this, Sunny provides [`SampledCorrelationsStatic`](@ref). Here, however, we +# will use [`print_wrapped_intensities`](@ref), which gives static intensities +# averaged over the individual Bravais sublattices (in effect, all ``𝐪`` +# intensities are periodically wrapped to the first Brillouin zone). + +print_wrapped_intensities(sys) + +# The known zero-field energy-minimizing magnetic structure of FeI₂ is a two-up, +# two-down order. It can be described as a generalized spiral with a single +# propagation wavevector ``𝐤``. Rotational symmetry allows for three equivalent +# orientations: ``±𝐤 = [0, -1/4, 1/4]``, ``[1/4, 0, 1/4]``, or +# ``[-1/4,1/4,1/4]``. Small systems can spontaneously break this symmetry, but +# for larger systems, defects and competing domains are to be expected. +# Nonetheless, `print_wrapped_intensities` shows large intensities consistent +# with a subset of the known ordering wavevectors. +# +# Let's break the three-fold symmetry by hand. The function +# [`suggest_magnetic_supercell`](@ref) takes one or more ``𝐤`` modes, and +# suggests a magnetic cell shape that is commensurate. + +suggest_magnetic_supercell([[0, -1/4, 1/4]]) + +# Calling [`reshape_supercell`](@ref) yields a much smaller system, making it +# much easier to find the global energy minimum. Plot the system again, now +# including "ghost" spins out to 12Å, to verify that the magnetic order is +# consistent with FeI₂. + +sys_min = reshape_supercell(sys, [1 0 0; 0 2 1; 0 -2 1]) +randomize_spins!(sys_min) +minimize_energy!(sys_min); +plot_spins(sys_min; color=[s[3] for s in sys_min.dipoles], ghost_radius=12) + +# ### Spin wave theory +# +# Now that the system has been relaxed to an energy minimized ground state, we +# can calculate the spin wave spectrum. Because we are working with a system of +# SU(3) coherent states, this calculation will require a multi-flavor boson +# generalization of the usual spin wave theory. + +# Calculate and plot the spectrum along a momentum-space path that connects a +# sequence of high-symmetry ``𝐪``-points. These Sunny commands are identical to +# those described in [`our previous CoRh₂O₄ tutorial`](@ref "1. Spin wave +# simulations of CoRh₂O₄"). + +qs = [[0,0,0], [1,0,0], [0,1,0], [1/2,0,0], [0,1,0], [0,0,0]] +path = q_space_path(cryst, qs, 500) +swt = SpinWaveTheory(sys_min; measure=ssf_perp(sys_min)) +res = intensities_bands(swt, path) +plot_intensities(res; units) + +# To make direct comparison with inelastic neutron scattering (INS) data, we +# must account for empirical broadening of the data. Model this using a +# [`lorentzian`](@ref) kernel, with a full-width at half-maximum of 0.3 meV. + +kernel = lorentzian(fwhm=0.3) +energies = range(0, 10, 300); # 0 < ω < 10 (meV) + +# Also, a real FeI₂ sample will exhibit competing magnetic domains. We use +# [`domain_average`](@ref) to average over the three possible domain +# orientations. This involves 120° rotations about the axis ``\hat{z} = [0, 0, +# 1]`` in global Cartesian coordinates. + +rotations = [([0,0,1], n*(2π/3)) for n in 0:2] +weights = [1, 1, 1] +res = domain_average(cryst, path; rotations, weights) do path_rotated + intensities(swt, path_rotated; energies, kernel) +end +plot_intensities(res; units, colormap=:viridis) + +# This result can be directly compared to experimental neutron scattering data +# from [Bai et al.](https://doi.org/10.1038/s41567-020-01110-1) +# ```@raw html +# +# ``` +# +# (The publication figure used a non-standard coordinate system to label the +# wave vectors.) +# +# To get this agreement, the theory of SU(3) coherent states is essential. The +# lower band has large quadrupolar character, and arises from the strong +# easy-axis anisotropy of FeI₂. +# +# An interesting exercise is to repeat the same study, but using `:dipole` mode +# instead of `:SUN`. That alternative choice would constrain the coherent state +# dynamics to the space of dipoles only, and the flat band of single-ion bound +# states would be missing. diff --git a/examples/04_GSD_FeI2.jl b/examples/04_GSD_FeI2.jl index 6b052dc25..d6f86ec74 100644 --- a/examples/04_GSD_FeI2.jl +++ b/examples/04_GSD_FeI2.jl @@ -1,167 +1,129 @@ # # 4. Generalized spin dynamics of FeI₂ at finite *T* -using Sunny, LinearAlgebra, GLMakie - -# In the [previous FeI₂ tutorial](@ref "1. Multi-flavor spin wave simulations of -# FeI₂ (Showcase)"), we used multi-flavor spin wave theory to calculate the -# dynamical structure factor. Here, we perform a similar calculation using a -# [generalized classical spin dynamics](https://arxiv.org/abs/2209.01265) that -# captures the coupled dynamics of spin dipoles and quadrupoles for -# configurations sampled at finite temperature. +# The [previous FeI₂ tutorial](@ref "3. Multi-flavor spin wave simulations of +# FeI₂") used multi-flavor spin wave theory to calculate the dynamical spin +# structure factor. Here we perform an analogous calculation at finite +# temperature using the [classical dynamics of SU(_N_) coherent +# states](https://doi.org/10.1103/PhysRevB.106.054423). # -# Compared to spin wave theory, simulations using classical dynamics will be -# slower and limited in ``𝐪``-space resolution. However, they make it is -# possible to study [temperature driven phase -# transitions](https://arxiv.org/abs/2310.19905). They may also be used to study -# out-of-equilibrium systems (e.g., relaxation of spin glasses), or systems with -# quenched inhomogeneities that require large simulation volumes. +# Compared to spin wave theory, classical spin dynamics in real-space is +# typically much slower, and is limited in ``𝐪``-space resolution. The +# approach, however, allows for thermal fluctuations, can be used to explore +# [finite temperature phases](https://doi.org/10.1103/PhysRevB.109.014427), and +# enables the study of [highly non-equilibrium +# processes](https://doi.org/10.1103/PhysRevB.106.235154). # -# In this tutorial, we show how to study the finite temperature dynamics of FeI₂ -# using the classical approach. It is important to stress that the estimation of -# ``S(𝐪,ω)`` with classical dynamics is fundamentally a Monte Carlo -# calculation: sample spin configurations are drawn from thermal equilibrium and -# used as initial conditions for generating dissipationless trajectories. The -# correlations of these trajectories are then averaged and used to calculate -# scattering intensities. It is therefore important to ensure that the initial -# spin configurations are sampled appropriately and that sufficient statistics -# are collected. We will demonstrate one approach here. +# The structure of this tutorial largely follows our [previous study of CoRh₂O₄ +# at finite *T*](@ref "2. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T*"). +# The main difference is that CoRh₂O₄ can be well described with `:dipole` mode, +# whereas FeI₂ has a strong easy-axis anisotropy that introduces a single-ion +# bound state and necessitates the use of `:SUN` mode. # -# As an overview, we will: -# -# 1. Identify the ground state. -# 2. Measure correlation data describing the excitations around that ground -# state. -# 3. Use the correlation data to compute scattering intensities. -# -# To begin, please follow our [previous tutorial](@ref "1. Multi-flavor spin -# wave simulations of FeI₂ (Showcase)") to initialize a FeI₂ `sys` with lattice -# dimensions ``4×4×4``. - -a = b = 4.05012#hide -c = 6.75214#hide -latvecs = lattice_vectors(a, b, c, 90, 90, 120)#hide -positions = [[0,0,0], [1/3, 2/3, 1/4], [2/3, 1/3, 3/4]]#hide -types = ["Fe", "I", "I"]#hide -FeI2 = Crystal(latvecs, positions; types)#hide -cryst = subcrystal(FeI2, "Fe")#hide -units = Units(:meV, :angstrom)#hide -sys = System(cryst, (4,4,4), [SpinInfo(1,S=1,g=2)], :SUN, seed=2)#hide -J1pm = -0.236#hide -J1pmpm = -0.161#hide -J1zpm = -0.261#hide -J2pm = 0.026#hide -J3pm = 0.166#hide -J′0pm = 0.037#hide -J′1pm = 0.013#hide -J′2apm = 0.068#hide -J1zz = -0.236#hide -J2zz = 0.113#hide -J3zz = 0.211#hide -J′0zz = -0.036#hide -J′1zz = 0.051#hide -J′2azz = 0.073#hide -J1xx = J1pm + J1pmpm#hide -J1yy = J1pm - J1pmpm#hide -J1yz = J1zpm#hide -set_exchange!(sys, [J1xx 0.0 0.0; 0.0 J1yy J1yz; 0.0 J1yz J1zz], Bond(1,1,[1,0,0]))#hide -set_exchange!(sys, [J2pm 0.0 0.0; 0.0 J2pm 0.0; 0.0 0.0 J2zz], Bond(1,1,[1,2,0]))#hide -set_exchange!(sys, [J3pm 0.0 0.0; 0.0 J3pm 0.0; 0.0 0.0 J3zz], Bond(1,1,[2,0,0]))#hide -set_exchange!(sys, [J′0pm 0.0 0.0; 0.0 J′0pm 0.0; 0.0 0.0 J′0zz], Bond(1,1,[0,0,1]))#hide -set_exchange!(sys, [J′1pm 0.0 0.0; 0.0 J′1pm 0.0; 0.0 0.0 J′1zz], Bond(1,1,[1,0,1]))#hide -set_exchange!(sys, [J′2apm 0.0 0.0; 0.0 J′2apm 0.0; 0.0 0.0 J′2azz], Bond(1,1,[1,2,1]))#hide -D = 2.165#hide -set_onsite_coupling!(sys, S -> -D*S[3]^2, 1)#hide -sys - -# ## Finding a ground state - -# As [previously observed](@ref "1. Multi-flavor spin wave simulations of FeI₂ -# (Showcase)"), direct energy minimization is susceptible to trapping in a local -# energy minimum. - -randomize_spins!(sys) -minimize_energy!(sys) -plot_spins(sys; color=[s[3] for s in sys.dipoles]) - -# Alternatively, one can search for the ordered state by sampling spin -# configurations from thermal equilibrium. Sunny supports this via a -# [`Langevin`](@ref) dynamics of SU(_N_) coherent states. This dynamics involves -# a dimensionless `damping` magnitude and target temperature `kT` for thermal -# fluctuations. - -kT = 0.2 # Temperature in meV -langevin = Langevin(; damping=0.2, kT) - -# Use [`suggest_timestep`](@ref) to select an integration timestep for the given -# error tolerance, e.g. `tol=1e-2`. The spin configuration in `sys` should -# ideally be relaxed into thermal equilibrium, but the current, energy-minimized -# configuration will also work reasonably well. - -suggest_timestep(sys, langevin; tol=1e-2) -langevin.dt = 0.027; +# Construct the FeI₂ system as described in the [previous tutorial](@ref "3. +# Multi-flavor spin wave simulations of FeI₂"). + +using Sunny, GLMakie + +units = Units(:meV, :angstrom) +a = b = 4.05012 +c = 6.75214 +latvecs = lattice_vectors(a, b, c, 90, 90, 120) +cryst = Crystal(latvecs, [[0,0,0]], 164; types=["Fe"]) + +sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=1, g=2)], :SUN) +J1pm = -0.236 +J1pmpm = -0.161 +J1zpm = -0.261 +J2pm = 0.026 +J3pm = 0.166 +J′0pm = 0.037 +J′1pm = 0.013 +J′2apm = 0.068 +J1zz = -0.236 +J2zz = 0.113 +J3zz = 0.211 +J′0zz = -0.036 +J′1zz = 0.051 +J′2azz = 0.073 +J1xx = J1pm + J1pmpm +J1yy = J1pm - J1pmpm +J1yz = J1zpm +set_exchange!(sys, [J1xx 0.0 0.0; 0.0 J1yy J1yz; 0.0 J1yz J1zz], Bond(1,1,[1,0,0])) +set_exchange!(sys, [J2pm 0.0 0.0; 0.0 J2pm 0.0; 0.0 0.0 J2zz], Bond(1,1,[1,2,0])) +set_exchange!(sys, [J3pm 0.0 0.0; 0.0 J3pm 0.0; 0.0 0.0 J3zz], Bond(1,1,[2,0,0])) +set_exchange!(sys, [J′0pm 0.0 0.0; 0.0 J′0pm 0.0; 0.0 0.0 J′0zz], Bond(1,1,[0,0,1])) +set_exchange!(sys, [J′1pm 0.0 0.0; 0.0 J′1pm 0.0; 0.0 0.0 J′1zz], Bond(1,1,[1,0,1])) +set_exchange!(sys, [J′2apm 0.0 0.0; 0.0 J′2apm 0.0; 0.0 0.0 J′2azz], Bond(1,1,[1,2,1])) +D = 2.165 +set_onsite_coupling!(sys, S -> -D*S[3]^2, 1) + +# ### Thermalization + +# To study thermal fluctuations in real-space, use a large system size with +# 16×16×4 copies of the chemical cell. + +sys_large = resize_supercell(sys, (16, 16, 4)) + +# Previously, we used [`minimize_energy!`](@ref) to find a local energy minimum. +# Here, we will instead use [`Langevin`](@ref) dynamics to relax the system into +# thermal equilibrium. The temperature 2.3 K ≈ 0.2 meV is within the ordered +# phase, but large enough so that the dynamics can overcome local energy +# barriers and annihilate defects. + +langevin = Langevin(; damping=0.2, kT=2.3*units.K) + +# Use [`suggest_timestep`](@ref) to select an integration timestep for the error +# tolerance `tol=1e-2`. Initializing `sys` to some low-energy configuration +# usually works well. + +randomize_spins!(sys_large) +minimize_energy!(sys_large; maxiters=10) +suggest_timestep(sys_large, langevin; tol=1e-2) +langevin.dt = 0.03; -# Sample spin configurations using Langevin dynamics. We have carefully selected -# a temperature of 0.2 eV that is below the ordering temperature, but large -# enough to that the dynamics can overcome local energy barriers and annihilate -# defects. +# Run a Langevin trajectory for 10,000 time-steps and plot the spins. At this +# angle, it is difficult to discern the magnetic order. for _ in 1:10_000 - step!(sys, langevin) + step!(sys_large, langevin) end +plot_spins(sys_large; color=[s[3] for s in sys_large.dipoles]) -# Calling [`suggest_timestep`](@ref) shows that thermalization has not -# substantially altered the suggested `dt`. - -suggest_timestep(sys, langevin; tol=1e-2) - -# Although thermal fluctuations are present, the correct antiferromagnetic order -# (2 up, 2 down) has been found. +# The two-up, two-down spiral order can be verified by calling +# [`print_wrapped_intensities`](@ref). A single propagation wavevector ``±𝐤`` +# provides most of the static intensity in ``\mathcal{S}(𝐪)``. A smaller amount +# of intensity is spread among many other wavevectors due to thermal +# fluctuations. -plot_spins(sys; color=[s[3] for s in sys.dipoles]) +print_wrapped_intensities(sys_large) -# For other phases, it can be much harder to find thermal equilibrium, and more -# complicated sampling procedures may be necessary. +# Calling [`suggest_timestep`](@ref) shows that thermalization has not +# substantially altered the suggested `dt`. -# ## Calculating Thermal-Averaged Correlations $⟨S^{αβ}(𝐪,ω)⟩$ -# -# Our aim is to study the classical spin dynamics for states sampled in thermal -# equilibrium. To minimize finite size effects, and achieve sufficient momentum -# space resolution, we should significantly enlarge the system volume. The -# function [`resize_supercell`](@ref) takes new dimensions as multiples of the -# unit cell lattice vectors. +suggest_timestep(sys_large, langevin; tol=1e-2) -sys_large = resize_supercell(sys, (16,16,4)) # 16x16x4 copies of the original unit cell -plot_spins(sys_large; color=[s[3] for s in sys_large.dipoles]) +# ### Structure factor in the paramagnetic phase -# Now we will re-thermalize the system to a configuration just above the -# ordering temperature. +# Now we will re-thermalize the system to a temperature of 3.5 K ≈ 0.30 meV, +# which is in the paramagnetic phase. -kT = 3.5 * units.K # 3.5 K ≈ 0.30 meV -langevin.kT = kT +langevin.kT = 3.5 * units.K for _ in 1:10_000 step!(sys_large, langevin) end -# With this increase in temperature, the suggested timestep has increased slightly. +# At this higher temperature, the suggested timestep has increased slightly. suggest_timestep(sys_large, langevin; tol=1e-2) langevin.dt = 0.040; -# The next step is to collect correlation data ``S^{αβ}`` into a +# Now collect dynamical spin structure factor data using a # [`SampledCorrelations`](@ref) object. This will involve sampling spin -# configurations from thermal equilibrium, and then integrating [an -# energy-conserving generalized classical spin -# dynamics](https://arxiv.org/abs/2204.07563) to collect Fourier-space -# information about normal modes. Quantization of these modes yields the -# magnons, and the associated dynamical spin-spin correlations can be compared -# with neutron scattering intensities ``S^{αβ}(q,ω)``. Because this is a -# real-space calculation, data is only available for discrete ``q`` modes (the -# resolution scales like inverse system size). -# -# The `SampledCorrelations` object requires specification of an integration -# timestep `dt`, an energy range for the intensity measurements, and a choice of -# pair correlation measurement. A rule of thumb is that the `dt` can be twice -# larger than the Langevin integration timestep. +# configurations from thermal equilibrium and integrating a [classical spin +# dynamics for SU(_N_) coherent states](https://arxiv.org/abs/2204.07563). +# Normal modes appearing in the classical dynamics can be quantized to yield +# magnetic excitations. The associated structure factor intensities +# ``S^{αβ}(q,ω)`` can be compared with inelastic neutron scattering data . dt = 2*langevin.dt energies = range(0, 7.5, 120) @@ -183,31 +145,25 @@ for _ in 1:2 add_sample!(sc, sys_large) end -# Now, `sc` has more samples included: - -sc - -# ## Computing Scattering Intensities +# Measure intensities along a path connecting high-symmetry ``𝐪``-points, +# specified in reciprocal lattice units (RLU). Incorporate the +# [`FormFactor`](@ref) appropriate to Fe²⁺. A classical-to-quantum rescaling of +# normal mode occupations will be performed according to the temperature `kT`. +# The large statistical noise could be reduced by averaging over more thermal +# samples. -# Extract the thermally-averaged correlation data ``⟨S^{αβ}(q,ω)⟩`` for a given -# set of q-points using [`intensities`](@ref). A classical-to-quantum rescaling -# of normal mode occupations will be performed according to the temperature -# `kT`. - -qs = Sunny.QPoints([[0, 0, 0], [0.5, 0.5, 0.5]]) -res = intensities(sc, qs; kT, energies) - -fig = lines(res.energies, res.data[:,1]; axis=(xlabel="meV", ylabel="Intensity"), label="(0,0,0)") -lines!(res.energies, res.data[:,2]; label="(π,π,π)") +formfactors = [FormFactor("Fe2"; g_lande=3/2)] +res = intensities(sc, [[0, 0, 0], [0.5, 0.5, 0.5]]; energies, formfactors, langevin.kT) +fig = lines(res.energies, res.data[:, 1]; axis=(xlabel="meV", ylabel="Intensity"), label="(0,0,0)") +lines!(res.energies, res.data[:, 2]; label="(π,π,π)") axislegend() fig -# Significant fluctuations indicate a large stochastic error. This data can be -# made smooth by iterating over more calls to `add_sample!` iterations. -# -# Next, we will measure intensities along the [`q_space_path`](@ref) that -# connects high symmetry points. Here we will also apply a [`FormFactor`](@ref) -# appropriate to Fe²⁺. +# Next, we will measure intensities along a [`q_space_path`](@ref) that connects +# high symmetry points. Because this is a real-space calculation, data is only +# available for discrete ``𝐪`` modes, with resolution that scales inversely to +# linear system size. Intensities at ``ω = 0`` dominate, so to enhance +# visibility, we restrict the color range empirically. qs = [[0, 0, 0], # List of wave vectors that define a path [1, 0, 0], @@ -216,22 +172,13 @@ qs = [[0, 0, 0], # List of wave vectors that define a path [0, 1, 0], [0, 0, 0]] qpath = q_space_path(cryst, qs, 500) -formfactors = [FormFactor("Fe2"; g_lande=3/2)] -res = intensities(sc, qpath; kT, energies, formfactors) +res = intensities(sc, qpath; energies, formfactors, langevin.kT) plot_intensities(res; colorrange=(0.0, 1.0)) -# On can also view the intensity along a ``𝐪``-space slice at a fixed energy -# value. +# One can also view the intensity along a [`q_space_grid`](@ref) for a fixed +# energy value. Alternatively, use [`intensities_static`](@ref) to integrate +# over all available energies. grid = q_space_grid(cryst, [1, 0, 0], range(-1.5, 1.5, 300), [0, 1, 0], (-1.5, 1.5); orthogonalize=true) -res = intensities(sc, grid; energies=[3.88], kT) -plot_intensities(res; colorrange=(0.0, 0.3)) - -# Significant ``𝐪``-space discretization is apparent. This is a consequence of -# the finite system size used to simulate the classical dynamics. The basic -# intensity data stored inside `SampledCorrelations` is only available on a -# discrete grid of ``𝐪``-points, with resolution that scales inversely to -# linear system size. Sunny has an internal function to investigate this grid as -# a 3D array (subject to change): - -size(Sunny.available_wave_vectors(sc)) +res = intensities(sc, grid; energies=[3.88], formfactors, langevin.kT) +plot_intensities(res) diff --git a/examples/05_MC_Ising.jl b/examples/05_MC_Ising.jl index 5143f1a7f..f2a32be20 100644 --- a/examples/05_MC_Ising.jl +++ b/examples/05_MC_Ising.jl @@ -4,34 +4,30 @@ using Sunny, GLMakie -# Sunny expects a 3D [`Crystal`](@ref) unit cell. To model a square lattice, we -# create an orthogonal unit cell where the ``z``-spacing is distinct from the -# ``x`` and ``y`` spacing. +# [`Crystal`](@ref) unit cell are always 3D. To model a square lattice, we +# create a tetragonal cell with one atom and an elongated lattice constant +# ``c``. + a = 1 -latvecs = lattice_vectors(a,a,10a,90,90,90) -crystal = Crystal(latvecs, [[0,0,0]]) - -# Create a [`System`](@ref) of spins with linear size `L` in the ``x`` and ``y`` -# directions, and only one layer in the ``z`` direction. The option `:dipole` -# means that the system will store Heisenberg spins, as opposed to SU(``N``) -# coherent states. Polarize the initial spin configuration using -# [`polarize_spins!`](@ref). Following the Ising convention, we will restrict -# these spins to the ``z``-axis and give them magnitude ``S=1``. -# -# By default, Sunny expects the magnetic field in tesla. Selecting -# [`Units.theory`](@ref Units) instead allows for dimensionless units. Following -# Ising conventions, select ``g=-1`` so that the Zeeman coupling between -# external field ``𝐁`` and spin dipole ``𝐬`` is ``-𝐁⋅𝐬``. +latvecs = lattice_vectors(a, a, 10a, 90, 90, 90) +crystal = Crystal(latvecs, [[0, 0, 0]]) + +# Create a [`System`](@ref) of spin dipoles. Following the Ising convention, we +# will restrict the dipoles to ``±1`` along the global ``\hat{z}``-axis. Select +# ``g=-1`` so that the Zeeman coupling between external field ``𝐁`` and spin +# dipole ``𝐬`` is ``-𝐁⋅𝐬``. The initial supercell size is ``L×L``. + L = 128 -sys = System(crystal, (L,L,1), [SpinInfo(1, S=1, g=-1)], :dipole, seed=0) -polarize_spins!(sys, (0,0,1)) +sys = System(crystal, (1, 1, 1), [SpinInfo(1, S=1, g=-1)], :dipole, seed=0) +sys = resize_supercell(sys, (L, L, 1)) +polarize_spins!(sys, (0, 0, 1)) # Use [`set_exchange!`](@ref) to include a ferromagnetic Heisenberg interaction # along nearest-neighbor bonds. The [`Bond`](@ref) below connects two spins -# displaced by one lattice constant in the ``x``-direction. This interaction -# will be propagated to all nearest-neighbors bonds in the system, consistent -# with the symmetries of the square lattice. -set_exchange!(sys, -1.0, Bond(1,1,(1,0,0))) +# displaced by the lattice vector ``𝐚₁``. This interaction will be propagated +# to all nearest-neighbors bonds in the system, consistent with the symmetries +# of the square lattice. +set_exchange!(sys, -1.0, Bond(1, 1, (1, 0, 0))) # If an external field is desired, it can be set using [`set_field!`](@ref). B = 0 @@ -52,4 +48,4 @@ for i in 1:nsweeps end # Plot the Ising spins by extracting the ``z``-component of the dipoles -heatmap(reshape([s.z for s in sys.dipoles], (L,L))) +heatmap(reshape([s[3] for s in sys.dipoles], (L, L))) diff --git a/src/Optimization.jl b/src/Optimization.jl index 7557970c8..28051022b 100644 --- a/src/Optimization.jl +++ b/src/Optimization.jl @@ -140,9 +140,8 @@ function minimize_energy!(sys::System{N}; maxiters=1000, subiters=10, method=Opt αs .*= 0 end - res_str = number_to_simple_string(g_res; digits=2) - tol_str = number_to_simple_string(g_tol; digits=2) - - @warn "Optimization failed to converge within $maxiters iterations ($res_str ≰ $tol_str)" + # res_str = number_to_simple_string(g_res; digits=2) + # tol_str = number_to_simple_string(g_tol; digits=2) + # @warn "Optimization failed to converge within $maxiters iterations ($res_str ≰ $tol_str)" return -1 end diff --git a/src/System/System.jl b/src/System/System.jl index 8b0e8d7e8..7ddbe9e30 100644 --- a/src/System/System.jl +++ b/src/System/System.jl @@ -21,7 +21,7 @@ see the documentation page: [Interaction Strength Renormalization](@ref). An optional `seed` may be provided to achieve reproducible random number generation. -All spins are initially polarized in the ``z``-direction. +All spins are initially polarized in the global ``z``-direction. """ function System(crystal::Crystal, latsize::NTuple{3,Int}, infos::Vector{SpinInfo}, mode::Symbol; seed=nothing, units=nothing) @@ -93,8 +93,7 @@ function mode_to_str(sys::System{N}) where N end function lattice_to_str(sys::System) - lat_str = isnothing(sys.origin) ? "Lattice" : "Reshaped lattice" - return lat_str * " ($(join(sys.latsize, "×")))×$(natoms(sys.crystal))" + return "Lattice (" * join(sys.latsize, "×") * ")×" * string(natoms(sys.crystal)) end function energy_to_str(sys::System) @@ -112,7 +111,7 @@ end function Base.show(io::IO, ::MIME"text/plain", sys::System{N}) where N printstyled(io, "System $(mode_to_str(sys))\n"; bold=true, color=:underline) println(io, lattice_to_str(sys)) - if !isnothing(sys.origin) + if !isnothing(sys.origin) && cell_shape(sys) != cell_shape(sys.origin) shape = number_to_math_string.(cell_shape(sys)) println(io, formatted_matrix(shape; prefix="Reshaped cell ")) end