Entangled units updates (#324)

* Fixes for external fields and more tests
* Additional function aliasing
* Plumbing for setting entangled unit states from dipoles

src/EntangledUnits/EntangledReshaping.jl

@@ -85,5 +85,19 @@ function repeat_periodically(esys, counts)
     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 resize_supercell(esys::EntangledSystem, dims::NTuple{3,Int})
+    (; sys, sys_origin, contraction_info) = esys
+
+    # Resize both entangled and original system periodically
+    sys_new = reshape_supercell(sys, diagm(Vec3(dims)))
+    sys_origin_new = reshape_supercell(sys_origin, diagm(Vec3(dims)))
+
+    # 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

src/EntangledUnits/EntangledSpinWaveTheory.jl

@@ -255,6 +255,53 @@ function excitations!(T, tmp, swt::EntangledSpinWaveTheory, q)
     end
 end
 
+# No changes
+function excitations(swt::EntangledSpinWaveTheory, q)
+    L = nbands(swt)
+    T = zeros(ComplexF64, 2L, 2L)
+    H = zeros(ComplexF64, 2L, 2L)
+    energies = excitations!(T, copy(H), swt, q)
+    return (energies, T)
+end
+
+# No changes
+function dispersion(swt::EntangledSpinWaveTheory, qpts)
+    L = nbands(swt)
+    qpts = convert(AbstractQPoints, qpts)
+    disp = zeros(L, length(qpts.qs))
+    for (iq, q) in enumerate(qpts.qs)
+        view(disp, :, iq) .= view(excitations(swt, q)[1], 1:L)
+    end
+    return reshape(disp, L, size(qpts.qs)...)
+end
+
+# No changes
+function energy_per_site_lswt_correction(swt::EntangledSpinWaveTheory; opts...)
+    any(in(keys(opts)), (:rtol, :atol, :maxevals)) || error("Must specify one of `rtol`, `atol`, or `maxevals` to control momentum-space integration.")
+
+    (; sys) = swt
+    Natoms = natoms(sys.crystal)
+    L = nbands(swt)
+    H = zeros(ComplexF64, 2L, 2L)
+    V = zeros(ComplexF64, 2L, 2L)
+
+    # The uniform correction to the classical energy (trace of the (1,1)-block
+    # of the spin-wave Hamiltonian)
+    dynamical_matrix!(H, swt, zero(Vec3))
+    δE₁ = -real(tr(view(H, 1:L, 1:L))) / Natoms
+
+    # Integrate zero-point energy over the first Brillouin zone 𝐪 ∈ [0, 1]³ for
+    # magnetic cell in reshaped RLU
+    δE₂ = hcubature((0,0,0), (1,1,1); opts...) do q_reshaped
+        dynamical_matrix!(H, swt, q_reshaped)
+        ωs = bogoliubov!(V, H)
+        return sum(view(ωs, 1:L)) / 2Natoms
+    end
+
+    # Error bars in δE₂[2] are discarded
+    return δE₁ + δE₂[1]
+end
+
 # Only change is no Ewald
 function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::EntangledSpinWaveTheory, q_reshaped::Vec3)
     (; sys) = swt

src/EntangledUnits/EntangledUnits.jl

@@ -244,6 +244,10 @@ function entangle_system(sys::System{M}, units) where M
     # Construct contracted crystal
     contracted_crystal, contraction_info = contract_crystal(sys.crystal, units)
 
+    # Make sure we have a uniform external field
+    @assert allequal(@view sys.extfield[:,:,:,:]) "`EntangledSystems` requires a uniform applied field." 
+    B = sys.extfield[1,1,1,1]
+
     # 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)
@@ -254,6 +258,9 @@ function entangle_system(sys::System{M}, units) where M
     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)
 
+    # Transfer rng from origin system to entangled system
+    copy!(sys_entangled.rng, sys.rng)
+
     # TODO: Extend to inhomogenous systems
     # For each contracted site, scan original interactions and reconstruct as necessary.
     new_pair_data = Tuple{Bond, Matrix{ComplexF64}}[]
@@ -267,9 +274,14 @@ function entangle_system(sys::System{M}, units) where M
         # Pair interactions that become within-unit interactions
         original_interactions = sys.interactions_union[relevant_sites] 
         for (site, interaction) in zip(relevant_sites, original_interactions)
+            # Onsite anisotropy portion
             onsite_original = interaction.onsite
             unit_index = contraction_info.forward[site][2]
             unit_operator += local_op_to_product_space(onsite_original, unit_index, Ns)
+
+            # Zeeman portion
+            S = [local_op_to_product_space(S, unit_index, Ns) for S in spin_matrices((Ns[unit_index]-1)/2)]
+            unit_operator += Hermitian((sys.gs[1, 1, 1, site] * B)' * S)
         end
 
         # Sort all PairCouplings in couplings that will be within a unit and couplings that will be between units

src/EntangledUnits/TypesAndAliasing.jl

@@ -26,6 +26,10 @@ struct CrystalContractionInfo
     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
 
+function Base.copy(cci::CrystalContractionInfo)
+    return CrystalContractionInfo(copy(cci.forward), copy(cci.inverse))
+end
+
 
 ################################################################################
 # System 
@@ -65,25 +69,23 @@ function EntangledSystem(sys::System{N}, units) where {N}
     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
+    # Generate observable field. This observable field has as many entries 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)
+    (; 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
@@ -97,7 +99,7 @@ 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)))")
+    print(io, "EntangledSystem($(mode_to_str(esys.sys)), $(supercell_to_str(esys.sys_origin.dims, esys.sys_origin.crystal)), $(energy_to_str(esys.sys)))")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", esys::EntangledSystem)
@@ -116,6 +118,16 @@ 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 clone_system(esys::EntangledSystem)
+    sys = clone_system(esys.sys)
+    sys_origin = clone_system(esys.sys_origin)
+    contraction_info = copy(esys.contraction_info)
+    dipole_operators = copy(esys.dipole_operators)
+    source_idcs = copy(esys.source_idcs)
+
+    return EntangledSystem(sys, sys_origin, contraction_info, dipole_operators, source_idcs)
+end
+
 function set_field!(esys::EntangledSystem, B)
     (; sys, sys_origin, dipole_operators, source_idcs) = esys 
     B_old = sys_origin.extfield[1,1,1,1] 
@@ -127,16 +139,54 @@ function set_field!(esys::EntangledSystem, B)
         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)
+        sys.interactions_union[unit].onsite += Hermitian(ΔB' * S)
     end
 end
 
+# TODO: Actually, we could give a well-defined meaning to this procedure. Implement this.
 function set_field_at!(::EntangledSystem, _, _)
-    error("`EntangledSystem`s do not support inhomogenous external fields. Use `set_field!(sys, B).")
+    error("`EntangledSystem`s do not support inhomogenous external fields. Use `set_field!(sys, B) to set a uniform field.")
 end
 
 
-set_dipole!(esys::EntangledSystem, dipole, site; kwargs...) = error("Setting dipoles of an EntangledSystem not well defined.") 
+function set_dipole!(::EntangledSystem, dipole, site)
+    error("`set_dipole!` operation for `EntangledSystem` not well defined. Consider using `set_coherent!` to set the state of each entangled unit.")
+end
+
+# Find the unique coherent state corresponding to a set of fully-polarized
+# dipoles on each site inside a specified entangled unit.
+function coherent_state_from_dipoles(esys::EntangledSystem, dipoles, unit)
+    (; sys_origin, contraction_info) = esys
+
+    # Find the atom indices (of original system) that lie in the specified unit
+    # (of contracted system).
+    atoms = [id.site for id in contraction_info.inverse[unit]]
+
+    # Test that the number of specified dipoles is equal to the number of atoms
+    # inside the entangled unit.
+    @assert length(dipoles) == length(atoms) "Invalid number of dipoles for specified unit."
+
+    # Retrieve the dimensions of the local Hilbert spaces corresponding to those
+    # atoms.
+    Ns = Ns_in_units(sys_origin, contraction_info)[unit]
+
+    # Generate a list of coherent states corresponding to given dipoles _in each
+    # local Hilbert space_.
+    coherents = []
+    for (dipole, N) in zip(dipoles, Ns)
+        # Get the spin matrices in the appropriate representation for the site.
+        S = spin_matrices((N-1)/2)
+
+        # Find a local coherent state representation of the dipole
+        coherent = eigvecs(S' * dipole)[:,N] # Retrieve highest-weight eigenvector
+        push!(coherents, coherent)
+    end
+
+    # Return the tensor product of each of these local coherent states to get
+    # the coherent state for the entangled unit.
+    return kron(coherents...)
+end
+
 
 # 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