Skip to content

Commit

Permalink
Allow reshaping an inhomogeneous system to a single unit cell (#191)
Browse files Browse the repository at this point in the history 
Before this commit, reshaping an inhomogeneous system was disallowed. Now, it is allowed in certain specialized contexts. The primary motivation is enable LSWT-KPM for large, inhomogeneous systems. To support this, we need to allow reshaping of an inhomogeneous system into a single, enlarged "crystal" unit cell. With the refactors, it also becomes possible to call `repeat_periodically` on an inhomogeneous system. As a design decision, `reshape_supercell` is still disallowed for inhomogeneous systems, because there is not yet an obvious use case for this.

Co-authored-by: David A. Dahlbom <34151623+ddahlbom@users.noreply.github.com>
  • Loading branch information
@kbarros @ddahlbom
kbarros and ddahlbom authored Nov 15, 2023
1 parent 863761f commit 72b35a6
Show file tree
Hide file tree
Showing 6 changed files with 110 additions and 116 deletions.
140 changes: 63 additions & 77 deletions src/Reshaping.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,43 +8,49 @@ periodicity of a requested supercell. The columns of the ``3×3`` integer matrix
original crystal lattice vectors.
"""
function reshape_supercell(sys::System{N}, shape) where N
is_homogeneous(sys) || error("Cannot reshape system with inhomogeneous interactions.")

orig = orig_crystal(sys)
check_shape_commensurate(orig, shape)

supervecs = orig.latvecs * shape
check_latvecs_commensurate(orig, supervecs)

# Convert `shape` to multiples of primitive lattice vectors
if !isnothing(orig.prim_latvecs)
prim_shape = orig.prim_latvecs \ supervecs
prim_shape′ = round.(Int, prim_shape)
@assert prim_shape′ ≈ prim_shape
new_latsize = NTuple{3, Int}(gcd.(eachcol(prim_shape′)))
else
shape′ = round.(Int, shape)
@assert shape′ ≈ shape
new_latsize = NTuple{3, Int}(gcd.(eachcol(shape′)))
end
primvecs = @something orig.prim_latvecs orig.latvecs
prim_shape = primvecs \ orig.latvecs * shape
prim_shape′ = round.(Int, prim_shape)
@assert prim_shape′ ≈ prim_shape

# Unit cell for new system, in units of original unit cell.
new_cell_shape = Mat3(shape * diagm(collect(inv.(new_latsize))))
enlarged_latsize = NTuple{3, Int}(gcd.(eachcol(prim_shape′)))
reduced_shape = Mat3(shape * diagm(collect(inv.(enlarged_latsize))))

return reshape_supercell_aux(sys, new_latsize, new_cell_shape)
return reshape_supercell_aux(sys, enlarged_latsize, reduced_shape)
end


# Transfer homogeneous interactions from `sys.origin` to reshaped `sys`
function set_interactions_from_origin!(sys::System{N}) where N
origin = sys.origin
new_ints = interactions_homog(sys)
orig_ints = interactions_homog(origin)
# Transfer interactions from `src` to reshaped `sys`.
#
# Frequently `src` will refer to `sys.origin`, associated with the original
# crystal, which will make symmetry analysis available. In this case, the
# process to set a new interaction is to first modify `sys.origin`, and then to
# call this function.
function transfer_interactions!(sys::System{N}, src::System{N}) where N
new_ints = interactions_homog(sys)

for new_i in 1:natoms(sys.crystal)
i = map_atom_to_crystal(sys.crystal, new_i, origin.crystal)
new_ints[new_i].onsite = orig_ints[i].onsite
# Find `src` interaction either through an atom index or a site index
if is_homogeneous(src)
i = map_atom_to_other_crystal(sys.crystal, new_i, src.crystal)
else
i = map_atom_to_other_system(sys.crystal, new_i, src)
end
src_int = src.interactions_union[i]

# Copy onsite couplings
new_ints[new_i].onsite = src_int.onsite

# Copy pair couplings
new_pc = PairCoupling[]
for pc in orig_ints[i].pair
new_bond = transform_bond(sys.crystal, new_i, origin.crystal, pc.bond)
for pc in src_int.pair
new_bond = transform_bond(sys.crystal, new_i, src.crystal, pc.bond)
isculled = bond_parity(new_bond)
push!(new_pc, PairCoupling(isculled, new_bond, pc.scalar, pc.bilin, pc.biquad, pc.general))
end
Expand All @@ -55,56 +61,33 @@ end


function reshape_supercell_aux(sys::System{N}, new_latsize::NTuple{3, Int}, new_cell_shape::Mat3) where N
is_homogeneous(sys) || error("Cannot reshape system with inhomogeneous interactions.")

# `origin` describes the unit cell of the original system. For sequential
# reshapings, `sys.origin` keeps its original meaning. Make a deep copy so
# that the new system fully owns `origin`, and mutable updates to the
# previous system won't affect this one.
origin = clone_system(isnothing(sys.origin) ? sys : sys.origin)

# If `new_cell_shape == I`, we can effectively restore the unit cell of
# `origin`. Otherwise, we will need to reshape the crystal, map the
# interactions, and keep a reference to the original system.
if new_cell_shape ≈ I
new_cryst = origin.crystal

new_na = natoms(new_cryst)
new_Ns = zeros(Int, new_latsize..., new_na)
new_κs = zeros(Float64, new_latsize..., new_na)
new_gs = zeros(Mat3, new_latsize..., new_na)
new_extfield = zeros(Vec3, new_latsize..., new_na)
new_dipoles = zeros(Vec3, new_latsize..., new_na)
new_coherents = zeros(CVec{N}, new_latsize..., new_na)

new_dipole_buffers = Array{Vec3, 4}[]
new_coherent_buffers = Array{CVec{N}, 4}[]

new_ints = interactions_homog(origin)

new_sys = System(nothing, origin.mode, new_cryst, new_latsize, new_Ns, new_κs, new_gs, new_ints, nothing,
new_extfield, new_dipoles, new_coherents, new_dipole_buffers, new_coherent_buffers, origin.units, copy(sys.rng))
else
new_cryst = reshape_crystal(origin.crystal, new_cell_shape)

new_na = natoms(new_cryst)
new_Ns = zeros(Int, new_latsize..., new_na)
new_κs = zeros(Float64, new_latsize..., new_na)
new_gs = zeros(Mat3, new_latsize..., new_na)
new_extfield = zeros(Vec3, new_latsize..., new_na)
new_dipoles = zeros(Vec3, new_latsize..., new_na)
new_coherents = zeros(CVec{N}, new_latsize..., new_na)

new_dipole_buffers = Array{Vec3, 4}[]
new_coherent_buffers = Array{CVec{N}, 4}[]

new_ints = empty_interactions(origin.mode, new_na, N)

new_sys = System(origin, origin.mode, new_cryst, new_latsize, new_Ns, new_κs, new_gs, new_ints, nothing,
new_extfield, new_dipoles, new_coherents, new_dipole_buffers, new_coherent_buffers, origin.units, copy(sys.rng))

set_interactions_from_origin!(new_sys)
end
# Reshape the crystal
new_cryst = reshape_crystal(orig_crystal(sys), new_cell_shape)
new_na = natoms(new_cryst)

# Allocate data for new system, but with an empty list of interactions
new_Ns = zeros(Int, new_latsize..., new_na)
new_κs = zeros(Float64, new_latsize..., new_na)
new_gs = zeros(Mat3, new_latsize..., new_na)
new_ints = empty_interactions(sys.mode, new_na, N)
new_ewald = nothing
new_extfield = zeros(Vec3, new_latsize..., new_na)
new_dipoles = zeros(Vec3, new_latsize..., new_na)
new_coherents = zeros(CVec{N}, new_latsize..., new_na)
new_dipole_buffers = Array{Vec3, 4}[]
new_coherent_buffers = Array{CVec{N}, 4}[]

# Preserve origin for repeated reshapings. In other words, ensure that
# new_sys.origin === sys.origin
orig_sys = @something sys.origin sys

new_sys = System(orig_sys, sys.mode, new_cryst, new_latsize, new_Ns, new_κs, new_gs, new_ints, new_ewald,
new_extfield, new_dipoles, new_coherents, new_dipole_buffers, new_coherent_buffers, sys.units, copy(sys.rng))

# Transfer interactions. In the case of an inhomogeneous system, the
# interactions in `sys` have detached from `orig`, so we use the latest
# ones.
transfer_interactions!(new_sys, sys)

# Copy per-site quantities
for new_site in eachsite(new_sys)
Expand All @@ -117,14 +100,16 @@ function reshape_supercell_aux(sys::System{N}, new_latsize::NTuple{3, Int}, new_
new_sys.coherents[new_site] = sys.coherents[site]
end

# Restore dipole-dipole interactions if present
# Restore dipole-dipole interactions if present. This involves pre-computing
# an interaction matrix that depends on `new_latsize`.
if !isnothing(sys.ewald)
enable_dipole_dipole!(new_sys)
end

return new_sys
end


# Shape of a possibly reshaped unit cell, given in multiples of the original
# unit cell.
function cell_shape(sys)
Expand Down Expand Up @@ -158,7 +143,8 @@ times in each dimension, specified by the tuple `counts`.
See also [`reshape_supercell`](@ref).
"""
function repeat_periodically(sys::System{N}, counts::NTuple{3,Int}) where N
@assert all(>=(1), counts)
all(>=(1), counts) || error("Require at least one count in each direction.")

# Scale each column by `counts` and reshape
return reshape_supercell_aux(sys, counts .* sys.latsize, cell_shape(sys))
end
2 changes: 1 addition & 1 deletion src/SpinWaveTheory/SpinWaveTheory.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ function SpinWaveTheory(sys::System{N}; energy_ϵ::Float64=1e-8, observables=not

# Reshape into single unit cell. A clone will always be performed, even if
# no reshaping happens.
cellsize_mag = cell_shape(sys) * diagm(SVector(sys.latsize))
cellsize_mag = cell_shape(sys) * diagm(Vec3(sys.latsize))
sys = reshape_supercell_aux(sys, (1,1,1), cellsize_mag)

# Rotate local operators to quantization axis
Expand Down
55 changes: 27 additions & 28 deletions src/Symmetry/Crystal.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ cryst = Crystal(latvecs, positions, 227; setting="1")
See also [`lattice_vectors`](@ref).
"""
struct Crystal
origin :: Union{Nothing, Crystal} # Origin crystal (invariant under `subcrystal` and reshaping)
root :: Union{Nothing, Crystal} # Root crystal (invariant under `subcrystal` and reshaping)
latvecs :: Mat3 # Lattice vectors as columns
prim_latvecs :: Mat3 # Primitive lattice vectors
recipvecs :: Mat3 # Reciprocal lattice vectors (conventional)
Expand Down Expand Up @@ -408,50 +408,49 @@ inputs to [`reshape_supercell`](@ref).
```
"""
function primitive_cell_shape(cryst::Crystal)
origin = @something cryst.origin cryst
if isnothing(origin.prim_latvecs)
root = @something cryst.root cryst
if isnothing(root.prim_latvecs)
error("Primitive lattice vectors not available.")
end
return origin.latvecs \ origin.prim_latvecs
return root.latvecs \ root.prim_latvecs
end


function check_latvecs_commensurate(origin, new_latvecs)
# Work in multiples of the primitive lattice vectors, if they exist.
prim_latvecs = @something origin.prim_latvecs origin.latvecs
cell_shape = prim_latvecs \ new_latvecs
function check_shape_commensurate(cryst, shape)
root = @something cryst.root cryst
if !isnothing(root.prim_latvecs)
shape = root.prim_latvecs \ root.latvecs * shape
end

# Each of these components must be integer
if !(round.(cell_shape) ≈ cell_shape)
if issomething(origin.prim_latvecs)
error("Cell shape must be integer multiples of primitive lattice vectors. Calculated $cell_shape.")
if !(round.(shape) ≈ shape)
if !isnothing(root.prim_latvecs)
error("Cell shape must be integer multiples of primitive lattice vectors. Calculated $shape.")
else
error("Cell shape must be a 3×3 matrix of integers. Received $cell_shape.")
error("Cell shape must be a 3×3 matrix of integers. Received $shape.")
end
end
end

function reshape_crystal(cryst::Crystal, new_cell_shape::Mat3)
# Find the original crystal (unreshaped, prior to subcrystal calls, etc.).
# The field `origin.latvecs` gives meaning to the "fractional" coordinate
# system of `new_cell_shape`. Note that `cryst` may be formed as a
# subcrystal of `origin`. For reshaping purposes, therefore, we must use
# `cryst.positions` and not `origin.positions`, etc.
origin = @something cryst.origin cryst
function reshape_crystal(cryst::Crystal, new_shape::Mat3)
# Check that desired lattice vectors are commensurate with root cell
check_shape_commensurate(cryst, new_shape)

# The `root.latvecs` defines the fractional coordinate system, but note that
# `cryst` may be formed as a subcrystal of `root`. For reshaping purposes,
# therefore, we must use `cryst.positions` and not `root.positions`, etc.
root = @something cryst.root cryst

# Lattice vectors of the new unit cell in global coordinates
new_latvecs = origin.latvecs * new_cell_shape
new_latvecs = root.latvecs * new_shape

# Return this crystal if possible
new_latvecs ≈ cryst.latvecs && return cryst

# Symmetry precision needs to be rescaled for the new unit cell. Ideally we
# would have three separate rescalings (one per lattice vector), but we're
# forced to pick just one. Scale according to volume change.
new_symprec = origin.symprec / cbrt(abs(det(new_cell_shape)))

# Check that desired lattice vectors are commensurate with origin cell
check_latvecs_commensurate(origin, new_latvecs)
new_symprec = root.symprec / cbrt(abs(det(new_shape)))

# This matrix defines a mapping from fractional coordinates `x` in `cryst`
# to fractional coordinates `y` in the new unit cell.
Expand Down Expand Up @@ -503,8 +502,8 @@ function reshape_crystal(cryst::Crystal, new_cell_shape::Mat3)
# lost with the resizing procedure.
new_symops = SymOp[]

return Crystal(origin, new_latvecs, origin.prim_latvecs, new_recipvecs, new_positions, new_types, new_classes,
new_sitesyms, new_symops, origin.spacegroup, new_symprec)
return Crystal(root, new_latvecs, root.prim_latvecs, new_recipvecs, new_positions, new_types, new_classes,
new_sitesyms, new_symops, root.spacegroup, new_symprec)
end


Expand All @@ -531,7 +530,7 @@ function subcrystal(cryst::Crystal, types::Vararg{String, N}) where N
end

function subcrystal(cryst::Crystal, classes::Vararg{Int, N}) where N
origin = @something cryst.origin cryst
root = @something cryst.root cryst
for c in classes
if !(c in cryst.classes)
error("Class '$c' is not present in crystal.")
Expand All @@ -548,7 +547,7 @@ function subcrystal(cryst::Crystal, classes::Vararg{Int, N}) where N
@info "Atoms have been renumbered in subcrystal."
end

ret = Crystal(origin, cryst.latvecs, cryst.prim_latvecs, cryst.recipvecs, new_positions, new_types,
ret = Crystal(root, cryst.latvecs, cryst.prim_latvecs, cryst.recipvecs, new_positions, new_types,
new_classes, new_sitesyms, cryst.symops, cryst.spacegroup, cryst.symprec)
return ret
end
Expand Down
2 changes: 1 addition & 1 deletion src/System/OnsiteCoupling.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,7 @@ function set_onsite_coupling!(sys::System, op, i::Int)
# contains full symmetry information.
if !isnothing(sys.origin)
set_onsite_coupling!(sys.origin, op, i)
set_interactions_from_origin!(sys)
transfer_interactions!(sys, sys.origin)
return
end

Expand Down
2 changes: 1 addition & 1 deletion src/System/PairExchange.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,7 +152,7 @@ function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float
# contains full symmetry information.
if !isnothing(sys.origin)
set_pair_coupling_aux!(sys.origin, scalar, bilin, biquad, tensordec, bond)
set_interactions_from_origin!(sys)
transfer_interactions!(sys, sys.origin)
return
end

Expand Down
25 changes: 17 additions & 8 deletions src/System/System.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,8 @@ function System(crystal::Crystal, latsize::NTuple{3,Int}, infos::Vector{SpinInfo

# The lattice vectors of `crystal` must be conventional (`crystal` cannot be
# reshaped).
if !isnothing(crystal.origin)
@assert crystal.latvecs == crystal.origin.latvecs
if !isnothing(crystal.root)
@assert crystal.latvecs == crystal.root.latvecs
end

na = natoms(crystal)
Expand Down Expand Up @@ -232,8 +232,9 @@ magnetic_moment(sys::System, site) = sys.units.μB * sys.gs[site] * sys.dipoles[
# Total volume of system
volume(sys::System) = cell_volume(sys.crystal) * prod(sys.latsize)

# The original crystal for a system. It is guaranteed to be un-reshaped, and its
# lattice vectors define the "conventional" unit cell.
# The crystal originally used to construct a system. It is guaranteed to be
# un-reshaped, and its lattice vectors define the "conventional" unit cell. It
# may, however, be a subcrystal of `orig_crystal(sys).root`.
orig_crystal(sys) = something(sys.origin, sys).crystal

# Position of a site in fractional coordinates of the original crystal
Expand Down Expand Up @@ -279,12 +280,20 @@ function site_to_atom(sys::System{N}, site) where N
end

# Maps atom `i` in `cryst` to the corresponding atom in `orig_cryst`
function map_atom_to_crystal(cryst, i, orig_cryst)
r = cryst.positions[i]
orig_r = orig_cryst.latvecs \ cryst.latvecs * r
function map_atom_to_other_crystal(cryst::Crystal, i, orig_cryst::Crystal)
global_r = cryst.latvecs * cryst.positions[i]
orig_r = orig_cryst.latvecs \ global_r
return position_to_atom(orig_cryst, orig_r)
end

# Maps atom `i` in `cryst` to the corresponding site in `orig_sys`
function map_atom_to_other_system(cryst::Crystal, i, orig_sys::System)
global_r = cryst.latvecs * cryst.positions[i]
orig_r = orig_crystal(orig_sys).latvecs \ global_r
return position_to_site(orig_sys, orig_r)
end


# Given a `bond` for `cryst`, return a corresponding new bond for the reshaped
# `new_cryst`. The new bond will begin at atom `new_i`.
function transform_bond(new_cryst::Crystal, new_i::Int, cryst::Crystal, bond::Bond)
Expand Down Expand Up @@ -323,7 +332,7 @@ function symmetry_equivalent_bonds(sys::System, bond::Bond)

for new_i in 1:natoms(sys.crystal)
# atom index in original crystal
i = map_atom_to_crystal(sys.crystal, new_i, orig_crystal(sys))
i = map_atom_to_other_crystal(sys.crystal, new_i, orig_crystal(sys))

# loop over symmetry equivalent bonds in original crystal
for bond′ in all_symmetry_related_bonds_for_atom(orig_crystal(sys), i, bond)
Expand Down

0 comments on commit 72b35a6

Please sign in to comment.