Reference

Extensions for making Sunny.jl usable

B_isotropic_from_phonons(model, Φ; T, cryst, qgrid=(12,12,12), q_cell=:primitive, eps_meV=1e-6)

Return isotropic Debye–Waller B-factors (Ų) per site from first principles: Bs(T) = 8π² ⟨us^2⟩(T).

assemble_force_constants!(model)

Builds real-space force constants Φ as a dictionary of 3×3 blocks keyed by (i,j,R). Per-bond block: K = kL*(êêᵀ) + kT*(I - êêᵀ) with ê = r0/‖r0‖, applied to (uj@R - ui@0). Conservation laws are enforced later by enforce_asr!.

build_model(crystal; mass, neighbors_sunny=nothing, neighbors=nothing, cutoff=nothing, use_sunny_radius=true, kL=1.0, kT=1.0, supercell=(1,1,1))

Create a Model from a Sunny-like crystal object. Minimal interface expected:

• Either Sunny-style fields: crystal.latvecs (3×3), crystal.positions (fractional), crystal.types (Strings)
• or generic fields: crystal.lattice (3×3), crystal.positions (fractional), crystal.species (Symbols/Strings)

Bond assembly precedence:

1. If neighbors_sunny (e.g. Vector{Sunny.Bond}) is provided → fast path via Sunny geometry.
2. Else if neighbors (list of tuples as documented in neighbor_bonds_from_sunny) is provided → fast tuple path.
3. Else if use_sunny_radius && cutoff!=nothing && Sunny loaded → build via Sunny.all_bonds_for_atom.
4. Else → fallback to portable cutoff scan neighbor_bonds_cutoff.

dynamical_matrix(model, Φ, q_cart)

Mass-weighted dynamical matrix D(q) from real-space Φ blocks. q_cart is Cartesian (Å⁻¹). Returns Hermitian matrix.

enforce_asr!(Φ, N)

Acoustic sum rule: for each atom i, ∑{j,R} Φ{i j}(R) = 0. Adjust on-site blocks to satisfy translational invariance.

msd_from_phonons(model, Φ; T, cryst, qgrid=(12,12,12), q_cell=:primitive, eps_meV=1e-6)

Compute per-site mean-square displacements ⟨us^2⟩(T) in Ų directly from the phonon spectrum: ⟨us^2⟩ = (1/Nq) * Σ{q,ν} (ħ² / (2 Ms EJ)) * coth(EJ / (2 kB T)) * Σα |ephys(s,α;q,ν)|² where EJ = EmeV ⋅ meVtoJ and ephys are physical (non mass-weighted) polarization vectors. Modes with E ≤ eps_meV are skipped (Γ handling); refine qgrid for convergence.

neighbor_bonds_cutoff(lattice, fracpos; cutoff, kL, kT, supercell=(1,1,1))

Construct bonds using a radial cutoff by scanning a small supercell (portable path). kL and kT can be numbers or functions (i,j,rij)->stiffness for non-uniform springs. Returns Vector{Bond}.

neighbor_bonds_from_sunny(cryst, bonds; kL=1.0, kT=1.0)

Fast-path bond construction using Sunny's bonds and geometry. bonds may be a Vector of:

• Sunny.Bond (fields: i,j,n)
• Tuples (i,j,n::SVector{3,Int})
• Tuples (i,j,Δ::SVector{3,Float64}) # direct Cartesian Δ

kL/kT can be numbers or (i,j,r0)->value. Returns Vector{Bond}.

neighbor_bonds_radius(cryst; cutoff, kL, kT)

Sunny-powered radius query: builds bonds using Sunny.all_bonds_for_atom. Requires using Sunny. Returns Vector{Bond}.

onephonon_dsf(model, Φ, q, Evals; T=300.0, bcoh=nothing, B=nothing, η=0.5, mass_unit=:amu, q_basis=:cart, q_cell=:primitive, cryst=nothing)

Computes one-phonon coherent S(q,E) on energy grid Evals (meV).

• Debye–Waller is intrinsic: full anisotropic tensors Us(T) are computed from the phonon spectrum (no ΘD).
• bcoh optional Vector of coherent scattering lengths (length N)
• η Gaussian HWHM for energy broadening (meV)
• mass_unit :amu or :kg (Model.mass is assumed to be :amu typically)
• q_basis/q_cell/cryst follow phonons conventions

Returns S(E) (arbitrary but consistent units).

onephonon_dsf_4d(model, Φ, q1, q2, q3, Evals;
                 q_basis=:rlu, q_cell=:primitive, cryst=nothing,
                 T=300.0, η=0.5, mass_unit=:amu, bcoh=nothing,
                 threads=true)

Compute the one-phonon coherent dynamic structure factor on a 4D grid:

• q1, q2, q3 are 1D axes defining a rectilinear grid in q-space.
  • If q_basis=:rlu, they are fractional coordinates along the three reciprocal basis vectors.
  • If q_basis=:cart, they are Cartesian components (Å⁻¹).
• Evals is the 1D energy axis (meV). Result has shape (length(q1), length(q2), length(q3), length(Evals)).

This routine is thread-parallel over the q-grid.

phonons(model, Φ, q; q_basis=:cart, q_cell=:primitive)

Diagonalize D(q). Returns (E, Evec) with energies in meV and 3N×3N eigenvectors. q may be in Cartesian Å⁻¹ (q_basis=:cart) or RLU (q_basis=:rlu with q_cell=:primitive or :conventional). If passing RLU, also pass the underlying crystal cryst via keyword.

q_cartesian(cryst, q; basis=:cart, cell=:primitive)

Convert a wavevector to Cartesian (Å⁻¹). q may be:

• basis=:cart → already Cartesian (Å⁻¹)
• basis=:rlu → reduced lattice units (h,k,l). cell chooses :primitive or :conventional.

Uses Sunny.prim_recipvecs(cryst) if Sunny is loaded and cell=:primitive, otherwise falls back to 2π * inv(L)' where L is the corresponding direct lattice.