Add field cryst.recipvecs

examples/fei2_tutorial.jl

@@ -246,12 +246,18 @@ 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 space. Use [`reciprocal_vectors`](@ref) convert
-# from reciprocal lattice units (RLU) to absolute units in inverse angstroms.
+# 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.
+
+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`.
-rv = reciprocal_vectors(cryst)
-points = [rv*p for p in [[0,0,0], [1,0,0], [0,1,0], [1/2,0,0], [0,1,0], [0,0,0]]]
+# provided $q$-points with given `density`. Also keep track of location and
+# names of the special $𝐪$-points for plotting purposes.
+
 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]"])

src/Intensities/Binning.jl

@@ -131,6 +131,7 @@ function binning_parameters_aabb(params)
     return lower_aabb_q, upper_aabb_q
 end
 
+# FIXME: various conversions go away. Always absolute units.
 """
 If `params` expects to bin values `(k,ω)` in absolute units, then calling
 
@@ -247,6 +248,7 @@ function unit_resolution_binning_parameters(ωvals::AbstractVector{Float64})
     return ωstart, ωend, ωbinwidth
 end
 
+# FIXME: remove `units` arg -- always absolute
 """
     slice_2D_binning_parameter(sc::SampledCorrelations, cut_from_q, cut_to_q, cut_bins::Int64, cut_width::Float64; plane_normal = [0,0,1],cut_height = cutwidth, units = :absolute)
 
@@ -408,7 +410,6 @@ function intensities_binned(sc::SampledCorrelations, params::BinningParameters;
     output_intensities = zeros(Float64,numbins...)
     output_counts = zeros(Float64,numbins...)
     ωvals = ωs(sc)
-    recip_vecs = reciprocal_vectors(sc.crystal)
 
     # Find an axis-aligned bounding box containing the histogram.
     # The AABB needs to be in q-space because that's where we index
@@ -452,7 +453,7 @@ function intensities_binned(sc::SampledCorrelations, params::BinningParameters;
         # Which is the analog of this scattering mode in the first BZ?
         base_cell = CartesianIndex(mod1.(cell.I,Ls)...)
         q .= ((cell.I .- 1) ./ Ls) # q is in R.L.U.
-        k .= recip_vecs * q # But binning is done in absolute units
+        k .= sc.crystal.recipvecs * q # FIXME -- no scaling factor needed, q will already in absolute units.
         for (iω,ω) in enumerate(ωvals)
             if isnothing(integrated_kernel) # `Delta-function energy' logic
                 # Figure out which bin this goes in

src/Intensities/Interpolation.jl

@@ -106,13 +106,12 @@ function pruned_stencil_info(sc::SampledCorrelations, qs, interp::InterpolationS
     @assert sum(counts) == length(m_info)
 
     # Calculate corresponding q (RLU) and k (global) vectors
-    recip_vecs = reciprocal_vectors(sc.crystal) # Note, qs will be in terms of sc.crystal by this point, not origin_crystal
     qs_all = map(ms_all) do ms
        map(m -> m ./ sc.latsize, ms) 
     end
 
     ks_all = map(qs_all) do qs
-        map(q -> recip_vecs * q, qs)
+        map(q -> sc.crystal.recipvecs * q, qs) # FIXME: should q still by in sc.crystal, or in absolute units?
     end
 
     return (; qs_all, ks_all, idcs_all, counts)

src/Intensities/PowderAveraging.jl

@@ -77,12 +77,11 @@ function powder_average_binned(sc::SampledCorrelations, radial_binning_parameter
     output_intensities = zeros(Float64,r_bin_count,ω_bin_count)
     output_counts = zeros(Float64,r_bin_count,ω_bin_count)
     ωvals = ωs(sc)
-    recip_vecs = reciprocal_vectors(sc.crystal)
 
     # Loop over every scattering vector
     Ls = sc.latsize
     if isnothing(bzsize)
-        bzsize = (1,1,1) .* ceil(Int64,rend/eigmin(recip_vecs))
+        bzsize = (1,1,1) .* ceil(Int64,rend/eigmin(sc.crystal.recipvecs)) # TODO: ceil(Int64, a/b) -> div(a, b, RoundUp)
     end
     for cell in CartesianIndices(Ls .* bzsize)
         base_cell = CartesianIndex(mod1.(cell.I,Ls)...)
@@ -93,7 +92,7 @@ function powder_average_binned(sc::SampledCorrelations, radial_binning_parameter
 
             # Figure out which radial bin this scattering vector goes in
             # The spheres are surfaces of fixed |k|, with k in absolute units
-            k = recip_vecs * q
+            k = sc.crystal.recipvecs * q # FIXME: q should already be in absolute units
             r_coordinate = norm(k) 
 
             # Check if the radius falls within the histogram

src/SampledCorrelations/SampledCorrelations.jl

@@ -9,7 +9,7 @@ struct SampledCorrelations{N}
     # 𝒮^{αβ}(q,ω) data and metadata
     data           :: Array{ComplexF64, 7}   # Raw SF data for 1st BZ (numcorrelations × natoms × natoms × latsize × energy)
     crystal        :: Crystal                # Crystal for interpretation of q indices in `data`
-    origin_crystal :: Union{Nothing,Crystal} # Original user-specified crystal (if different from above)
+    origin_crystal :: Union{Nothing,Crystal} # Original user-specified crystal (if different from above)  # FIXME: Eliminate
     Δω             :: Float64                # Energy step size (could make this a virtual property)
 
     # Correlation info (αβ indices of 𝒮^{αβ}(q,ω))

src/SpinWaveTheory/SpinWaveTheory.jl

@@ -156,6 +156,5 @@ function SpinWaveTheory(sys::System{N}; energy_ϵ::Float64=1e-8, energy_tol::Flo
     end
 
     positions = [global_position(sys, site) for site in all_sites(sys)][:]
-    recipvecs = reciprocal_vectors(sys.crystal)
-    return SpinWaveTheory(sys, s̃_mat, T̃_mat, Q̃_mat, c′_coef, R_mat, positions, recipvecs, energy_ϵ, energy_tol)
+    return SpinWaveTheory(sys, s̃_mat, T̃_mat, Q̃_mat, c′_coef, R_mat, positions, sys.crystal.recipvecs, energy_ϵ, energy_tol)
 end

src/Sunny.jl

@@ -49,8 +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, reciprocal_vectors,
-    Bond, reference_bonds, print_site, print_bond, print_symmetry_table,
+export Crystal, subcrystal, lattice_vectors, lattice_params, Bond,
+    reference_bonds, print_site, print_bond, print_symmetry_table,
     print_suggested_frame
 
 include("Units.jl")

src/Symmetry/Crystal.jl

@@ -69,8 +69,9 @@ cryst = Crystal(latvecs, positions, 227; setting="1")
 See also [`lattice_vectors`](@ref).
 """
 struct Crystal
-    latvecs        :: Mat3                                 # Lattice vectors as columns
+    latvecs        :: Mat3                                 # Lattice vectors as columns (conventional)
     prim_latvecs   :: Mat3                                 # Primitive lattice vectors
+    recipvecs      :: Mat3                                 # Reciprocal lattice vectors (conventional)
     positions      :: Vector{Vec3}                         # Positions in fractional coords
     types          :: Vector{String}                       # Types
     classes        :: Vector{Int}                          # Class indices
@@ -94,18 +95,6 @@ Volume of the crystal unit cell.
 """
 cell_volume(cryst::Crystal) = abs(det(cryst.latvecs))
 
-"""
-    reciprocal_vectors(cryst::Crystal)
-
-A 3×3 matrix with columns that are the reciprocal lattice vectors of the
-crystallographic unit cell. Defined such that:
-
-```julia
-@assert reciprocal_vectors(cryst)' * cryst.latvecs ≈ 2π*I
-```
-"""
-reciprocal_vectors(cryst::Crystal) = 2π*Mat3(inv(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.
@@ -204,6 +193,7 @@ function crystal_from_inferred_symmetry(latvecs::Mat3, positions::Vector{Vec3},
         end
     end
 
+    recipvecs = 2π*Mat3(inv(latvecs)')
     positions = wrap_to_unit_cell.(positions; symprec)
 
     cell = Spglib.Cell(latvecs, positions, types)
@@ -231,7 +221,7 @@ function crystal_from_inferred_symmetry(latvecs::Mat3, positions::Vector{Vec3},
 
     sitesyms = SiteSymmetry.(d.site_symmetry_symbols, multiplicities, d.wyckoffs)
 
-    ret = Crystal(latvecs, d.primitive_lattice, positions, types, classes, sitesyms, symops, spacegroup, symprec)
+    ret = Crystal(latvecs, d.primitive_lattice, recipvecs, positions, types, classes, sitesyms, symops, spacegroup, symprec)
     validate(ret)
     return ret
 end
@@ -393,7 +383,8 @@ function crystal_from_symops(latvecs::Mat3, positions::Vector{Vec3}, types::Vect
         inferred
     else
         prim_latvecs = latvecs
-        Crystal(latvecs, prim_latvecs, all_positions, all_types, classes, nothing, symops, spacegroup, symprec)
+        recipvecs = 2π*Mat3(inv(latvecs)')
+        Crystal(latvecs, prim_latvecs, recipvecs, all_positions, all_types, classes, nothing, symops, spacegroup, symprec)
     end
     sort_sites!(ret)
     validate(ret)
@@ -411,6 +402,9 @@ function reshape_crystal(cryst::Crystal, new_cell_size::Mat3)
     # Lattice vectors of the new unit cell in global coordinates
     new_latvecs = cryst.latvecs * new_cell_size
 
+    # Reciprocal lattice vectors of the new unit cell
+    new_recipvecs = 2π * Mat3(inv(new_latvecs)')
+
     # These don't change because both are in global coordinates
     prim_latvecs = cryst.prim_latvecs
 
@@ -463,8 +457,8 @@ function reshape_crystal(cryst::Crystal, new_cell_size::Mat3)
     # lost with the resizing procedure.
     new_symops = SymOp[]
 
-    return Crystal(new_latvecs, prim_latvecs, new_positions, new_types, new_classes, new_sitesyms,
-                new_symops, cryst.spacegroup, new_symprec)
+    return Crystal(new_latvecs, prim_latvecs, new_recipvecs, new_positions, new_types, new_classes,
+                   new_sitesyms, new_symops, cryst.spacegroup, new_symprec)
 end
 
 

src/Symmetry/SymmetryAnalysis.jl

@@ -159,12 +159,10 @@ function all_bonds_for_atom(cryst::Crystal, i::Int, max_dist; min_dist=0.0)
     max_dist += 4 * cryst.symprec * ℓ
     min_dist -= 4 * cryst.symprec * ℓ
 
-    recip_vecs = reciprocal_vectors(cryst)
-
     # box_lengths[i] represents the perpendicular distance between two parallel
     # boundary planes spanned by lattice vectors a_j and a_k (where indices j
     # and k differ from i)
-    box_lengths = [a⋅b/norm(b) for (a,b) = zip(eachcol(cryst.latvecs), eachcol(recip_vecs))]
+    box_lengths = [a⋅b/norm(b) for (a,b) = zip(eachcol(cryst.latvecs), eachcol(cryst.recipvecs))]
     n_max = round.(Int, max_dist ./ box_lengths, RoundUp)
 
     bonds = Bond[]

src/System/Ewald.jl

@@ -42,12 +42,12 @@ function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) :: Arra
 
     # Superlattice vectors that describe periodicity of system and their inverse
     latscale = diagm(Vec3(latsize))
-    supervecs = collect(eachcol(cryst.latvecs * latscale))
-    recipvecs = collect(eachcol(latscale \ reciprocal_vectors(cryst)))
+    slatvecs = collect(eachcol(cryst.latvecs * latscale))
+    srecipvecs = collect(eachcol(latscale \ cryst.recipvecs))
 
     # Precalculate constants
     I₃ = Mat3(I)
-    V = (supervecs[1] × supervecs[2]) ⋅ supervecs[3]
+    V = (slatvecs[1] × slatvecs[2]) ⋅ slatvecs[3]
     L = cbrt(V)
     # Roughly balances the real and Fourier space costs
     σ = L/3
@@ -57,10 +57,10 @@ function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) :: Arra
     rmax = 6√2 * σ
     kmax = 6√2 / σ
 
-    nmax = map(supervecs, recipvecs) do a, b
+    nmax = map(slatvecs, srecipvecs) do a, b
         round(Int, rmax / (a⋅normalize(b)) + 1e-6) + 1
     end
-    mmax = map(supervecs, recipvecs) do a, b
+    mmax = map(slatvecs, srecipvecs) do a, b
         round(Int, kmax / (b⋅normalize(a)) + 1e-6)
     end
 
@@ -76,7 +76,7 @@ function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) :: Arra
         ## Real space part
         for n1 = -nmax[1]:nmax[1], n2 = -nmax[2]:nmax[2], n3 = -nmax[3]:nmax[3]
             n = Vec3(n1, n2, n3)
-            rvec = Δr + n' * supervecs
+            rvec = Δr + n' * slatvecs
             r² = rvec⋅rvec
             if 0 < r² <= rmax*rmax
                 r = √r²
@@ -91,7 +91,7 @@ function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) :: Arra
         #####################################################
         ## Fourier space part
         for m1 = -mmax[1]:mmax[1], m2 = -mmax[2]:mmax[2], m3 = -mmax[3]:mmax[3]
-            k = Vec3(m1, m2, m3)' * recipvecs
+            k = Vec3(m1, m2, m3)' * srecipvecs
             k² = k⋅k
             if 0 < k² <= kmax*kmax
                 # Replace exp(-ikr) -> cos(kr). It's valid to drop the imaginary

test/test_lswt.jl

@@ -33,7 +33,7 @@ end
     randomize_spins!(sys)
     @test minimize_energy!(sys) > 0
 
-    k = reciprocal_vectors(cryst) * rand(Float64, 3)
+    k = cryst.recipvecs * rand(Float64, 3)
     swt = SpinWaveTheory(sys)
     ωk_num = dispersion(swt, [k])[1, :]
 
@@ -89,7 +89,7 @@ end
         set_dipole!(sys, (-1, -1, 1), position_to_site(sys, (1/2, 0, 1/2)))
         set_dipole!(sys, (-1, 1, -1), position_to_site(sys, (0, 1/2, 1/2)))
         swt = SpinWaveTheory(sys)
-        k = reciprocal_vectors(fcc) * [0.8, 0.6, 0.1]
+        k = fcc.recipvecs * [0.8, 0.6, 0.1]
         _, Sαβs =  Sunny.dssf(swt, [k])
 
         sunny_trace = [real(tr(Sαβs[1,a])) for a in axes(Sαβs)[2]]
@@ -149,7 +149,7 @@ end
         latvecs = lattice_vectors(a, a, 10a, 90, 90, 90)
         positions = [[0, 0, 0]]
         cryst = Crystal(latvecs, positions)
-        q = reciprocal_vectors(cryst) * [0.12, 0.23, 0.34]
+        q = cryst.recipvecs * [0.12, 0.23, 0.34]
 
         dims = (2, 2, 1)
         sys = System(cryst, dims, [SpinInfo(1; S, g=1)], :dipole; units=Units.theory)