Progress

src/EntangledUnits/EntangledUnits.jl

@@ -117,6 +117,11 @@ end
 # Pull out original indices of sites in entangled unit
 sites_in_unit(contraction_info, i) = [inverse_data.site for inverse_data in contraction_info.inverse[i]] 
 
+function original_units(esys::EntangledSystem)
+    (; sys, contraction_info) = esys
+    return [Tuple(sites_in_unit(contraction_info, unit)) for unit in axes(sys.dipoles, 4)]
+end
+
 # List of all pair-wise bonds in a unit. The resulting bonds are to be
 # interpreted in terms of the original crystal.
 function bonds_in_unit(contraction_info, i)
@@ -466,7 +471,12 @@ function expected_dipoles_of_entangled_system!(dipole_buf, esys::EntangledSystem
         end
     end
 end
-expected_dipoles_of_entangled_system!(sys_original, sys_entangled) = expected_dipoles_of_entangled_system!(sys_original.dipoles, sys_entangled)
+
+function sync_dipoles(esys::EntangledSystem)
+    esys.synced && return true
+    expected_dipoles_of_entangled_system!(esys.sys_origin.dipoles, esys)
+    esys.synced = true
+end
 
 
 # TODO: Write this function:

src/EntangledUnits/Types.jl

@@ -28,47 +28,116 @@ end
 #     Ns_unit           :: Vector{Vector{Int64}}   # This eliminates need to carry original system around in many places
 # end
 
-struct EntangledSystem
-    sys               :: System
-    sys_origin        :: System
-    contraction_info  :: CrystalContractionInfo
-    Ns_unit           :: Vector{Vector{Int64}}   # This eliminates need to carry original system around in many places -- now that system is there, possibly eliminate
+mutable struct EntangledSystem
+    const sys               :: System
+    const sys_origin        :: System
+    synced                  :: Bool
+    const contraction_info  :: CrystalContractionInfo
+    const Ns_unit           :: Vector{Vector{Int64}}   # This eliminates need to carry original system around in many places -- now that system is there, possibly eliminate
 end
 
 function EntangledSystem(sys, units)
     (; sys_entangled, contraction_info, Ns_unit) = entangle_system(sys, units)
     sys_origin = clone_system(sys)
-    EntangledSystem(sys_entangled, sys_origin, contraction_info, Ns_unit)
+    EntangledSystem(sys_entangled, sys_origin, false, contraction_info, Ns_unit)
 end
 
 
 ################################################################################
 # Aliasing and field forwarding
 ################################################################################
 
-# 
+# Functions to access System fields of EntangledSystem
 randomize_spins!(esys::EntangledSystem; kwargs...) = randomize_spins!(esys.sys; kwargs...)
 minimize_energy!(esys::EntangledSystem; kwargs...) = minimize_energy!(esys.sys; kwargs...)
 energy(esys::EntangledSystem; kwargs...) = energy(esys.sys; kwargs...)
 set_coherent!(esys::EntangledSystem, coherent, site; kwargs...) = set_coherent!(esys.sys, coherent, site; kwargs...)
-eachsite(esys::EntangledSystem) = eachsite(esys.sys)
+eachsite(esys::EntangledSystem) = eachsite(esys.sys) # Not sure that we want this
 
-# Temporary until better observable approach for classical
+# Functions acting on original System of an EntangledSystem 
+set_dipole!(esys::EntangledSystem, dipole, site; kwargs...) = error("Setting dipoles of an EntangledSystem not well defined.") # Could replicate behavior of normal SU(N) system
+function magnetic_moment(esys::EntangledSystem, site; kwargs...) 
+    magnetic_moment(esys.sys_origin, site; kwargs...)
+end
 function plot_spins(esys::EntangledSystem; kwargs...)
-    expected_dipoles_of_entangled_system!(esys.sys_origin.dipoles, esys)
+    sync_dipoles(esys)
     plot_spins(esys.sys_origin; kwargs...)
 end
 
+# Functions acting on both systems
+function reshape_supercell(esys::EntangledSystem, shape)
+    (; sys, sys_origin, contraction_info) = esys
+    new_sys_origin = reshape_supercell(sys_origin, shape)
+
+    new_crystal = sys_origin_new.crystal
+    new_na = natoms(new_crystal)
+
+    units = original_units(esys)
+    new_units = []
+
+    new_atoms = 1:new_na
+    new_units = []
+    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) # Just work with first unit cell
+
+        # Find corresponding original atom number
+        site = position_to_site(sys_origin, position(new_sys_origin, new_site))
+        atom = site[4]
+
+        # Find the unit to which this original atom belongs
+        unit = findfirst(unit -> 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(new_sys_origin, position) for position in unit_positions]
+        new_unit = Int64[]
+        for new_site in new_unit_sites
+            i, j, k, a = new_site
+            if !(i == j == k == 1)
+                error("Specified reshaping incomptable with specified entangled units. (Unit split between unit cells.)")
+            end
+            push!(new_unit, a)
+        end
+        push!(new_units, Tuple(new_unit))
+
+        idcs = findall(atom -> atom in new_unit, new_atoms)
+        deleteat!(new_atoms, idcs)
+    end
+
+    # Construct new ContractionInfo for reshaped system
+    _, contraction_info = contract_crystal(new_sys_origin.crystal, new_units)
+
+    # Reshape the entangled system as well
+    new_sys = reshape_supercell(sys, shape)
+
+    # new_esys = EntangledSystem()
+
+end
+
+function repeat_periodically(esys, counts)
+    sync_dipoles(esys)
+    sys_new = repeat_periodically(esys.sys, counts)
+    sys_origin_new = repeat_periodically(esys.sys_origin, counts)
+    return EntangledSystem(sys_new, sys_origin_new, true, esys.contraction_info, esys.Ns_unit)
+end
+
+
 # Forward field requests to internal system
 function Base.getproperty(value::EntangledSystem, name::Symbol)
-    if name ∉ [:sys, :sys_origin, :contraction_info, :Ns_unit] 
+    if name == :coherents
         return getfield(value.sys, name)
+    elseif name == :dipoles
+        return getfield(value.sys_origin, name)
     end
     return getfield(value, name)
 end
 function Base.setproperty!(value::EntangledSystem, name::Symbol, x)
-    if name ∉ [:sys, :sys_origin, :contraction_info, :Ns_unit] 
+    if name == :coherents 
         return setfield!(value.sys, name, convert(fieldtype(System, name), x))
+    elseif name == :dipoles
+        error("Cannot set `dipoles` of EntangledSystem directly.")
     end
     return setfield!(value, name, convert(fieldtype(EntangledSystem, name), x))
 end