diff --git a/src/EntangledUnits/EntangledReshaping.jl b/src/EntangledUnits/EntangledReshaping.jl new file mode 100644 index 000000000..d4337f6da --- /dev/null +++ b/src/EntangledUnits/EntangledReshaping.jl @@ -0,0 +1,89 @@ +#TODO: Test this rigorously -- this is key to reshaping EntangledSystems and +# making EntangledSpinWaveTheorys. + +# An entangled system can be specified with a list of +# tuples, e.g. [(1,2), (3,4)], which will group the original sites 1 and 2 into +# one unit and 3 and 4 into another. Given an EntangledSystem constructed from a +# sys_origin and and a reshaped sys_origin, this function returns a list of +# tuples for specifying the corresponding entanglement of the reshaped system. +function units_for_reshaped_system(reshaped_sys_origin, esys) + (; sys_origin) = esys + units = original_unit_spec(esys) + new_crystal = reshaped_sys_origin.crystal + new_atoms = collect(1:natoms(new_crystal)) + new_units = [] + + # Take a list of all the new atoms in the reshaped system. Pick the first. + # Map it back to the original system to determine what unit it belongs to. + # Then map all members of the unit forward to define the unit in terms of + # the atoms of the reshaped system. Remove these atoms from the list of + # sites left to be "entangled" and repeat until list of new atoms is + # exhausted. + while length(new_atoms) > 0 + # Pick any site from list of new sites + new_atom = new_atoms[1] + new_site = CartesianIndex(1, 1, 1, new_atom) # Only need to define for a single unit cell, may as well be the first + new_position = position(reshaped_sys_origin, new_site) + + # Find corresponding original atom number. + site = position_to_site(sys_origin, position(reshaped_sys_origin, new_site)) + original_atom = site[4] + position_of_corresponding_atom = position(sys_origin, (1, 1, 1, original_atom)) + offset = new_position - position_of_corresponding_atom + + # Find the unit to which this original atom belongs. + unit = units[findfirst(unit -> original_atom in unit, units)] + + # Find positions of all atoms in the unit, find corresponding sites in reshape system, and define unit for reshaped system. + unit_positions = [position(sys_origin, CartesianIndex(1, 1, 1, atom)) for atom in unit] + new_unit_sites = [position_to_site(reshaped_sys_origin, position + offset) for position in unit_positions] + new_unit = Int64[] + for new_site in new_unit_sites + i, j, k, a = new_site.I + if !(i == j == k == 1) + error("Specified reshaping incompatible with specified entangled units. (Unit split between crystalographic unit cells.)") + end + push!(new_unit, a) + end + push!(new_units, Tuple(new_unit)) + + # Delete members of newly defined unit from list of all atoms in the + # reshaped system. + idcs = findall(atom -> atom in new_unit, new_atoms) + deleteat!(new_atoms, idcs) + end + + return new_units +end + +function reshape_supercell(esys::EntangledSystem, shape) + (; sys, sys_origin) = esys + + # Reshape the origin System. + sys_origin_new = reshape_supercell(sys_origin, shape) + + # Reshape the the underlying "entangled" System. + units_new = units_for_reshaped_system(sys_origin_new, esys) + _, contraction_info = contract_crystal(sys_origin_new.crystal, units_new) + sys_new = reshape_supercell(sys, shape) + + # Construct dipole operator field for reshaped EntangledSystem + dipole_operators_origin = all_dipole_observables(sys_origin_new; apply_g=false) + (; observables, source_idcs) = observables_to_product_space(dipole_operators_origin, sys_origin_new, contraction_info) + + return EntangledSystem(sys_new, sys_origin_new, contraction_info, observables, source_idcs) +end + +function repeat_periodically(esys, counts) + (; sys, sys_origin, contraction_info) = esys + + # Repeat both entangled and original system periodically + sys_new = repeat_periodically(sys, counts) + sys_origin_new = repeat_periodically(sys_origin, counts) + + # Construct dipole operator field for reshaped EntangledSystem + dipole_operators_origin = all_dipole_observables(sys_origin_new; apply_g=false) + (; observables, source_idcs) = observables_to_product_space(dipole_operators_origin, sys_origin_new, contraction_info) + + return EntangledSystem(sys_new, sys_origin_new, contraction_info, observables, source_idcs) +end \ No newline at end of file diff --git a/src/EntangledUnits/EntangledSampledCorrelations.jl b/src/EntangledUnits/EntangledSampledCorrelations.jl new file mode 100644 index 000000000..cd3628512 --- /dev/null +++ b/src/EntangledUnits/EntangledSampledCorrelations.jl @@ -0,0 +1,170 @@ +struct EntangledSampledCorrelations + sc::SampledCorrelations # Parent SampledCorrelations + esys::EntangledSystem # Probably don't need to carry around the whole thing -- defeats spirit of original design for SC +end + +struct EntangledSampledCorrelationsStatic + sc::SampledCorrelationsStatic # Parent SampledCorrelations + esys::EntangledSystem # Probably don't need to carry around the whole thing -- defeats spirit of original design for SC +end + +function Base.setproperty!(sc::T, sym::Symbol, val) where {T<:Union{EntangledSampledCorrelations, EntangledSampledCorrelationsStatic}} + sc = sc.sc + if sym == :measure + @assert sc.measure.observables ≈ val.observables "New MeasureSpec must contain identical observables." + @assert all(x -> x == 1, sc.measure.corr_pairs .== val.corr_pairs) "New MeasureSpec must contain identical correlation pairs." + setfield!(sc, :measure, val) + else + setfield!(sc, sym, val) + end +end + +function Base.show(io::IO, ::EntangledSampledCorrelations) + print(io, "EntangledSampledCorrelations") +end + +function Base.show(io::IO, ::EntangledSampledCorrelationsStatic) + print(io, "EntangledSampledCorrelationsStatic") +end + +function Base.show(io::IO, ::MIME"text/plain", esc::EntangledSampledCorrelations) + (; crystal, sys_dims, nsamples) = esc.sc + printstyled(io, "EntangledSampledCorrelations"; bold=true, color=:underline) + println(io," ($(Base.format_bytes(Base.summarysize(esc))))") + print(io,"[") + printstyled(io,"S(q,ω)"; bold=true) + print(io," | nω = $(round(Int, size(esc.sc.data)[7]/2 + 1)), Δω = $(round(esc.sc.Δω, digits=4))") + println(io," | $nsamples $(nsamples > 1 ? "samples" : "sample")]") + println(io,"Lattice: $sys_dims × $(natoms(crystal))") +end + +function Base.show(io::IO, ::MIME"text/plain", esc::EntangledSampledCorrelationsStatic) + (; crystal, sys_dims, nsamples) = esc.sc.parent + printstyled(io, "SampledCorrelationsStatic"; bold=true, color=:underline) + println(io," ($(Base.format_bytes(Base.summarysize(esc))))") + print(io,"[") + printstyled(io,"S(q)"; bold=true) + println(io," | $nsamples $(nsamples > 1 ? "samples" : "sample")]") + println(io,"Lattice: $sys_dims × $(natoms(crystal))") +end + +function Base.setproperty!(esc::EntangledSampledCorrelations, sym::Symbol, val) + if sym == :measure + measure = val + (; observables) = observables_to_product_space(measure.observables, esc.esys.sys_origin, esc.esys.contraction_info) + @assert esc.sc.measure.observables ≈ observables "New MeasureSpec must contain identical observables." + @assert all(x -> x == 1, esc.sc.measure.corr_pairs .== measure.corr_pairs) "New MeasureSpec must contain identical correlation pairs." + setfield!(esc.sc, :measure, val) # Sets new combiner + else + setfield!(sc, sym, val) + end +end + +Base.setproperty!(esc::EntangledSampledCorrelationsStatic, sym::Symbol, val) = setproperty!(esc.parent, sym, val) + + +# Take observables specified in terms or original system and transform them into +# a field of observables in the tensor product space, together with mapping +# information for populating the expectation values of these operators with +# respect to the entangled system. +function observables_to_product_space(observables, sys_origin, contraction_info) + Ns_per_unit = Ns_in_units(sys_origin, contraction_info) + Ns_all = [prod(Ns) for Ns in Ns_per_unit] + N = Ns_all[1] + @assert all(==(N), Ns_all) "All entangled units must have the same dimension Hilbert space." + + observables_new = fill(zeros(ComplexF64, N, N), size(observables)) + positions = zeros(Vec3, size(observables)[2:end]) + source_idcs = zeros(Int64, size(observables)[2:end]) + + for site in eachsite(sys_origin) + atom = site.I[4] + positions[site] = sys_origin.crystal.positions[atom] + source_idcs[site] = contraction_info.forward[atom][1] + for μ in axes(observables, 1) + obs = observables[μ, site] + unit, unitsub = contraction_info.forward[atom] + Ns = Ns_per_unit[unit] + observables_new[μ, site] = local_op_to_product_space(obs, unitsub, Ns) + end + end + observables = observables_new + + return (; observables, positions, source_idcs) +end + + +# Configures an ordinary SampledCorrelations for an entangled system. The +# measure is assumed to correspond to the sites of the original "unentangled" +# system. +function SampledCorrelations(esys::EntangledSystem; measure, energies, dt, calculate_errors=false) + # Convert observables on different sites to "multiposition" observables in + # tensor product spaces. With entangled units, the position of an operator + # cannot be uniquely determined from the `atom` index of a `Site`. Instead, + # the position is recorded independently, and the index of the relevant + # coherent state (which may now be used for operators corresponding to + # multiple positions) is recorded in `source_idcs`. + (; observables, positions, source_idcs) = observables_to_product_space(measure.observables, esys.sys_origin, esys.contraction_info) + + # Make a sampled correlations for the esys. + sc = SampledCorrelations(esys.sys; measure, energies, dt, calculate_errors, positions) + + # Replace relevant fields or the resulting SampledCorrelations. Note use of + # undocumented `positions` keyword. This can be eliminated if positions are + # migrated into the MeasureSpec. + crystal = esys.sys_origin.crystal + origin_crystal = orig_crystal(esys.sys_origin) + sc_new = SampledCorrelations(sc.data, sc.M, crystal, origin_crystal, sc.Δω, measure, observables, positions, source_idcs, measure.corr_pairs, + sc.measperiod, sc.dt, sc.nsamples, sc.samplebuf, sc.corrbuf, sc.space_fft!, sc.time_fft!, sc.corr_fft!, sc.corr_ifft!) + + return EntangledSampledCorrelations(sc_new, esys) +end + +function SampledCorrelationsStatic(esys::EntangledSystem; measure, calculate_errors=false) + (; observables, positions, source_idcs) = observables_to_product_space(measure.observables, esys.sys_origin, esys.contraction_info) + sc = SampledCorrelations(esys.sys; measure, energies=nothing, dt=NaN, calculate_errors, positions) + + # Replace relevant fields + crystal = esys.sys_origin.crystal + origin_crystal = orig_crystal(esys.sys_origin) + parent = SampledCorrelations(sc.data, sc.M, crystal, origin_crystal, sc.Δω, measure, observables, positions, source_idcs, measure.corr_pairs, + sc.measperiod, sc.dt, sc.nsamples, sc.samplebuf, sc.corrbuf, sc.space_fft!, sc.time_fft!, sc.corr_fft!, sc.corr_ifft!) + + ssc = SampledCorrelationsStatic(parent) + return EntangledSampledCorrelationsStatic(ssc, esys) +end + + +function add_sample!(esc::EntangledSampledCorrelations, esys::EntangledSystem; window=:cosine) + new_sample!(esc.sc, esys.sys) + accum_sample!(esc.sc; window) +end + +function add_sample!(esc::EntangledSampledCorrelationsStatic, esys::EntangledSystem; window=:cosine) + add_sample!(esc.sc, esys.sys) +end + +available_energies(esc::EntangledSampledCorrelations) = available_energies(esc.sc) +available_wave_vectors(esc::EntangledSampledCorrelations) = available_wave_vectors(esc.sc) + +function clone_correlations(esc::EntangledSampledCorrelations; kwargs...) + sc = clone_correlations(esc.sc) + EntangledSampledCorrelations(sc, esc.esys) +end + +function merge_correlations(escs::Vector{EntangledSampledCorrelations}; kwargs...) + sc_merged = merge_correlations([esc.sc for esc in escs]) + EntangledSampledCorrelations(sc_merged, escs[1].esys) +end + +function intensities(esc::EntangledSampledCorrelations, qpts; kwargs...) + intensities(esc.sc, qpts; kwargs...) +end + +function intensities_static(esc::EntangledSampledCorrelations, qpts; kwargs...) + intensities_static(esc.sc, qpts; kwargs...) +end + +function intensities_static(esc::EntangledSampledCorrelationsStatic, qpts; kwargs...) + intensities_static(esc.sc, qpts; kwargs...) +end \ No newline at end of file diff --git a/src/EntangledUnits/EntangledSpinWaveTheory.jl b/src/EntangledUnits/EntangledSpinWaveTheory.jl new file mode 100644 index 000000000..0c5b134d1 --- /dev/null +++ b/src/EntangledUnits/EntangledSpinWaveTheory.jl @@ -0,0 +1,335 @@ +struct SWTDataEntangled + local_unitaries :: Vector{Matrix{ComplexF64}} # Aligns to quantization axis on each site + observables_localized :: Array{HermitianC64, 3} # Observables in local frame for each subsite (for intensity calcs) + observable_buf :: Array{ComplexF64, 2} # Buffer for use while constructing boson rep of observables +end + +struct EntangledSpinWaveTheory <: AbstractSpinWaveTheory + sys :: System + data :: SWTDataEntangled + measure :: MeasureSpec + regularization :: Float64 + + crystal_origin :: Crystal + contraction_info :: CrystalContractionInfo + Ns_unit :: Vector{Vector{Int64}} +end + + +function SpinWaveTheory(esys::EntangledSystem; measure, regularization=1e-8) + (; sys, sys_origin) = esys + crystal_origin = orig_crystal(sys_origin) + if !isnothing(sys.ewald) + error("EntangledSpinWaveTheory does not yet support long-range dipole-dipole interactions.") + end + + measure = @something measure empty_measurespec(sys) + if length(eachsite(sys_origin)) != prod(size(measure.observables)[2:5]) + error("Size mismatch. Check that measure is built using consistent system.") + end + + # Reshape (flatten) both sys and sys_origin in corresponding ways + sys, sys_origin = map([sys, sys_origin]) do sys + new_shape = cell_shape(sys) * diagm(Vec3(sys.dims)) + new_cryst = reshape_crystal(orig_crystal(sys), new_shape) + + # Sort crystal positions so that their order matches sites in sys. Quadratic + # scaling in system size. + global_positions = global_position.(Ref(sys), vec(eachsite(sys))) + p = map(new_cryst.positions) do r + pos = new_cryst.latvecs * r + findfirst(global_positions) do refpos + isapprox(pos, refpos, atol=new_cryst.symprec) + end + end + @assert allunique(p) + permute_sites!(new_cryst, p) + + # Create a new system with dims (1,1,1). A clone happens in all cases. + return reshape_supercell_aux(sys, new_cryst, (1,1,1)) + end + + # Construct the new contraction_info (i.e. reconstruct maps between + # entangled and unentangled systems) + new_units = units_for_reshaped_system(sys_origin, esys) + _, new_contraction_info = contract_crystal(sys_origin.crystal, new_units) + new_Ns_unit = Ns_in_units(sys_origin, new_contraction_info) + + data = swt_data_entangled(sys, sys_origin, new_Ns_unit, new_contraction_info, measure) + + return EntangledSpinWaveTheory(sys, data, measure, regularization, crystal_origin, new_contraction_info, new_Ns_unit) +end + +function Base.show(io::IO, ::MIME"text/plain", swt::EntangledSpinWaveTheory) + printstyled(io, "EntangledSpinWaveTheory\n"; bold=true, color=:underline) + println(io, "Entangled units in magnetic supercell: $(natoms(swt.sys.crystal))") + # Add observable info? +end + +nbands(swt::EntangledSpinWaveTheory) = (swt.sys.Ns[1]-1) * natoms(swt.sys.crystal) + +function to_reshaped_rlu(esys::EntangledSystem, q) + (; sys, sys_origin) = esys + return sys.crystal.recipvecs \ (orig_crystal(sys_origin).recipvecs * q) +end + +# obs are observables _given in terms of `sys_original`_ +function swt_data_entangled(sys, sys_origin, Ns_unit, contraction_info, measure) + + # Calculate transformation matrices into local reference frames + N = sys.Ns[1] # Assume uniform contraction for now + nunits = length(eachsite(sys)) + natoms = length(eachsite(sys_origin)) + nobs = size(measure.observables, 1) + observables = reshape(measure.observables, nobs, natoms) + + # Check to make sure all "units" are the same -- this can be generalized later. + atoms_per_unit = [length(info) for info in contraction_info.inverse] + @assert allequal(atoms_per_unit) "All units must have the same number of interior sites" + atoms_per_unit = atoms_per_unit[1] + + Ns_contracted = map(Ns -> prod(Ns), Ns_unit) + @assert allequal(Ns_contracted) "All units must have the same dimension local Hilbert space" + @assert Ns_contracted[1] == N "Unit dimension inconsistent with system" # Sanity check. This should never happen. + + # Preallocate buffers for local unitaries and observables. + local_unitaries = Vector{Matrix{ComplexF64}}(undef, nunits) + observables_localized = Array{HermitianC64}(undef, atoms_per_unit, nobs, nunits) + observable_buf = zeros(ComplexF64, N, N) + + # Find local unitaries (for units) + for i in 1:nunits + # Create unitary that rotates [0, ..., 0, 1] into ground state direction + # Z that defines quantization axis + Z = sys.coherents[i] + U = hcat(nullspace(Z'), Z) + local_unitaries[i] = U + + # Rotate observables into local reference frames + for μ in 1:num_observables(measure) + A = observables[μ, i] + for ui in 1:atoms_per_unit + A_prod = local_op_to_product_space(A, ui, Ns_unit[i]) + observables_localized[ui, μ, i] = Hermitian(U' * convert(Matrix, A_prod) * U) + end + end + end + + # Rotate all observables in each unit into local reference frame. + for unit in 1:nunits + U = local_unitaries[unit] + + # Rotate interactions into local reference frames + int = sys.interactions_union[unit] + + # Rotate onsite anisotropy (not that, for entangled units, onsite already includes Zeeman) + int.onsite = Hermitian(U' * int.onsite * U) + + # Transform pair couplings into tensor decomposition and rotate. + pair_new = PairCoupling[] + for pc in int.pair + # Convert PairCoupling to a purely general (tensor decomposed) interaction. + pc_general = as_general_pair_coupling(pc, sys) + + # Rotate tensor decomposition into local frame. + bond = pc.bond + @assert bond.i == unit + U′ = local_unitaries[bond.j] + pc_rotated = rotate_general_coupling_into_local_frame(pc_general, U, U′) + + push!(pair_new, pc_rotated) + end + int.pair = pair_new + end + + return SWTDataEntangled( + local_unitaries, + observables_localized, + observable_buf, + ) +end + + +function intensities_bands(swt::EntangledSpinWaveTheory, qpts; kT=0) + (; sys, contraction_info, data, measure, crystal_origin) = swt + isempty(measure.observables) && error("No observables! Construct SpinWaveTheory with a `measure` argument.") + + qpts = convert(AbstractQPoints, qpts) + cryst = orig_crystal(sys) + + # Number of atoms in magnetic cell + @assert sys.dims == (1,1,1) + nunits = length(eachsite(sys)) + # Number of chemical cells in magnetic cell + # Ncells = Na / natoms(cryst) # TODO: Pass information about natoms in unreshaped, uncontracted system + Ncells = 1 + # Number of quasiparticle modes + L = nbands(swt) + # Number of wavevectors + Nq = length(qpts.qs) + + # Preallocation + T = zeros(ComplexF64, 2L, 2L) + H = zeros(ComplexF64, 2L, 2L) + Avec_pref = zeros(ComplexF64, nunits) + disp = zeros(Float64, L, Nq) + intensity = zeros(eltype(measure), L, Nq) + + # Temporary storage for pair correlations + Ncorr = length(measure.corr_pairs) + corrbuf = zeros(ComplexF64, Ncorr) + + Nobs = size(measure.observables, 1) + + for (iq, q) in enumerate(qpts.qs) + q_global = cryst.recipvecs * q + q_reshaped = cryst.recipvecs \ (crystal_origin.recipvecs * q) + view(disp, :, iq) .= view(excitations!(T, H, swt, q), 1:L) + + for i in 1:nunits + r_global = global_position(sys, (1,1,1,i)) + Avec_pref[i] = exp(- im * dot(q_global, r_global)) + end + + Avec = zeros(ComplexF64, Nobs) + + # Fill `intensity` array + O = data.observable_buf + for band = 1:L + fill!(Avec, 0) + N = sys.Ns[1] + t = reshape(view(T, :, band), N-1, nunits, 2) + for i in 1:nunits + inverse_infos = contraction_info.inverse[i] + for μ in 1:Nobs + fill!(O, 0) + for (subsite, inverse_info) in enumerate(inverse_infos) + site_original, offset = inverse_info.site, inverse_info.offset + ff = get_swt_formfactor(measure, 1, site_original) + prefactor = exp(-2π*im*(q_reshaped ⋅ offset)) * compute_form_factor(ff, norm2(q_global)) + O += prefactor * data.observables_localized[subsite, μ, i] + end + for α in 1:N-1 + Avec[μ] += Avec_pref[i] * (O[α, N] * t[α, i, 2] + O[N, α] * t[α, i, 1]) + end + end + end + + map!(corrbuf, measure.corr_pairs) do (α, β) + Avec[α] * conj(Avec[β]) / Ncells + end + intensity[band, iq] = thermal_prefactor(disp[band]; kT) * measure.combiner(q_global, corrbuf) + end + end + + disp = reshape(disp, L, size(qpts.qs)...) + intensity = reshape(intensity, L, size(qpts.qs)...) + return BandIntensities(cryst, qpts, disp, intensity) +end + +################################################################################ +# To be simplified, but requires editing existing code +################################################################################ +# No changes +function dynamical_matrix!(H, swt::EntangledSpinWaveTheory, q_reshaped) + if swt.sys.mode == :SUN + swt_hamiltonian_SUN!(H, swt, q_reshaped) + else + @assert swt.sys.mode in (:dipole, :dipole_large_S) + swt_hamiltonian_dipole!(H, swt, q_reshaped) + end +end + +# No changes +function excitations!(T, tmp, swt::EntangledSpinWaveTheory, q) + L = nbands(swt) + size(T) == size(tmp) == (2L, 2L) || error("Arguments T and tmp must be $(2L)×$(2L) matrices") + + q_reshaped = to_reshaped_rlu(swt.sys, q) + dynamical_matrix!(tmp, swt, q_reshaped) + + try + return bogoliubov!(T, tmp) + catch _ + error("Instability at wavevector q = $q") + end +end + +# Only change is no Ewald +function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::EntangledSpinWaveTheory, q_reshaped::Vec3) + (; sys) = swt + + N = sys.Ns[1] + Na = natoms(sys.crystal) + L = (N-1) * Na + + # Clear the Hamiltonian + @assert size(H) == (2L, 2L) + H .= 0 + blockdims = (N-1, Na, N-1, Na) + H11 = reshape(view(H, 1:L, 1:L), blockdims) + H12 = reshape(view(H, 1:L, L+1:2L), blockdims) + H21 = reshape(view(H, L+1:2L, 1:L), blockdims) + H22 = reshape(view(H, L+1:2L, L+1:2L), blockdims) + + for (i, int) in enumerate(sys.interactions_union) + + # Onsite coupling, including Zeeman. Note that op has already been + # transformed according to the local frame of sublattice i. + op = int.onsite + for m in 1:N-1 + for n in 1:N-1 + c = 0.5 * (op[m, n] - δ(m, n) * op[N, N]) + H11[m, i, n, i] += c + H22[n, i, m, i] += c + end + end + + for coupling in int.pair + (; isculled, bond) = coupling + isculled && break + (; i, j) = bond + phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping + + # Set "general" pair interactions of the form Aᵢ⊗Bⱼ. Note that Aᵢ + # and Bᵢ have already been transformed according to the local frames + # of sublattice i and j, respectively. + for (Ai, Bj) in coupling.general.data + for m in 1:N-1, n in 1:N-1 + c = 0.5 * (Ai[m,n] - δ(m,n)*Ai[N,N]) * (Bj[N,N]) + H11[m, i, n, i] += c + H22[n, i, m, i] += c + + c = 0.5 * Ai[N,N] * (Bj[m,n] - δ(m,n)*Bj[N,N]) + H11[m, j, n, j] += c + H22[n, j, m, j] += c + + c = 0.5 * Ai[m,N] * Bj[N,n] + H11[m, i, n, j] += c * phase + H22[n, j, m, i] += c * conj(phase) + + c = 0.5 * Ai[N,m] * Bj[n,N] + H11[n, j, m, i] += c * conj(phase) + H22[m, i, n, j] += c * phase + + c = 0.5 * Ai[m,N] * Bj[n,N] + H12[m, i, n, j] += c * phase + H12[n, j, m, i] += c * conj(phase) + H21[n, j, m, i] += conj(c) * conj(phase) + H21[m, i, n, j] += conj(c) * phase + end + end + end + end + + # H must be hermitian up to round-off errors + @assert diffnorm2(H, H') < 1e-12 + + # Make H exactly hermitian + hermitianpart!(H) + + # Add small constant shift for positive-definiteness + for i in 1:2L + H[i,i] += swt.regularization + end +end \ No newline at end of file diff --git a/src/EntangledUnits/EntangledUnits.jl b/src/EntangledUnits/EntangledUnits.jl new file mode 100644 index 000000000..0ebf76b12 --- /dev/null +++ b/src/EntangledUnits/EntangledUnits.jl @@ -0,0 +1,340 @@ +################################################################################ +# Crystal contraction logic +################################################################################ + +# Takes a crystal and a list of integer tuples. The tuples indicate which sites +# in the original crystal are to be grouped, i.e., contracted into a single site +# of a new crystal. +function contract_crystal(crystal, units) + + # Determine which sites are entangled and which are not. + unentangled_sites = 1:natoms(crystal) + entangled_sites = Int64[] + for unit in units, site in unit + push!(entangled_sites, site) + unentangled_sites = filter(!=(site), unentangled_sites) + end + + # Check that unit definitions are valid. + @assert length(entangled_sites) + length(unentangled_sites) == natoms(crystal) "Invalid entangled unit specification." + @assert isempty(intersect(entangled_sites, unentangled_sites)) "Invalid entangled unit specification." # Sanity check. Should be true by construction. Remove later + + # Prepare mapping dictionaries and initialize iteration variable. + forward_map = Dict{Int64, Tuple{Int64, Int64}}() + inverse_map = Dict{Tuple{Int64, Int64}, InverseData}() + new_site_current = 1 + + # Add sites that are *not* encapsulated in any unit as individual sites in + # the new crystal with no associated displacement. + new_positions = [] + for site in unentangled_sites + push!(new_positions, crystal.positions[site]) + new_pair = (new_site_current, 1) + forward_map[site] = new_pair + inverse_map[new_pair] = InverseData(site, Vec3(0, 0, 0)) + new_site_current += 1 + end + + # Assign entangled units to single site in new crystal and record mapping + # information. + for unit in units + # Find new position by averaging location of entangled positions. + old_positions = [crystal.positions[i] for i in unit] + new_position = sum(old_positions) / length(old_positions) + push!(new_positions, new_position) + + # Record forward and inverse mapping information, including + # the displacement data from the new unit position. + for (j, site) in enumerate(unit) + new_pair = (new_site_current, j) + offset = crystal.positions[site] - new_position + + forward_map[site] = new_pair + inverse_map[new_pair] = InverseData(site, offset) + end + + new_site_current += 1 + end + + nsites_new = new_site_current - 1 + + # Sorted list version of forward map + forward_keys = collect(keys(forward_map)) + forward_vals = collect(values(forward_map)) + idcs = sortperm(forward_keys) + forward_list = [forward_vals[i] for i in idcs] + + # Sorted list version of inverse map + inverse_keys = collect(keys(inverse_map)) + inverse_vals = collect(values(inverse_map)) + idcs = sortperm(inverse_keys) + inverse_keys = inverse_keys[idcs] + inverse_vals = inverse_vals[idcs] + + inverse_list = [InverseData[] for _ in 1:nsites_new] + for (n, key) in enumerate(inverse_keys) + new_site, _ = key # `key` is a tuple (global_site, local_site_in_unit) + push!(inverse_list[new_site], inverse_vals[n]) + end + + # Generate a new contracted crystal and information to invert contraction. + # Space group must be set to 1 to allow "invalid" anisotropies -- these are + # generated by PairCouplings that become onsite couplings. NB: Add option to + # override to unconventional cell warnings and errors? (Probably yes) + new_crystal = Crystal(crystal.latvecs, new_positions, 1) + contraction_info = CrystalContractionInfo(forward_list, inverse_list) + + return new_crystal, contraction_info +end + + +# Reconstruct original crystal from contracted Crystal and a CrystalContractionInfo +function expand_crystal(contracted_crystal, contraction_info) + (; forward, inverse) = contraction_info + contracted_positions = contracted_crystal.positions + nsites_expanded = length(forward) + expanded_positions = [Vec3(0, 0, 0) for _ in 1:nsites_expanded] + for (contracted_idx, original_site_data) in enumerate(inverse) + for (; site, offset) in original_site_data + expanded_positions[site] = contracted_positions[contracted_idx] + offset + end + end + Crystal(contracted_crystal.latvecs, expanded_positions) +end + + +# Returns a list of length equal to the number of "units" in the a contracted +# crystal. Each list element is itself a list of integers, each of which +# corresponds to the N of the corresponding site of the original system. The +# order is consistent with that given by the `inverse` field of a +# `CyrstalContractionInfo`. +function Ns_in_units(sys_original, contraction_info) + Ns = [Int64[] for _ in 1:length(contraction_info.inverse)] + for (n, contracted_sites) in enumerate(contraction_info.inverse) + for (; site) in contracted_sites + push!(Ns[n], sys_original.Ns[site]) + end + end + Ns +end + + +# Pull out original indices of sites in entangled unit +atoms_in_unit(contraction_info, i) = [inverse_data.site for inverse_data in contraction_info.inverse[i]] + +# Get the list of tuples specifying the units in terms of the uncontracted +# system. +function original_unit_spec(esys::EntangledSystem) + (; sys, contraction_info) = esys + return [Tuple(atoms_in_unit(contraction_info, unit)) for unit in axes(sys.dipoles, 4)] +end + + +# Interpreted in terms of the original crystal. +function bonds_in_unit(contraction_info, i) + sites = atoms_in_unit(contraction_info, i) + nsites = length(sites) + bonds = Bond[] + for i in 1:nsites, j in i+1:nsites + push!(bonds, Bond(i, j, [0, 0, 0])) + end + return bonds +end + + +# Checks whether a bond given in terms of the original crystal lies inside a +# single unit of the contracted system. +function bond_is_in_unit(bond::Bond, ci::CrystalContractionInfo) + (ci.forward[bond.i][1] == ci.forward[bond.j][1]) && (bond.n == [0, 0, 0]) +end + + +# Converts what was a pair coupling between two different sites in the original system into a single +# on-bond operator (an onsite operator in terms of the "units".) +function accum_pair_coupling_into_bond_operator_in_unit!(op, pc, sys, contracted_site, contraction_info) + (; bond, scalar, bilin, biquad, general, isculled) = pc + isculled && return + + Ns_all = Ns_in_units(sys, contraction_info) + Ns_unit = Ns_all[contracted_site] + I_unit = I(prod(Ns_unit)) + + # Collect local Hilbert space information + i, j = bond.i, bond.j + @assert contraction_info.forward[i][1] == contracted_site "Sanity check -- remove later" + @assert contraction_info.forward[j][1] == contracted_site "Sanity check -- remove later" + i_unit = contraction_info.forward[i][2] + j_unit = contraction_info.forward[j][2] + Ni = sys.Ns[1, 1, 1, i] + Nj = sys.Ns[1, 1, 1, j] + + # Add scalar part + op .+= scalar*I_unit + + # Add bilinear part + J = bilin isa Float64 ? bilin*I(3) : bilin + Si = [local_op_to_product_space(Sa, i_unit, Ns_unit) for Sa in spin_matrices((Ni-1)/2)] + Sj = [local_op_to_product_space(Sb, j_unit, Ns_unit) for Sb in spin_matrices((Nj-1)/2)] + op .+= Si' * J * Sj + + # Add biquadratic part + K = biquad isa Float64 ? diagm(biquad * Sunny.scalar_biquad_metric) : biquad + Oi = [local_op_to_product_space(Oa, i_unit, Ns_unit) for Oa in stevens_matrices_of_dim(2; N=Ni)] + Oj = [local_op_to_product_space(Ob, j_unit, Ns_unit) for Ob in stevens_matrices_of_dim(2; N=Nj)] + op .+= Oi' * K * Oj + + # Add general part + for (A, B) in general.data + op .+= local_op_to_product_space(A, i_unit, Ns_unit) * local_op_to_product_space(B, j_unit, Ns_unit) + end +end + +# Converts a pair coupling from the original system into a pair coupling between +# units in the new system. +function pair_coupling_into_bond_operator_between_units(pc, sys, contraction_info) + (; bond, scalar, bilin, biquad, general) = pc + (; i, j, n) = bond + unit1, unitsub1 = contraction_info.forward[i] + unit2, unitsub2 = contraction_info.forward[j] + + Ns_local = Ns_in_units(sys, contraction_info) + Ns_contracted = map(Ns -> prod(Ns), Ns_local) + Ns1 = Ns_local[unit1] + Ns2 = Ns_local[unit2] + N1 = sys.Ns[1, 1, 1, i] + N2 = sys.Ns[1, 1, 1, j] + N = Ns_contracted[unit1] * Ns_contracted[unit2] + newbond = Bond(unit1, unit2, n) + + bond_operator = zeros(ComplexF64, N, N) + + # Add scalar part + bond_operator += scalar*I(N) + + # Add bilinear part + J = bilin isa Float64 ? bilin*I(3) : bilin + Si = [kron(local_op_to_product_space(Sa, unitsub1, Ns1), I(Ns_contracted[unit1])) for Sa in spin_matrices((N1-1)/2)] + Sj = [kron(I(Ns_contracted[unit2]), local_op_to_product_space(Sa, unitsub2, Ns2)) for Sa in spin_matrices((N2-1)/2)] + bond_operator += Si' * J * Sj + + # Add biquadratic part + K = biquad isa Float64 ? diagm(biquad * Sunny.scalar_biquad_metric) : biquad + Oi = [kron(local_op_to_product_space(Oa, unitsub1, Ns1), I(Ns_contracted[unit1])) for Oa in stevens_matrices_of_dim(2; N=N1)] + Oj = [kron(I(Ns_contracted[unit2]), local_op_to_product_space(Ob, unitsub2, Ns2)) for Ob in stevens_matrices_of_dim(2; N=N2)] + bond_operator += Oi' * K * Oj + + # Add general part + for (A, B) in general.data + bond_operator += kron( local_op_to_product_space(A, unitsub1, Ns1), I(Ns_contracted[unit1]) ) * + kron( I(Ns_contracted[unit2]), local_op_to_product_space(B, unitsub2, Ns2) ) + end + + return (; newbond, bond_operator) +end + +function entangle_system(::System{0}, _) + error("Cannot contract a dipole system. Use :SUN mode.") +end + + +# More battle-tested approach. Iterates over over contracted system. Move to +# iteration scheme of original system moving forward for clarity and support of +# inhomogenous interactioninteractions +function entangle_system(sys::System{M}, units) where M + # Construct contracted crystal + contracted_crystal, contraction_info = contract_crystal(sys.crystal, units) + + # Determine Ns for local Hilbert spaces (all must be equal). (TODO: Determine if alternative behavior preferable in mixed case.) + Ns_unit = Ns_in_units(sys, contraction_info) + Ns_contracted = map(Ns -> prod(Ns), Ns_unit) + @assert allequal(Ns_contracted) "After contraction, the dimensions of the local Hilbert spaces on each generalized site must all be equal." + + # Construct empty contracted system + dims = size(sys.dipoles)[1:3] + spin_infos = [i => Moment(s=(N-1)/2, g=1.0) for (i, N) in enumerate(Ns_contracted)] # TODO: Decisions about g-factor + sys_entangled = System(contracted_crystal, spin_infos, :SUN; dims) + + # TODO: Extend to inhomogenous systems + # For each contracted site, scan original interactions and reconstruct as necessary. + new_pair_data = Tuple{Bond, Matrix{ComplexF64}}[] + for (contracted_site, N) in zip(1:natoms(contracted_crystal), Ns_contracted) + Ns = Ns_unit[contracted_site] + + ## Onsite portion of interaction + relevant_sites = atoms_in_unit(contraction_info, contracted_site) + unit_operator = zeros(ComplexF64, N, N) + + # Pair interactions that become within-unit interactions + original_interactions = sys.interactions_union[relevant_sites] + for (site, interaction) in zip(relevant_sites, original_interactions) + onsite_original = interaction.onsite + unit_index = contraction_info.forward[site][2] + unit_operator += local_op_to_product_space(onsite_original, unit_index, Ns) + end + + # Sort all PairCouplings in couplings that will be within a unit and couplings that will be between units + pcs_intra = PairCoupling[] + pcs_inter = PairCoupling[] + for interaction in original_interactions, pc in interaction.pair + (; bond) = pc + if bond_is_in_unit(bond, contraction_info) + push!(pcs_intra, pc) + else + push!(pcs_inter, pc) + end + end + + # Convert intra-unit PairCouplings to onsite couplings + for pc in pcs_intra + accum_pair_coupling_into_bond_operator_in_unit!(unit_operator, pc, sys, contracted_site, contraction_info) + end + set_onsite_coupling!(sys_entangled, unit_operator, contracted_site) + + ## Convert inter-unit PairCouplings into new pair couplings + for pc in pcs_inter + (; newbond, bond_operator) = pair_coupling_into_bond_operator_between_units(pc, sys, contraction_info) + push!(new_pair_data, (newbond, bond_operator)) + end + end + + # Now have list of bonds and bond operators. First we must find individual + # exemplars of each symmetry class of bonds in terms of the *new* crystal. + all_bonds_with_interactions = [data[1] for data in new_pair_data] + exemplars = Bond[] + while length(all_bonds_with_interactions) > 0 + exemplar = all_bonds_with_interactions[1] + all_bonds_with_interactions = filter(all_bonds_with_interactions) do bond + !is_related_by_symmetry(contracted_crystal, bond, exemplar) + end + push!(exemplars, exemplar) + end + + # We collected all bond_operators associated with a particular exemplar, sum + # them, and set the interaction + for bond in exemplars + relevant_interactions = filter(data -> data[1] == bond, new_pair_data) + bond_operator = sum(data[2] for data in relevant_interactions) + set_pair_coupling!(sys_entangled, bond_operator, bond) + end + + return (; sys_entangled, contraction_info) +end + + +function set_expected_dipoles_of_entangled_system!(esys) + for site in eachsite(esys.sys_origin) + set_expected_dipole_of_entangled_system!(esys, site) + end +end + +function set_expected_dipole_of_entangled_system!(esys, site) + (; sys, sys_origin, dipole_operators, source_idcs) = esys + (; dipoles) = sys_origin + + a, b, c, _ = site.I + source_idx = source_idcs[site] + Z = sys.coherents[a, b, c, source_idx] + dipoles[site] = ntuple(i -> real(dot(Z, dipole_operators[i, site], Z)), 3) + + nothing +end diff --git a/src/EntangledUnits/TypesAndAliasing.jl b/src/EntangledUnits/TypesAndAliasing.jl new file mode 100644 index 000000000..9108fd0d7 --- /dev/null +++ b/src/EntangledUnits/TypesAndAliasing.jl @@ -0,0 +1,185 @@ +################################################################################ +# Crystal contraction +################################################################################ +# Data for mapping one site inside a unit back to the site of the original +# system. +struct InverseData + site :: Int64 # Atom index of original, uncontracted crystal + offset :: Vec3 # Position offset of original atom relative to center of unit +end + +# `forward` contains a list from sites of the original crystal to a site of the +# contracted crystal, including an extra index to keep track entangled units: +# `(contracted_crystal_site_index, intra_unit_site_index)`. If the first index +# refers to a site in the new crystal that does not contain multiple units, than +# the second index will always be 1. +# +# `inverse` contains a list of length equal to the number of sites in the +# contracted crystal (corresponding to `contracted_crystal_site_index` above). +# Each element of this list is another list of tuples, +# `(site_of_original_crystal, position_offset)`. The position offset is applied +# to the position of the contracted crystal to recover the corresponding +# location in the original crystal. The length of these sublists corresponds to +# the number of sites within the entangled unit. +struct CrystalContractionInfo + forward :: Vector{Tuple{Int64, Int64}} # Original site index -> full unit index (contracted crystal site index and unit subindex) + inverse :: Vector{Vector{InverseData}} # List ordered according to contracted crystal sites. Each element is itself a list containing original crystal site indices and corresponding offset information +end + + +################################################################################ +# System +################################################################################ +struct EntangledSystem + # Entangled System, original system, and mapping info between systems + sys :: System # System containing entangled units + sys_origin :: System # Original "uncontracted" system + contraction_info :: CrystalContractionInfo # Forward and inverse mapping data for sys <-> sys_origin + + # Observable field for dipoles and mapping information for loops + dipole_operators :: Array{Matrix{ComplexF64}, 5} # An observable field corresponding to dipoles of the original system. + source_idcs :: Array{Int64, 4} # Metadata for populating the original dipoles from entangled sites. +end + +""" + EntangledSystem(sys::System{N}, units) + +Create an `EntangledSystem` from an existing `System`. `units` is a list of +tuples specifying the atoms inside each unit cell that will be grouped into a +single "entangled unit." All entangled units must lie entirely inside a unit +cell. Currently this feature is only supported for systems that can be viewed as +a regular lattice of a single unit type (all dimers, all trimers, etc). Sunny +will use the SU(_N_) formalism to model each one of these units as a distinct +Hilbert space in which the full quantum mechanical structure is locally +preserved. + +Interactions must be specified for the original `System`. Sunny will +automatically reconstruct the appropriate interactions for the +`EntangledSystem`. +""" +function EntangledSystem(sys::System{N}, units) where {N} + # Since external field is stored in the onsite interactions in + # EntangledSystems and EntangledSystems are always homogenous (i.e., + # interactions are indexed by atom/unit, not site, external field + # information must also be tracked by atom/unit index only. + for atom in axes(sys.coherents, 4) + @assert allequal(@view sys.gs[:,:,:,atom]) "`EntangledSystem` require g-factors be uniform across unit cells" + end + @assert allequal(@view sys.extfield[:,:,:,:]) "`EntangledSystems` requires a uniform applied field." + + # Generate pair of contracted and uncontracted systems + (; sys_entangled, contraction_info) = entangle_system(sys, units) + sys_origin = clone_system(sys) + + # Generate observable field. This observable field has as many entres as the + # uncontracted system but contains operators in the local product spaces of + # the contracted system. `source_idcs` provides the unit index (of the + # contracted system) in terms of the atom index (of the uncontracted + # system). + dipole_operators_origin = all_dipole_observables(sys_origin; apply_g=false) + (; observables, source_idcs) = observables_to_product_space(dipole_operators_origin, sys_origin, contraction_info) + + esys = EntangledSystem(sys_entangled, sys_origin, contraction_info, observables, source_idcs) + + # Coordinate sys_entangled and sys_origin + set_expected_dipoles_of_entangled_system!(esys) + set_field!(esys, sys.extfield[1,1,1,1]) # Note external field checked to be uniform + + return esys +end + +function EntangledSystem(::System{0}, _) + error("Cannot create an `EntangledSystem` from a `:dipole`-mode `System`.") +end + + +################################################################################ +# Aliasing +################################################################################ +function Base.show(io::IO, esys::EntangledSystem) + print(io, "EntangledSystem($(mode_to_str(esys.sys)), $(lattice_to_str(esys.sys_origin)), $(energy_to_str(esys.sys)))") +end + +function Base.show(io::IO, ::MIME"text/plain", esys::EntangledSystem) + printstyled(io, "EntangledSystem $(mode_to_str(esys.sys))\n"; bold=true, color=:underline) + println(io, supercell_to_str(esys.sys_origin.dims, esys.sys_origin.crystal)) + if !isnothing(esys.sys_origin.origin) + shape = number_to_math_string.(cell_shape(esys.sys_origin)) + println(io, formatted_matrix(shape; prefix="Reshaped cell ")) + end + println(io, energy_to_str(esys.sys)) +end + +eachsite(esys::EntangledSystem) = eachsite(esys.sys_origin) +eachunit(esys::EntangledSystem) = eachsite(esys.sys) + +energy(esys::EntangledSystem) = energy(esys.sys) +energy_per_site(esys::EntangledSystem) = energy(esys.sys) / length(eachsite(esys.sys_origin)) + +function set_field!(esys::EntangledSystem, B) + (; sys, sys_origin, dipole_operators, source_idcs) = esys + B_old = sys_origin.extfield[1,1,1,1] + set_field!(sys_origin, B) + + # Iterate through atom of original system and adjust the onsite operator of + # corresponding unit of contracted system. + for atom in axes(sys_origin.coherents, 4) + unit = source_idcs[1, 1, 1, atom] + S = dipole_operators[:, 1, 1, 1, atom] + ΔB = sys_origin.gs[1, 1, 1, atom]' * (B - B_old) + sys.interactions_union[unit].onsite -= Hermitian(ΔB' * S) + end +end + +function set_field_at!(::EntangledSystem, _, _) + error("`EntangledSystem`s do not support inhomogenous external fields. Use `set_field!(sys, B).") +end + + +set_dipole!(esys::EntangledSystem, dipole, site; kwargs...) = error("Setting dipoles of an EntangledSystem not well defined.") + +# Sets the coherent state of a specified unit. The `site` refers to the +# contracted lattice (i.e., to a "unit"). The function then updates all dipoles +# in the uncontracted system that are determined by the coherent state. +function set_coherent!(esys::EntangledSystem, coherent, site) + set_coherent!(esys.sys, coherent, site) + a, b, c, unit = site.I + for atom in atoms_in_unit(esys.contraction_info, unit) + set_expected_dipole_of_entangled_system!(esys, CartesianIndex(a, b, c, atom)) + end +end + +function randomize_spins!(esys::EntangledSystem) + randomize_spins!(esys.sys) + set_expected_dipoles_of_entangled_system!(esys) +end + +function minimize_energy!(esys::EntangledSystem; kwargs...) + optout = minimize_energy!(esys.sys; kwargs...) + set_expected_dipoles_of_entangled_system!(esys) + return optout +end + +function magnetic_moment(esys::EntangledSystem, site; kwargs...) + magnetic_moment(esys.sys_origin, site; kwargs...) +end + +function plot_spins(esys::EntangledSystem; kwargs...) + plot_spins(esys.sys_origin; kwargs...) +end + +ssf_custom(f, esys::EntangledSystem; kwargs...) = ssf_custom(f, esys.sys_origin; kwargs...) +ssf_custom_bm(f, esys::EntangledSystem; kwargs...) = ssf_custom_bm(f, esys.sys_origin; kwargs...) +ssf_trace(esys::EntangledSystem; kwargs...) = ssf_trace(esys.sys_origin; kwargs...) +ssf_perp(esys::EntangledSystem; kwargs...) = ssf_perp(esys.sys_origin; kwargs...) + +# TODO: Note this simple wrapper makes everything work, but is not the most +# efficient solution. `step!` currently syncs the dipoles field of the +# EntangledSystem in a meaninless way. This field is ignored everywhere, but +# this step represents needless computation. +function step!(esys::EntangledSystem, integrator) + step!(esys.sys, integrator) + set_expected_dipoles_of_entangled_system!(esys) +end + +suggest_timestep(esys::EntangledSystem, integrator; kwargs...) = suggest_timestep(esys.sys, integrator; kwargs...) \ No newline at end of file diff --git a/src/Measurements/MeasureSpec.jl b/src/Measurements/MeasureSpec.jl index 2da0f9aae..2a21eb813 100644 --- a/src/Measurements/MeasureSpec.jl +++ b/src/Measurements/MeasureSpec.jl @@ -7,9 +7,10 @@ struct MeasureSpec{Op <: Union{Vec3, HermitianC64}, F, Ret} corr_pairs :: Vector{NTuple{2, Int}} # (ncorr) combiner :: F # (q::Vec3, obs) -> Ret formfactors :: Array{FormFactor, 2} # (nobs × natoms) + offsets :: Array{Vec3, 2} # (nobs × natoms) # TODO: Default combiner will be SVector? - function MeasureSpec(observables::Array{Op, 5}, corr_pairs, combiner::F, formfactors) where {Op, F} + function MeasureSpec(observables::Array{Op, 5}, corr_pairs, combiner::F, formfactors; offsets=nothing) where {Op, F} # Lift return type of combiner function to type-level Ret = only(Base.return_types(combiner, (Vec3, Vector{ComplexF64}))) @assert isbitstype(Ret) @@ -17,8 +18,11 @@ struct MeasureSpec{Op <: Union{Vec3, HermitianC64}, F, Ret} if isone(ndims(formfactors)) formfactors = [ff for _ in axes(observables, 1), ff in formfactors] end - @assert size(observables)[[1,5]] == size(formfactors) - return new{Op, F, Ret}(observables, corr_pairs, combiner, formfactors) + if isnothing(offsets) + offsets = zeros(Vec3, size(observables)[[1,5]]...) + end + @assert size(observables)[[1,5]] == size(formfactors) == size(offsets) + return new{Op, F, Ret}(observables, corr_pairs, combiner, formfactors, offsets) end end @@ -27,7 +31,7 @@ function Base.show(io::IO, ::MeasureSpec) end function Base.show(io::IO, ::MIME"text/plain", m::MeasureSpec) - nobs = size(m.observables, 1) + nobs = num_observables(m) ret = eltype(m) println(io, "MeasureSpec [$nobs observables, returns $ret]") end @@ -35,7 +39,6 @@ end Base.eltype(::MeasureSpec{Op, F, Ret}) where {Op, F, Ret} = Ret -# Replicate Sam's observables code for the moment num_observables(measure::MeasureSpec) = size(measure.observables, 1) num_correlations(measure::MeasureSpec) = length(measure.corr_pairs) diff --git a/src/Operators/TensorOperators.jl b/src/Operators/TensorOperators.jl index d91c11e15..0e2da57bb 100644 --- a/src/Operators/TensorOperators.jl +++ b/src/Operators/TensorOperators.jl @@ -82,6 +82,14 @@ function svd_tensor_expansion(D::Matrix{T}, N1, N2) where T return ret end +# Given a local operator, A, that lives within an entangled unit on local site +# i, construct I ⊗ … ⊗ I ⊗ A ⊗ I ⊗ … ⊗ I, where A is in the i position. +function local_op_to_product_space(op, i, Ns) + @assert size(op, 1) == Ns[i] "Given operator not consistent with dimension of local Hilbert space" + I1 = I(prod(Ns[begin:i-1])) + I2 = I(prod(Ns[i+1:end])) + return kron(I1, op, I2) +end """ to_product_space(A, B, ...) @@ -101,10 +109,8 @@ function to_product_space(A, B, rest...) end return map(enumerate(lists)) do (i, list) - I1 = I(prod(Ns[begin:i-1])) - I2 = I(prod(Ns[i+1:end])) return map(list) do op - kron(I1, op, I2) + local_op_to_product_space(op, i, Ns) end end -end +end \ No newline at end of file diff --git a/src/SampledCorrelations/CorrelationSampling.jl b/src/SampledCorrelations/CorrelationSampling.jl index 332df34a8..767165db9 100644 --- a/src/SampledCorrelations/CorrelationSampling.jl +++ b/src/SampledCorrelations/CorrelationSampling.jl @@ -1,4 +1,4 @@ -function observable_values!(buf, sys::System{N}, observables) where N +function observable_values!(buf, sys::System{N}, observables, atom_idcs) where N if N == 0 for i in axes(observables, 1) for site in eachsite(sys) @@ -9,48 +9,46 @@ function observable_values!(buf, sys::System{N}, observables) where N end else Zs = sys.coherents - for i in axes(observables, 1) - for site in eachsite(sys) - obs = observables[i, site] - buf[i, site] = dot(Zs[site], obs, Zs[site]) - end + for idx in CartesianIndices(observables) + _, la, lb, lc, pos = idx.I + atom = atom_idcs[la, lb, lc, pos] + buf[idx] = dot(Zs[la, lb, lc, atom], observables[idx], Zs[la, lb, lc, atom]) end end return nothing end -function trajectory!(buf, sys, dt, nsnaps, observables; measperiod = 1) +function trajectory!(buf, sys, dt, nsnaps, observables, atom_idcs; measperiod = 1) @assert size(observables, 1) == size(buf, 1) integrator = ImplicitMidpoint(dt) - observable_values!(@view(buf[:,:,:,:,:,1]), sys, observables) + observable_values!(@view(buf[:,:,:,:,:,1]), sys, observables, atom_idcs) for n in 2:nsnaps for _ in 1:measperiod step!(sys, integrator) end - observable_values!(@view(buf[:,:,:,:,:,n]), sys, observables) + observable_values!(@view(buf[:,:,:,:,:,n]), sys, observables, atom_idcs) end return nothing end function new_sample!(sc::SampledCorrelations, sys::System) - (; dt, samplebuf, measperiod, measure) = sc + (; dt, samplebuf, measperiod, observables, atom_idcs) = sc # Only fill the sample buffer half way; the rest is zero-padding buf_size = size(samplebuf, 6) nsnaps = (buf_size÷2) + 1 samplebuf[:,:,:,:,:,(nsnaps+1):end] .= 0 - @assert size(sys.dipoles) == size(samplebuf)[2:5] "`System` size not compatible with given `SampledCorrelations`" + # @assert size(sys.dipoles) == size(samplebuf)[2:5] "`System` size not compatible with given `SampledCorrelations`" - trajectory!(samplebuf, sys, dt, nsnaps, measure.observables; measperiod) + trajectory!(samplebuf, sys, dt, nsnaps, observables, atom_idcs; measperiod) return nothing end function accum_sample!(sc::SampledCorrelations; window) - (; data, M, measure, samplebuf, corrbuf, nsamples, space_fft!, time_fft!, corr_fft!, corr_ifft!) = sc - natoms = size(samplebuf)[5] - + (; data, M, corr_pairs, samplebuf, corrbuf, space_fft!, time_fft!, corr_fft!, corr_ifft!) = sc + npos = size(samplebuf)[5] num_time_offsets = size(samplebuf, 6) T = (num_time_offsets÷2) + 1 # Duration that each signal was recorded for @@ -83,7 +81,7 @@ function accum_sample!(sc::SampledCorrelations; window) count = sc.nsamples += 1 - for j in 1:natoms, i in 1:natoms, (c, (α, β)) in enumerate(measure.corr_pairs) + for j in 1:npos, i in 1:npos, (c, (α, β)) in enumerate(corr_pairs) # α, β = ci.I sample_α = @view samplebuf[α,:,:,:,i,:] diff --git a/src/SampledCorrelations/DataRetrieval.jl b/src/SampledCorrelations/DataRetrieval.jl index 1b89fd46d..89f32ee6d 100644 --- a/src/SampledCorrelations/DataRetrieval.jl +++ b/src/SampledCorrelations/DataRetrieval.jl @@ -100,10 +100,11 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT q_idx_info = pruned_wave_vector_info(sc, qs_reshaped) crystal = @something sc.origin_crystal sc.crystal NCorr = Val{size(sc.data, 1)}() - NAtoms = Val{size(sc.data, 2)}() + # NPos = Val{size(sc.data, 2)}() + NPos = Val{length(sc.crystal.positions)}() # Intensities calculation - intensities_rounded!(intensities, sc.data, sc.crystal, sc.measure, ffs, q_idx_info, ωidcs, NCorr, NAtoms) + intensities_rounded!(intensities, sc.data, sc.crystal, sc.positions, sc.measure.combiner, ffs, q_idx_info, ωidcs, NCorr, NPos) # Convert to a q-space density in original (not reshaped) RLU. intensities .*= det(sc.crystal.recipvecs) / det(crystal.recipvecs) @@ -133,19 +134,20 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT end end -function intensities_rounded!(intensities, data, crystal, measure::MeasureSpec{Op, F, Ret}, ff_atoms, q_idx_info, ωidcs, ::Val{NCorr}, ::Val{NAtoms}) where {Op, F, Ret, NCorr, NAtoms} +function intensities_rounded!(intensities, data, crystal, positions, combiner, ff_atoms, q_idx_info, ωidcs, ::Val{NCorr}, ::Val{NPos}) where {NCorr, NPos} (; qabs, idcs, counts) = q_idx_info + (; recipvecs) = crystal qidx = 1 for (qabs, idx, count) in zip(qabs, idcs, counts) - prefactors = prefactors_for_phase_averaging(qabs, crystal, ff_atoms, Val(NCorr), Val(NAtoms)) + prefactors = prefactors_for_phase_averaging(qabs, recipvecs, @view(positions[idx,:]), ff_atoms, Val(NCorr), Val(NPos)) # Perform phase-averaging over all omega for (n, iω) in enumerate(ωidcs) elems = zero(MVector{NCorr,ComplexF64}) - for j in 1:NAtoms, i in 1:NAtoms + for j in 1:NPos, i in 1:NPos elems .+= (prefactors[i] * conj(prefactors[j])) .* view(data, :, i, j, idx, iω) end - val = measure.combiner(qabs, SVector{NCorr, ComplexF64}(elems)) + val = combiner(qabs, SVector{NCorr, ComplexF64}(elems)) intensities[n, qidx] = val end diff --git a/src/SampledCorrelations/PhaseAveraging.jl b/src/SampledCorrelations/PhaseAveraging.jl index 3e952e12d..96f528ec1 100644 --- a/src/SampledCorrelations/PhaseAveraging.jl +++ b/src/SampledCorrelations/PhaseAveraging.jl @@ -16,13 +16,13 @@ function phase_averaged_elements(data, q_absolute::Vec3, crystal::Crystal, ff_at return SVector{NCorr,ComplexF64}(elems) end -function prefactors_for_phase_averaging(q_absolute::Vec3, crystal::Crystal, ff_atoms, ::Val{NCorr}, ::Val{NAtoms}) where {NCorr, NAtoms} +function prefactors_for_phase_averaging(q_absolute::Vec3, recipvecs, positions, ff_atoms, ::Val{NCorr}, ::Val{NAtoms}) where {NCorr, NAtoms} # Form factor ffs = ntuple(i -> compute_form_factor(ff_atoms[i], q_absolute⋅q_absolute), Val{NAtoms}()) # Overall phase factor for each site - q = crystal.recipvecs \ q_absolute - r = crystal.positions + q = recipvecs \ q_absolute + r = positions prefactor = ntuple(i -> ffs[i] * exp(- 2π*im * (q ⋅ r[i])), Val{NAtoms}()) return prefactor @@ -38,7 +38,7 @@ end # not used, but would be employed in a parallel intensities-like pipeline for # error propagation. function error_basis_reduction(data, q_absolute::Vec3, _::Crystal, ff_atoms, ::Val{NCorr}, ::Val{NAtoms}) where {NCorr, NAtoms} - elems = zero(MVector{NCorr,ComplexF64}) + elems = zero(MVector{NCorr, ComplexF64}) ffs = ntuple(i -> compute_form_factor(ff_atoms[i], q_absolute⋅q_absolute), NAtoms) for j in 1:NAtoms, i in 1:NAtoms diff --git a/src/SampledCorrelations/SampledCorrelations.jl b/src/SampledCorrelations/SampledCorrelations.jl index 284548b0d..0d8e9decf 100644 --- a/src/SampledCorrelations/SampledCorrelations.jl +++ b/src/SampledCorrelations/SampledCorrelations.jl @@ -4,8 +4,14 @@ mutable struct SampledCorrelations const M :: Union{Nothing, Array{Float64, 7}} # Running estimate of (nsamples - 1)*σ² (where σ² is the variance of intensities) const crystal :: Crystal # Crystal for interpretation of q indices in `data` const origin_crystal :: Union{Nothing,Crystal} # Original user-specified crystal (if different from above) - const Δω :: Float64 # Energy step size (could make this a virtual property) - measure :: MeasureSpec # Observable, correlation pairs, and combiner + const Δω :: Float64 # Energy step size + + # Observable information + measure :: MeasureSpec # Storehouse for combiner. Mutable so combiner can be changed. + const observables # :: Array{Op, 5} # (nobs × npos x latsize) -- note change of ordering relative to MeasureSpec. TODO: determine type strategy + const positions :: Array{Vec3, 4} # Position of each operator in fractional coordinates (latsize x npos) + const atom_idcs :: Array{Int64, 4} # Atom index corresponding to position of observable. + const corr_pairs :: Vector{NTuple{2, Int}} # (ncorr) # Trajectory specs const measperiod :: Int # Steps to skip between saving observables (i.e., downsampling factor for trajectories) @@ -49,7 +55,8 @@ function clone_correlations(sc::SampledCorrelations) corr_ifft! = FFTW.plan_ifft!(sc.corrbuf, 4) M = isnothing(sc.M) ? nothing : copy(sc.M) return SampledCorrelations( - copy(sc.data), M, sc.crystal, sc.origin_crystal, sc.Δω, deepcopy(sc.measure), + copy(sc.data), M, sc.crystal, sc.origin_crystal, sc.Δω, + deepcopy(sc.measure), copy(sc.observables), copy(sc.positions), copy(sc.atom_idcs), copy(sc.corr_pairs), sc.measperiod, sc.dt, sc.nsamples, copy(sc.samplebuf), copy(sc.corrbuf), space_fft!, time_fft!, corr_fft!, corr_ifft! ) @@ -126,7 +133,7 @@ can can then be extracted as pair-correlation [`intensities`](@ref) with appropriate classical-to-quantum correction factors. See also [`intensities_static`](@ref), which integrates over energy. """ -function SampledCorrelations(sys::System; measure, energies, dt, calculate_errors=false) +function SampledCorrelations(sys::System; measure, energies, dt, calculate_errors=false, positions=nothing) if isnothing(energies) n_all_ω = 1 measperiod = 1 @@ -143,17 +150,31 @@ function SampledCorrelations(sys::System; measure, energies, dt, calculate_error Δω = ωmax/(nω-1) end - # Preallocation - na = natoms(sys.crystal) + # Determine the positions of the observables in the MeasureSpec. By default, + # these will just be the atom indices. + positions = if isnothing(positions) + map(eachsite(sys)) do site + sys.crystal.positions[site.I[4]] + end + else + positions + end + + # Determine the number of positions. For an unentangled system, this will + # just be the number of atoms. + npos = size(positions, 4) + + # Determine which atom index is used to derive information about a given + # physical position. This becomes relevant for entangled units. + atom_idcs = map(site -> site.I[4], eachsite(sys)) - # The sample buffer holds n_non_neg_ω measurements, and the rest is a zero buffer measure = isnothing(measure) ? ssf_trace(sys) : measure num_observables(measure) - samplebuf = zeros(ComplexF64, num_observables(measure), sys.dims..., na, n_all_ω) + samplebuf = zeros(ComplexF64, num_observables(measure), sys.dims..., npos, n_all_ω) corrbuf = zeros(ComplexF64, sys.dims..., n_all_ω) # The output data has n_all_ω many (positive and negative and zero) frequencies - data = zeros(ComplexF64, num_correlations(measure), na, na, sys.dims..., n_all_ω) + data = zeros(ComplexF64, num_correlations(measure), npos, npos, sys.dims..., n_all_ω) M = calculate_errors ? zeros(Float64, size(data)...) : nothing # The normalization is defined so that the prod(sys.dims)-many estimates of @@ -173,7 +194,9 @@ function SampledCorrelations(sys::System; measure, energies, dt, calculate_error # Make Structure factor and add an initial sample origin_crystal = isnothing(sys.origin) ? nothing : sys.origin.crystal - sc = SampledCorrelations(data, M, sys.crystal, origin_crystal, Δω, measure, measperiod, dt, nsamples, + sc = SampledCorrelations(data, M, sys.crystal, origin_crystal, Δω, + measure, copy(measure.observables), positions, atom_idcs, copy(measure.corr_pairs), + measperiod, dt, nsamples, samplebuf, corrbuf, space_fft!, time_fft!, corr_fft!, corr_ifft!) return sc @@ -192,13 +215,12 @@ excitation spectrum cannot be performed. """ struct SampledCorrelationsStatic parent :: SampledCorrelations - - function SampledCorrelationsStatic(sys::System; measure, calculate_errors=false) - parent = SampledCorrelations(sys; measure, energies=nothing, dt=NaN, calculate_errors) - return new(parent) - end end +function SampledCorrelationsStatic(sys::System; measure, calculate_errors=false) + parent = SampledCorrelations(sys; measure, energies=nothing, dt=NaN, calculate_errors) + return SampledCorrelationsStatic(parent) +end function Base.show(io::IO, ::SampledCorrelations) print(io, "SampledCorrelations") @@ -211,24 +233,25 @@ end function Base.show(io::IO, ::MIME"text/plain", sc::SampledCorrelations) - (; crystal, sys_dims, nsamples) = sc - nω = Int(size(sc.data, 7)/2+1) - Δω = round(sc.Δω, digits=4) + (; crystal, nsamples) = sc + nω = round(Int, size(sc.data)[7]/2) + sys_dims = size(sc.data[4:6]) printstyled(io, "SampledCorrelations"; bold=true, color=:underline) println(io," ($(Base.format_bytes(Base.summarysize(sc))))") print(io,"[") printstyled(io,"S(q,ω)"; bold=true) - print(io," | nω = $nω, Δω = $Δω") + print(io," | nω = $nω, Δω = $(round(sc.Δω, digits=4))") println(io," | $nsamples $(nsamples > 1 ? "samples" : "sample")]") - println(io, supercell_to_str(sys_dims, crystal)) + println(io,"Lattice: $sys_dims × $(natoms(crystal))") end function Base.show(io::IO, ::MIME"text/plain", sc::SampledCorrelationsStatic) - (; crystal, sys_dims, nsamples) = sc.parent + (; crystal, nsamples) = sc.parent + sys_dims = size(sc.parent.data[4:6]) printstyled(io, "SampledCorrelationsStatic"; bold=true, color=:underline) println(io," ($(Base.format_bytes(Base.summarysize(sc))))") print(io,"[") printstyled(io,"S(q)"; bold=true) println(io," | $nsamples $(nsamples > 1 ? "samples" : "sample")]") - println(io, supercell_to_str(sys_dims, crystal)) -end + println(io,"Lattice: $sys_dims × $(natoms(crystal))") +end \ No newline at end of file diff --git a/src/SpinWaveTheory/DispersionAndIntensities.jl b/src/SpinWaveTheory/DispersionAndIntensities.jl index 6c2684b15..d0dc8c935 100644 --- a/src/SpinWaveTheory/DispersionAndIntensities.jl +++ b/src/SpinWaveTheory/DispersionAndIntensities.jl @@ -162,27 +162,28 @@ function intensities_bands(swt::SpinWaveTheory, qpts; kT=0, with_negative=false) # Number of wavevectors Nq = length(qpts.qs) + # Temporary storage for pair correlations + Nobs = num_observables(measure) + Ncorr = num_correlations(measure) + corrbuf = zeros(ComplexF64, Ncorr) + # Preallocation T = zeros(ComplexF64, 2L, 2L) H = zeros(ComplexF64, 2L, 2L) - Avec_pref = zeros(ComplexF64, Na) + Avec_pref = zeros(ComplexF64, Nobs, Na) disp = zeros(Float64, L, Nq) intensity = zeros(eltype(measure), L, Nq) - # Temporary storage for pair correlations - Nobs = size(measure.observables, 1) - Ncorr = length(measure.corr_pairs) - corrbuf = zeros(ComplexF64, Ncorr) for (iq, q) in enumerate(qpts.qs) q_global = cryst.recipvecs * q view(disp, :, iq) .= view(excitations!(T, H, swt, q), 1:L) - for i in 1:Na - r_global = global_position(sys, (1,1,1,i)) - ff = get_swt_formfactor(measure, 1, i) - Avec_pref[i] = exp(- im * dot(q_global, r_global)) - Avec_pref[i] *= compute_form_factor(ff, norm2(q_global)) + for i in 1:Na, μ in 1:Nobs + r_global = global_position(sys, (1,1,1,i)) # + offsets[μ,i] + ff = get_swt_formfactor(measure, μ, i) + Avec_pref[μ, i] = exp(- im * dot(q_global, r_global)) + Avec_pref[μ, i] *= compute_form_factor(ff, norm2(q_global)) end Avec = zeros(ComplexF64, Nobs) @@ -197,7 +198,7 @@ function intensities_bands(swt::SpinWaveTheory, qpts; kT=0, with_negative=false) for i in 1:Na, μ in 1:Nobs O = data.observables_localized[μ, i] for α in 1:N-1 - Avec[μ] += Avec_pref[i] * (O[α, N] * t[α, i, 2] + O[N, α] * t[α, i, 1]) + Avec[μ] += Avec_pref[μ, i] * (O[α, N] * t[α, i, 2] + O[N, α] * t[α, i, 1]) end end else @@ -211,7 +212,7 @@ function intensities_bands(swt::SpinWaveTheory, qpts; kT=0, with_negative=false) # local frame, z is longitudinal, and we are computing # the transverse part only, so the last entry is zero) displacement_local_frame = SA[t[i, 2] + t[i, 1], im * (t[i, 2] - t[i, 1]), 0.0] - Avec[μ] += Avec_pref[i] * (data.sqrtS[i]/sqrt(2)) * (O' * displacement_local_frame)[1] + Avec[μ] += Avec_pref[μ, i] * (data.sqrtS[i]/sqrt(2)) * (O' * displacement_local_frame)[1] end end diff --git a/src/Sunny.jl b/src/Sunny.jl index d0c07ea43..35d117f3b 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -111,6 +111,14 @@ include("SampledCorrelations/DataRetrieval.jl") export SampledCorrelations, SampledCorrelationsStatic, add_sample!, clone_correlations, merge_correlations +include("EntangledUnits/TypesAndAliasing.jl") +include("EntangledUnits/EntangledUnits.jl") +include("EntangledUnits/EntangledReshaping.jl") +include("EntangledUnits/EntangledSpinWaveTheory.jl") +include("EntangledUnits/EntangledSampledCorrelations.jl") +# export contract_crystal, EntangledSystem, set_expected_dipoles_of_entangled_system! +# export EntangledSpinWaveTheory, EntangledSampledCorrelations + include("MonteCarlo/Samplers.jl") include("MonteCarlo/BinnedArray.jl") include("MonteCarlo/ParallelTempering.jl") diff --git a/test/test_entangled_units.jl b/test/test_entangled_units.jl new file mode 100644 index 000000000..0116b679d --- /dev/null +++ b/test/test_entangled_units.jl @@ -0,0 +1,149 @@ +@testitem "Crystal contraction and expansion" begin + crystal = Sunny.diamond_crystal() + + # Specify various entangled units. + units_all = [ + [(1, 3), (2, 4), (5, 7), (6, 8)], + [(1, 3, 5, 6), (2, 4, 7, 8)], + [(1, 3)] + ] + + # Check that re-expansion of a contracted crystal matches original crystal + # in terms of site-ordering and positions. + for units in units_all + contracted_crystal, contraction_info = Sunny.contract_crystal(crystal, units) + expanded_crystal = Sunny.expand_crystal(contracted_crystal, contraction_info) + @test expanded_crystal.positions ≈ crystal.positions + end +end + +# TODO: Add reshapings tests. + +# TODO: Add test with magnetic unit cell larger than a single unit (i.e. not q=0 +# ordering). + +@testitem "Test Dimer Interactions" begin + J = 1.0 + J′ = 0.1 + latvecs = [ + 1 0 0 + 0 1 0 + 0 0 2 + ] + positions = [[0, 0, 0], [0.0, 0.5, 0.0]] + + crystal = Crystal(latvecs, positions, 1; types = ["A", "B"]) + sys = System(crystal, [1 => Moment(s=1/2, g=2), 2 => Moment(s=1/2, g=2)], :SUN) + set_exchange!(sys, J, Bond(1, 2, [0, 0, 0])) + set_exchange!(sys, J′, Bond(1, 1, [1, 0, 0])) + set_exchange!(sys, J′, Bond(2, 2, [1, 0, 0])) # Needed because we broke the symmetry equivalence of the two sites + + esys = Sunny.EntangledSystem(sys, [(1, 2)]) + interactions = esys.sys.interactions_union[1] + + # Test on-bond exchange + onsite_operator = interactions.onsite + S = spin_matrices(1/2) + Sl, Su = to_product_space(S, S) + onsite_ref = J * (Sl' * Su) + @test onsite_operator ≈ onsite_ref + + # Test external field works as expected + set_field!(esys, [0, 0, 1]) + onsite_operator = esys.sys.interactions_union[1].onsite + field_offset = -2*(Sl[3] + Su[3]) # 2 for g-factor + @test onsite_operator ≈ onsite_ref + field_offset + + # Test inter-bond exchange + pc = Sunny.as_general_pair_coupling(interactions.pair[1], esys.sys) + Sl1, Sl2 = to_product_space(Sl, Sl) + Su1, Su2 = to_product_space(Su, Su) + bond_operator = zeros(ComplexF64, 16, 16) + for (A, B) in pc.general.data + bond_operator .+= kron(A, B) + end + bond_ref = J′*((Sl2' * Sl1) .+ (Su2' * Su1)) + @test bond_operator ≈ bond_ref +end + +# @testitem "Ba3Mn2O8 Dispersion and Intensities" begin +# function J_of_Q(_, J1, J2, J3, J4, q) +# h, k, l = q +# ω1(h, k, l) = cos((2π/3)*(-h + k + l)) + cos((2π/3)*(-h - 2k + l)) + cos((2π/3)*(2h + k + l)) +# ω2(h, k, l) = cos(2π*k) + cos(2π*(h + k)) + cos(2π*h) +# ω4(h, k, l) = cos((2π/3)*(2h - 2k + l)) + cos((2π/3)*(2h + 4k + l)) + cos((2π/3)*(-4h - 2k + l)) +# -J1*ω1(h, k, l) + 2(J2 - J3)*ω2(h, k, l) - J4*ω4(h, k, l) +# end +# +# function disp0(q) +# J0, J1, J2, J3, J4, D = 1.642, 0.118, 0.256, 0.142, 0.037, -0.032 +# Δ0 = J0 + 2D/3 +# sqrt(Δ0^2 + (8/3)*Δ0*J_of_Q(J0, J1, J2, J3, J4, q) ) +# end +# +# function dispm(q) +# J0, J1, J2, J3, J4, D = 1.642, 0.118, 0.256, 0.142, 0.037, -0.032 +# Δm = J0 - D/3 +# sqrt(Δm^2 + (8/3)*Δm*J_of_Q(J0, J1, J2, J3, J4, q) ) +# end +# +# function Mn_crystal() +# latvecs = [5.710728 -2.855364 0.0; 0.0 4.945635522103099 0.0; 0.0 0.0 21.44383] +# positions = [[1/3, 2/3, 0.07374666666666664], [1/3, 2/3, 0.25958666666666663], [0.0, 0.0, 0.40708], [0.0, 0.0, 0.59292], [2/3, 1/3, 0.7404133333333334], [2/3, 1/3, 0.9262533333333333]] +# symprec = 0.01 +# return Crystal(latvecs, positions; symprec) +# end +# +# function Ba3Mn2O8(; dims=(1, 1, 1), g=1.0, D = -0.032, J0 = 1.642, J1 = 0.118, J2 = 0.256, J3 = 0.142, J4 = 0.037) +# sys = System(Mn_crystal(), dims, [SpinInfo(1; S=1, g)], :SUN) +# S = spin_matrices(1) +# S1, S2 = to_product_space(S, S) +# set_pair_coupling!(sys, J0 * (S1' * S2), Bond(1, 2, [0, 0, 0])) +# set_pair_coupling!(sys, J1 * (S1' * S2), Bond(3, 2, [0, 0, 0])) +# set_pair_coupling!(sys, J2 * (S1' * S2), Bond(1, 1, [1, 0, 0])) +# set_pair_coupling!(sys, J3 * (S1' * S2), Bond(1, 2, [1, 0, 0])) +# set_pair_coupling!(sys, J4 * (S1' * S2), Bond(3, 2, [1, 0, 0])) +# set_onsite_coupling!(sys, D*S[3]^2, 1) +# return sys +# end +# +# # Set up system and reshape into primitive cell +# crystal = Mn_crystal() +# sys = Ba3Mn2O8() +# shape = [-1 -1 2/3; 0 -1 1/3; 0 0 1/3] +# sys_reshaped = reshape_supercell(sys, shape) +# +# # Make entangled system and spin wave theory +# sys_entangled = EntangledSystem(sys_reshaped, [(1, 2)]) +# set_coherent!(sys_entangled.sys, Sunny.CVec{9}(0.0, 0.0, √3/3, 0.0, -√3/3, 0.0, √3/3, 0, 0), (1, 1, 1, 1)) +# swt = EntangledSpinWaveTheory(sys_entangled) +# +# # Calculation dispersion and intensities +# FWHM = 0.295 +# +# points = [ +# [0.175, 0.175, 1.5], +# [0.85, 0.85, 1.5], +# [0.85, 0.85, 3], +# [0.0, 0.0, 3], +# [0.0, 0.0, 8], +# ] +# qpts = q_space_path(crystal, points, 20) +# energies = range(0.0, 4.0, 10) +# +# is = intensities(swt, qpts; energies; kernel=gaussian(; fwhm=FWHM)) +# is_bands = intensities_bands(swt, qpts) +# +# d0_ref = disp0.(path) +# d0 = disp[:,8] +# band_err1 = √sum((d0_ref .- d0) .^ 2) +# @test band_err1 < 1e-7 +# +# dm_ref = dispm.(path) +# dm = disp[:,6] +# band_err2 = √sum((dm_ref .- dm) .^ 2) +# @test band_err2 < 1e-7 +# +# is_ref = [1.5318559039289553e-54 5.4757218095775535e-33 7.357403275541178e-17 3.9795517317276186e-6 0.9020338263369125 0.8314157146337926 2.911058313640675e-6 3.665251511398035e-17 1.6108014399323818e-33 1.2770658651600437e-42; 6.344775971921163e-16 1.8297903880834307e-5 2.186860140613253 1.125414819762803 2.372670614798904e-6 1.8884749639385228e-17 5.364594261005178e-34 5.293626361366888e-56 1.5507874038885882e-56 1.5731739245245433e-42; 1.916445283080627e-22 1.5112603523543528e-9 0.046321847672565045 5.889951708200119 0.0031331962838944093 6.5822336276008166e-12 5.1035843246680587e-26 1.398942793394379e-45 1.4068300720405523e-55 1.0612068475874395e-41; 1.0342024392072121e-16 6.295383373660083e-6 1.5314917741949632 1.5884648898847082 6.9292015386003456e-6 1.1783549059185498e-16 7.29390986501124e-33 1.5826736515812515e-54 7.205550912747091e-56 5.5831854222249985e-42; 7.816369005435853e-65 4.958135881093393e-41 1.1004603393444979e-22 8.837745291621732e-10 0.02729382866823639 3.4792844815914976 0.001870907449751048 4.029919472431617e-12 3.236721779435184e-26 1.258781799195952e-42; 4.1286086857912686e-32 4.109912848294968e-16 1.53865152859326e-5 2.312737091532787 1.4474386749428667 3.6605144080508623e-6 3.5053122616269065e-17 1.205982396647456e-33 1.888239550856115e-55 3.384617657031145e-42; 1.0956905921250642e-15 3.09807215254437e-5 3.6242252907637944 1.8260559377579642 3.77624592541727e-6 2.953182889984101e-17 8.250224084439066e-34 8.009410306835117e-56 4.543889223635611e-56 5.863106256174803e-42; 2.4641979663556544e-21 9.259498870544819e-9 0.1405047201368198 8.991804007003132 0.0023532144516871254 2.345225668006157e-12 8.408440394132043e-27 1.0520889028677954e-46 1.508356309761414e-55 1.1755786936605455e-41; 5.9006127773674436e-21 1.6083917235355553e-8 0.18138815381567364 8.612130243663842 0.0016307094799296721 1.1473990136958224e-12 2.86349056374958e-27 2.4776173823909797e-47 1.4600597514018359e-56 1.2635254358709617e-42; 3.534719908717514e-16 1.5534826875777495e-5 2.933172504328625 2.2501936681347714 6.4727258857098395e-6 6.621955526402235e-17 2.3498631523198714e-33 2.862774376331706e-55 9.224249931423226e-57 1.3867670560912468e-42; 2.9056707554206643e-34 1.1389749648037732e-17 1.8359864206642449e-6 1.1476609410649823 2.6195576292904894 2.104174627891465e-5 5.842573811256141e-16 5.564262837339542e-32 1.8126351339556568e-53 9.514020434069003e-43; 4.95669586498486e-51 4.090111065132001e-30 1.1512703563038524e-14 0.00011053986637697482 3.620428209362293 0.40448305917267885 1.541487919718299e-7 2.0039150465080363e-19 8.886236730407736e-37 8.03538987497683e-43; 1.6256804884099134e-104 1.626131554260588e-73 5.548499522448917e-48 6.457952119615799e-28 2.563973194366388e-13 0.0003472414267885618 1.6041647745239174 0.025279344649153018 1.358882914693268e-9 2.491710931378807e-22; 1.4445317508683965e-67 1.425287573157436e-43 4.79708370999404e-25 5.507461722013346e-12 0.00021568754619639152 0.02881363456352638 1.3130171839110966e-5 2.0409960671735094e-14 1.0822144233620588e-28 1.799465203023977e-45; 9.225789423317937e-80 5.21566646643929e-53 1.0058089054657094e-31 6.616381427047744e-16 1.4846530782284771e-5 1.1363937177899095 0.29670991367977784 2.64261978483415e-7 8.028542283600199e-19 8.320283660276395e-36; 1.0913958451883772e-93 2.04937342166082e-64 1.312681289078253e-40 2.868117121788508e-22 2.1376365239065535e-9 0.05434633134660532 4.713092846351799 0.0013942515956084163 1.4069403885772203e-12 4.8429461011708685e-27] +# @test isapprox(is, is_ref; atol=1e-12) +# end \ No newline at end of file