Skip to content

Commit

Permalink
connected_path_from_rlu
Browse files Browse the repository at this point in the history
  • Loading branch information
kbarros committed Aug 4, 2023
1 parent 3410dec commit 4962c5b
Show file tree
Hide file tree
Showing 10 changed files with 76 additions and 68 deletions.
6 changes: 3 additions & 3 deletions docs/src/structure-factor.md
Original file line number Diff line number Diff line change
Expand Up @@ -170,9 +170,9 @@ e.g., `bzsize=(2,2,2)`. The resulting list of wave vectors may then be passed to
Alternatively, [`intensities_binned`](@ref) can be used to place the exact data
into histogram bins for comparison with experiment.

The convenience function [`connected_path`](@ref) returns a list of wavevectors
sampled along a path that connects specified $𝐪$ points. This list can be used
as an input to `intensities`. Another convenience method,
The convenience function [`connected_path_from_rlu`](@ref) returns a list of
wavevectors sampled along a path that connects specified $𝐪$ points. This list
can be used as an input to `intensities`. Another convenience method,
[`spherical_shell`](@ref) will provide a list of wave vectors on a sphere of a
specified radius. This is useful for powder averaging.

Expand Down
32 changes: 11 additions & 21 deletions examples/fei2_tutorial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -245,29 +245,19 @@ plot_spins(sys_min; arrowlength=2.5, linewidth=0.75, arrowsize=1.5)

swt = SpinWaveTheory(sys_min);

# Next, select a sequence of wavevectors that will define a piecewise linear
# interpolation in reciprocal lattice units (RLU). Sunny expects wavevectors in
# absolute units, i.e., inverse angstroms. To convert from RLU, multiply by the
# $3×3$ matrix `cryst.recipvecs`, which provides the reciprocal lattice vectors
# for the conventional unit cell.
# Select a sequence of wavevectors that will define a piecewise linear
# interpolation in reciprocal lattice units (RLU).

points_rlu = [[0,0,0], [1,0,0], [0,1,0], [1/2,0,0], [0,1,0], [0,0,0]]
points = [cryst.recipvecs*p for p in points_rlu]

# The function [`connected_path`](@ref) will linearly sample between the
# provided $q$-points with given `density`. Also keep track of location and
# names of the special $𝐪$-points for plotting purposes.
# The function [`connected_path_from_rlu`](@ref) will linearly sample between
# the provided $q$-points with a given `density`. The `path` return value is a
# list of wavevectors in absolute units (inverse Å). The `xticks` return value
# keeps track of the locations of the special $𝐪$-points, and provides
# human-readable labels for use in plotting.

density = 50
path, numbers = connected_path(points, density);
xticks = (numbers, ["[0,0,0]", "[1,0,0]", "[0,1,0]", "[1/2,0,0]", "[0,1,0]", "[0,0,0]"])

# The dispersion relation can be calculated by providing `dispersion` with a
# `SpinWaveTheory` and a list of wave vectors. For each wave vector,
# `dispersion` will return a list of energies, one for each band. Frequently one
# wants dispersion information along lines that connect special wave vectors.
# The function [`connected_path`](@ref) linearly samples between provided
# $𝐪$-points, with a given sample density.
path, xticks = connected_path_from_rlu(cryst, points_rlu, density);

# The [`dispersion`](@ref) function defines the quasiparticle excitation
# energies $ω_i(𝐪)$ for each point $𝐪$ along the reciprocal space path.
Expand Down Expand Up @@ -332,13 +322,13 @@ heatmap!(ax, 1:size(is_averaged, 1), energies, is_averaged)
fig

# This result can be directly compared to experimental neutron scattering data
# [Bai et al.](https://doi.org/10.1038/s41567-020-01110-1)
# from [Bai et al.](https://doi.org/10.1038/s41567-020-01110-1)
# ```@raw html
# <img src="https://raw.githubusercontent.com/SunnySuite/Sunny.jl/main/docs/src/assets/FeI2_intensity.jpg">
# ```
#
# (Note that the publication figure uses a different coordinate system to label
# the wave vectors.)
# (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
Expand Down
8 changes: 4 additions & 4 deletions examples/longer_examples/CoRh2O4-tutorial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -108,7 +108,7 @@ end
# ### System Definition for CoRh<sub>2</sub>O<sub>4</sub>

# Define the crystal structure of CoRh$_2$O$_4$ in the conventional cell
xtal = Crystal(cif_path ;symprec=1e-4)
xtal = Crystal(cif_path; symprec=1e-4)
magxtal = subcrystal(xtal,"Co1")
view_crystal(magxtal,6.0)
print_symmetry_table(magxtal, 4.0)
Expand Down Expand Up @@ -214,7 +214,7 @@ symQpts = [[0.75, 0.75, 0.00], # List of wave vectors that define a path
[0.25, 1.00, 0.25],
[0.00, 1.00, 0.00],
[0.00,-4.00, 0.00]];
(Qpts, symQind) = connected_path(sc, symQpts,densQpts);
(Qpts, xticks) = connected_path_from_rlu(magxtal, symQpts, densQpts);
@time iqw = intensities(sc, Qpts, :perp; interpolation = :none, kT=integrator.kT, formfactors=ffs);

# If desired, broaden the sc in energy
Expand All @@ -227,10 +227,10 @@ iqwc = broaden_energy(sc, iqw, (ω, ω₀) -> lorentzian(ω-ω₀, η));
# Plot the resulting I(Q,W)
heatmap(1:size(iqwc, 1), ωs(sc), iqwc;
colorrange = (0, maximum(iqwc)/20000.0),
axis = (
axis = (;
xlabel="Momentum Transfer (r.l.u)",
ylabel="Energy Transfer (meV)",
xticks = (symQind, string.(symQpts)),
xticks,
xticklabelrotation=π/5,
aspect = 1.4,
)
Expand Down
4 changes: 2 additions & 2 deletions examples/longer_examples/MgCr2O4-tutorial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -288,7 +288,7 @@ axsqw = (xticks=-4:4, yticks=0:2:10, ylabel="E (meV)", ylabelsize=18, xlabelsize

qbs = 0.0:0.5:1.5 # Determine q_b for each slice
for (i, qb) in enumerate(qbs)
path, labels = connected_path(sc_pyro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40) # Generate a path of wave
path, _ = connected_path_from_rlu(xtal_pyro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40) # Generate a path of wave
## vectors through the BZ
Sqω_pyro = intensities(sc_pyro, path, :perp; kT) # Temperature keyword enables intensity rescaling

Expand All @@ -303,7 +303,7 @@ fig = Figure(; resolution=(1200,900))

qbs = 0.0:0.5:1.5
for (i, qb) in enumerate(qbs)
path, labels = connected_path(sc_mgcro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40) # Generate a path of wave
path, _ = connected_path_from_rlu(xtal_mgcro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40) # Generate a path of wave
## vectors through the BZ
Sqω_mgcro = intensities(sc_mgcro, path, :perp; kT) # Temperature keyword enables intensity rescaling

Expand Down
2 changes: 1 addition & 1 deletion examples/powder_averaging.jl
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ add_sample!(sc, sys)
# lorentzian(ω-ω₀, 0.1)`.

qpoints = [[0.0, 0.0, 0.0], [0.5, 0.0, 0.0], [0.5, 0.5, 0.0], [0.0, 0.0, 0.0]]
qs, markers = connected_path(sc, qpoints, 50)
qs, markers = connected_path_from_rlu(crystal, qpoints, 50)

is = intensities_interpolated(sc, qs; interpolation=:round, formula = intensity_formula(sc,:trace))
is_broad = broaden_energy(sc, is, (ω, ω₀) -> lorentzian-ω₀, 0.1))
Expand Down
35 changes: 23 additions & 12 deletions src/Intensities/Interpolation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -218,18 +218,24 @@ end


"""
connected_path(qs, density)
Returns two things: First, a list of wavevectors that linearly interpolates
between the points `qs` and second, a list of markers that index the locations
of each original interpolation point.
More sample-points will be returned between `q`-values that are further apart,
and this is controlled by the `density` parameter.
connected_path_from_rlu(cryst, qs_rlu, density)
Returns a pair `(path, xticks)`. The first return value is a path in reciprocal
space that samples linearly between the wavevectors in `qs`. The elements in
`qs` are defined in reciprocal lattice units (RLU) associated with the
[`reciprocal_lattice_vectors``](@ref) for `cryst`. The `density` parameter has
units of inverse length, and controls the number of samples between elements of
`qs`.
The second return value `xticks` can be used for plotting. The `xticks` object
is itself a pair `(numbers, labels)`, which give the locations of the
interpolating ``q``-points and a human-readable string.
"""
function connected_path(qs::Vector, density)
@assert length(qs) >= 2 "The list `qs` should include at least two wavevectors."
qs = Vec3.(qs)
function connected_path_from_rlu(cryst::Crystal, qs_rlu::Vector, density)
@assert length(qs_rlu) >= 2 "The list `qs` should include at least two wavevectors."
qs_rlu = Vec3.(qs_rlu)

qs = Ref(cryst.recipvecs) .* qs_rlu

path = Vec3[]
markers = Int[]
Expand All @@ -245,5 +251,10 @@ function connected_path(qs::Vector, density)
push!(markers, length(path)+1)
push!(path, qs[end])

return (path, markers)
labels = map(qs_rlu) do q_rlu
"[" * join(number_to_math_string.(q_rlu), ",") * "]"
end
xticks = (markers, labels)

return (path, xticks)
end
4 changes: 1 addition & 3 deletions src/SpinWaveTheory/SWTCalculations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -571,9 +571,7 @@ function intensity_formula(f::Function, swt::SpinWaveTheory, corr_ix; kernel::Un
#
# Smooth kernel --> I_of_ω = Intensity as a function of ω
#
calc_intensity = function(swt::SpinWaveTheory,q::Vec3)
_, qmag = chemical_to_magnetic(swt, q)

calc_intensity = function(swt::SpinWaveTheory, q::Vec3)
if sys.mode == :SUN
swt_hamiltonian_SUN!(swt, q, Hmat)
elseif sys.mode == :dipole
Expand Down
7 changes: 3 additions & 4 deletions src/Sunny.jl
Original file line number Diff line number Diff line change
Expand Up @@ -49,9 +49,8 @@ include("Symmetry/AllowedCouplings.jl")
include("Symmetry/AllowedAnisotropy.jl")
include("Symmetry/Parsing.jl")
include("Symmetry/Printing.jl")
export Crystal, subcrystal, lattice_vectors, lattice_params, Bond,
reference_bonds, print_site, print_bond, print_symmetry_table,
print_suggested_frame
export Crystal, subcrystal, lattice_vectors, lattice_params, reciprocal_lattice_vectors, Bond,
reference_bonds, print_site, print_bond, print_symmetry_table, print_suggested_frame

include("Units.jl")
export meV_per_K, Units
Expand Down Expand Up @@ -98,7 +97,7 @@ export SampledCorrelations, dynamical_correlations, instant_correlations, FormFa
include("Intensities/ElementContraction.jl")

include("Intensities/Interpolation.jl")
export intensities_interpolated, instant_intensities_interpolated, connected_path
export intensities_interpolated, instant_intensities_interpolated, connected_path_from_rlu

include("Intensities/Binning.jl")
export intensities_binned, BinningParameters, count_bins, integrate_axes!,
Expand Down
42 changes: 26 additions & 16 deletions src/Symmetry/Crystal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,20 +81,6 @@ struct Crystal
symprec :: Float64 # Tolerance to imperfections in symmetry
end

"""
natoms(cryst::Crystal)
Number of atoms in the unit cell, i.e., number of Bravais sublattices.
"""
@inline natoms(cryst::Crystal) = length(cryst.positions)

"""
cell_volume(cryst::Crystal)
Volume of the crystal unit cell.
"""
cell_volume(cryst::Crystal) = abs(det(cryst.latvecs))

# Constructs a crystal from the complete list of atom positions `positions`,
# representing fractions (between 0 and 1) of the lattice vectors `latvecs`.
# All symmetry information is automatically inferred.
Expand Down Expand Up @@ -145,6 +131,30 @@ function print_crystal_warnings(latvecs, positions)
end
end

"""
natoms(cryst::Crystal)
Number of atoms in the unit cell, i.e., number of Bravais sublattices.
"""
@inline natoms(cryst::Crystal) = length(cryst.positions)

"""
cell_volume(cryst::Crystal)
Volume of the crystal unit cell.
"""
cell_volume(cryst::Crystal) = abs(det(cryst.latvecs))

"""
reciprocal_lattice_vectors(cryst::Crystal)
Returns a ``3×3`` matrix, with columns that define the reciprocal lattice
vectors ``(𝐛₁,𝐛₂,𝐛₃)``. These are defined to satisfy ``𝐛ᵢ⋅𝐚ⱼ = 2πδᵢⱼ``,
where ``(𝐚₁,𝐚₂,𝐚₃)`` are the lattice vectors for `cryst` in real-space.
"""
reciprocal_lattice_vectors(cryst::Crystal) = cryst.recipvecs


function spacegroup_name(hall_number::Int)
# String representation of space group
sgt = Spglib.get_spacegroup_type(hall_number)
Expand Down Expand Up @@ -501,8 +511,8 @@ function subcrystal(cryst::Crystal, classes::Vararg{Int, N}) where N
@info "Atoms have been renumbered in subcrystal."
end

ret = Crystal(cryst.latvecs, cryst.prim_latvecs, new_positions, new_types, new_classes, new_sitesyms,
cryst.symops, cryst.spacegroup, cryst.symprec)
ret = Crystal(cryst.latvecs, cryst.prim_latvecs, cryst.recipvecs, new_positions, new_types,
new_classes, new_sitesyms, cryst.symops, cryst.spacegroup, cryst.symprec)
return ret
end

Expand Down
4 changes: 2 additions & 2 deletions test/test_correlation_sampling.jl

Large diffs are not rendered by default.

0 comments on commit 4962c5b

Please sign in to comment.